2021.10.18 -10.28基于mRNA表达表格做差异基因

2021-10-21 本文已影响0人千容安

rm(list=ls())
rm(list=ls())用于清除所有的变量，ls() 会返回当前环境里所有的objects的名字，list是rm() 里边的options之一，意思是“想要删除的变量的名字是。。”。rm(list=ls()) 意思是我要删除变量，删除变量的名字是ls()。

读数据的命令出错

BiocManager::install("data.table")
library(data.table)
用这个包读数据很快，而且不会有很多格式的问题

读数据的另一个命令
data <- read.table(choose.files(),header = TRUE,sep="\t")
data <- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F) 将数据储存为txt格式后再读
每次读文件之前，得告诉R你文件所在的位置，setwd("")
setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop")

看数据的话，是fix(data)，或者head(data[1:10,1:10])

fix(data)

data[1:10,1:10]

data<-fread("RNA-seq - 1.txt")
输入数据：data<-fread("文件名.txt")
将xlsx数据另存为txt格式。

火山图

data$color <- ifelse(data$padj<0.05 & abs(data$log2FoldChange)>= 1,ifelse(data$log2FoldChange > 1,'red','blue'),'gray')
# abs：绝对值，对于有生物学重复的样本，要求padj﹤0.05
color <- c(red = "red",gray = "gray",blue = "blue")

p <- ggplot(data, aes(log2FoldChange, -log10(padj), col = color)) +
 geom_point() +
  theme_bw() +
 scale_color_manual(values = color) +
labs(x="log2 (fold change)",y="-log10 (q-value)") +
  geom_hline(yintercept = -log10(0.05), lty=4,col="grey",lwd=0.6) +
  geom_vline(xintercept = c(-1, 1), lty=4,col="grey",lwd=0.6) +
  theme(legend.position = "none",
       panel.grid=element_blank(),
       axis.title = element_text(size = 16),
      axis.text = element_text(size = 14))
p

控制坐标轴范围：
p+scale_x_continuous(limits = c(-3,3))+scale_y_continuous(limits = c(0,5)) xlim
选最大值作为xlim的上下边界

x_lim <- max(logFC,-logFC)
y_lim<-max(-1*log10(padj))

确定x轴最大值：

colnames(data)
 [1] "gene_id"        "S10"            "S2NAC"         
 [4] "S7"             "S9"             "M3"            
 [7] "M7"             "M8"             "M9"            
[10] "baseMeanA"      "baseMeanB"      "baseMean"      
[13] "log2FoldChange" "lfcSE"          "stat"          
[16] "pvalue"         "padj"           "type"          
[19] "GO_BP"          "GO_MF"          "GO_CC"         
[22] "KEGG_PATHWAY"   "color"         
> max(abs(data[,13]))
[1] 23.22173

绘制火山图

library(ggplot2)
library(RColorBrewer)
theme_set(theme_bw())
p <- ggplot(data,aes(log2FoldChange,-1*log10(padj),
                   color = sig))+geom_point(aes(color=sig))+
   xlim(-23.22,22.23) +  labs(x="log2 Fold Change",y="-log10(pvalue)")
p <- p + scale_color_manual(values =c("#0072B5","grey","#BC3C28"))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(0),linetype=4)
p <- p +theme(axis.line = element_line(size=0))+ylim(0,y_lim)
p <- p +theme(axis.text=element_text(size=20),axis.title=element_text(size=20))
p

PCA图

根据基因表达值进行样本间的PCA分析，确定样本在PCA图中的位置，R语言中能够执行PCA分析的方法有很多，不过它们的算法都是统一的，随便使用任何一个R包就可以

# 读取数据
data <- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
# 提取第一列   (不用这一条也行）
symbol<-as.matrix(gene[,1])
# 选2-11列（自动将第一列当列名）
exp = as.matrix(data[,c(1:8)])
# 要用 log 转化后的基因表达值，降低不同基因表达水平数量级相差过大的问题
exp<-log2(exp+1)
#我们使用 FactoMineR 包中的方法，实现 PCA 分析和聚类添加
library(FactoMineR)
#样本中基因表达值的 PCA 分析
exp.pca <- PCA(exp, ncp = 2, scale.unit = TRUE, graph = FALSE)
plot(exp.pca)

读取数据遇到的一些问题：
当我查看数据时，前面都是没有内容的，还是RNA-seq的数据

代码如下：

解决如下：在运行标红的代码前，只要截取前8列的表达谱数据即可，因为做PCA图需要的数据只有前8列。

上一步获得了PCA分析结果，并观察到明显的组间差异。接下来使用ggplot2包绘制一个好看的PCA图。
需要事先在上述PCA分析结果中将关键信息提取出，例如样本点在PCA图中的位置信息，以及PCA轴的贡献度等。

#提取 PCA 中的样本位置，即x轴（pca1轴）和y轴（pca2轴）
pca_sample <- data.frame(exp.pca$ind$coord[ ,1:2])
head(pca_sample)
           Dim.1      Dim.2
Ciart   8.186812 -0.6273470
Fabp7   7.368796  0.5651271
Fmo2    4.965070 -1.0870700
Hif3a   7.467280 -0.7283460
Kirrel2 7.125757 -0.4987601
Plin4   6.743315 -0.8912442
#提取 PCA 前两轴的贡献度
pca_eig1 <- round(exp.pca$eig[1,2], 2)  
pca_eig2 <- round(exp.pca$eig[2,2],2 )

读取并合并样本分组信息
group <- read.delim('pca.xlsx', row.names = 1, sep = '\t', check.names = FALSE)
group <- group[rownames(pca_sample), ]
pca_sample <- cbind(pca_sample, group)

pca_sample$samples <- rownames(pca_sample)
head(pca_sample)  #作图数据中包含了样本坐标和分组信息

#ggplot2 绘制二维散点图
library(ggplot2)

p <- ggplot(data = pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group), size = 3) +  #根据样本坐标绘制二维散点图
scale_color_manual(values = c('orange', 'purple')) +  #自定义颜色
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
    legend.key = element_rect(fill = 'transparent')) +  #去除背景和网格线
labs(x =  paste('PCA1:', pca_eig1, '%'), y = paste('PCA2:', pca_eig2, '%'), color = '')  #将 PCA 轴贡献度添加到坐标轴标题中
p

若出现Error: unexpected symbol in "\u200b" ，将#注释去掉。

重新作图：

df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1)
iris.pca<-prcomp(t(df[,1:8]))
summary(iris.pca)
plot(iris.pca$x,col=rep(c("red","blue"),c(4,4)),pch=19)

PCA

用fread读数据，data.table=F，这样也是数值型

library("data.table")
da<-fread("RNA-seq - 1.txt",data.table=F)

getwd()
"C:/Users/Administrator.DESKTOP-4UQ3Q0K/Desktop"
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
iris.pca<-prcomp(df[,1:8])  #计算主成分
iris_pcs <-data.frame(iris.pca$x)
head(iris_pcs)
              PC1       PC2       PC3       PC4         PC5
Ciart   -4934.580 5217.2257 1760.1106 -248.6604  262.260206
Fabp7   -4472.953 1544.6017 -349.0215  595.1240   -8.948593
Fmo2    -6893.075  676.7536  266.6638 -213.0876  132.676879
Hif3a   -5814.847 3522.1495 1007.5449 -571.4270  489.318537
Kirrel2 -6008.083 2239.0950  625.8803   17.2822 -132.076437
Plin4   -6351.092 2396.1754  794.4072 -473.1440  365.604000
               PC6        PC7         PC8
Ciart    161.56826  103.28029   -4.916623
Fabp7   -417.51164 -230.57635  -16.571842
Fmo2      19.73212   70.27110  -46.352493
Hif3a    169.11638  163.44406 -118.575331
Kirrel2  -34.05482   36.89955   61.624917
Plin4    172.99939  251.82729   61.521027


summary(iris.pca)
Importance of components:
                             PC1       PC2       PC3   PC4 PC5   PC6
Standard deviation     1.640e+05 1.255e+03 639.23739 338.9 135 77.25
Proportion of Variance 9.999e-01 6.000e-05   0.00002   0.0   0  0.00
Cumulative Proportion  9.999e-01 1.000e+00   0.99999   1.0   1  1.00
                         PC7   PC8
Standard deviation     62.57 29.65
Proportion of Variance  0.00  0.00
Cumulative Proportion   1.00  1.00
# 做碎石图
library(ggplot2)
library(factoextra)
fviz_eig(iris.pca)
fviz_eig(iris.pca,addlabels = T)

碎石图

plot(exp.pca$x,cex = 2,main = "PCA analysis", 
     col = c(rep("red",3),rep("blue",3)),
     pch = c(rep(16,3),rep(17,3)))
plot(exp.pca)

PCA analysis

一些报错的解决办法
data.pca <- prcomp(data)
Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
解决办法：
data<-as.numeric(as.matrix(data))
这个错误是在使用read.table()函数读取文件遇到的，经过摸索，发现read.table()读取的数据类型为data.frame, 通过将y的类型由data.frame转换为numeric可解决报错。
一些命令含义
将基因表达值矩阵作个转置，使行为样本，列为基因

dim(df) # 查看行和列数
[1] 603  21
tmpp<-apply(tmpp,2,as.numeric)  # 2表示按列转数值类型
dim(tmpp)
[1]   8 603

按教程画一遍：

library(ggplot2)
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1)
head(df)
df_pca <- prcomp(df)
# Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
mode(df)
# "list"
df<-apply(df,2,as.numeric)
mode(df)
# "numeric"
df_pca<-prcomp(t(df[,1:8]))
df_pcs <-data.frame(df_pca$x, feature = rep(c("control","case"),c(4,4)))   
# 用rep构建标签

ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+ geom_point()

ggplot

去掉背景及网格线

ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+ 
+     geom_point()+ 
+     theme_bw() +
+     theme(panel.border=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.line= element_line(colour = "black"))

image.png

添加PC1 PC2的百分比

percentage<-paste(colnames(df_pcs),"(", paste(as.character(percentage), "%", ")", sep=""))
ggplot(df_pcs,aes(x=PC1,y=PC2,color=feature))+
+     geom_point(size=3)
+     xlab(percentage[1]) 
+     ylab(percentage[2])

image.png

使用ggplot2包绘制PCA图

library(ggplot2)
df<- read.table("RNA-seq - 1.txt",header = T,row.names = 1,sep="\t",check.names = F)
df[,1:8]<-log2(df[1:8]+1) 
head(df)
df<-apply(df,2,as.numeric)
mode(df)
df_pca<-prcomp(t(df[,1:8]))
data.pca<-prcomp(t(df[,1:8]))
summary(data.pca)
pca.scores <- as.data.frame(data.pca$scores)
head(pca.scores)
library(ggplot2)
ggplot(pca.scores,aes(PC1,PC2,col=rownames(pca.scores))) 
+     geom_point(size=3) 
+     geom_text(aes(label=rownames(pca.scores)),vjust = "outward") 
+     geom_hline(yintercept = 0,lty=2,col="red") 
+     geom_vline(xintercept = 0,lty=2,col="blue",lwd=1) 
+     theme_bw() + theme(legend.position = "none") 
+     theme(plot.title = element_text(hjust = 0.5)) 
+     labs(x="PCA_1",y="PCA_2",title = "PCA analysis")

2021.10.18 -10.28基于mRNA表达表格做差异基因

火山图

PCA图

猜你喜欢

热点阅读