表达矩阵log后进行差异分析

关于limma包差异分析结果的logFC解释

http://www.bio-info-trainee.com/1209.html

首先，我们要明白，limma接受的输入参数就是一个表达矩阵，而且是log后的表达矩阵（以2为底）。

那么最后计算得到的logFC这一列的值，其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值，那么就会有正负之分，代表了case相当于control组来说，该基因是上调还是下调。

我之前总是有疑问，明明是case一组的平均表达量和control一组的平均表达量差值呀，跟log foldchange没有什么关系呀。

后来，我终于想通了，因为我们输入的是log后的表达矩阵，那么case一组的平均表达量和control一组的平均表达量都是log了的，那么它们的差值其实就是log的foldchange

首先，我们要理解foldchange的意义，如果case是平均表达量是8，control是2，那么foldchange就是4，logFC就是2咯

那么在limma包里面，输入的时候case的平均表达量被log后是3，control是1，那么差值是2，就是说logFC就是2。

这不是巧合，只是一个很简单的数学公式log(x)-log(y)=log(x/y)

即FC（foldchange）=x/y
进行log处理后再做一次

> ###log处理
> exp_log <- log2(exp+1)
> ###挑出数据
> exp1 <-exp_log[,7:11]#即三组control,两组EGF 12
> ###查看数据
> dim(exp1)
[1] 31099     5
> dat1 <- as.data.frame(t(exp1))###转置列为基因名，行为样本
> dat1[1:4,1:4]
         1367452_at 1367453_at 1367454_at 1367455_at
GSM75272   14.07837   13.20202   14.19491   14.29942
GSM75273   14.11096   13.09994   13.97388   14.36622
GSM75274   14.27705   12.81412   14.30077   14.30550
GSM75275   13.82473   13.70316   14.69665   14.11910
> dim(exp1)
[1] 31099     5
> dim(dat1)
[1]     5 31099
> dat1 <- cbind(dat1,group_list_1)###将分组信息加入
> dim(dat1)
[1]     5 31100
> library(dplyr)
> ###PCA的统一操作走起
> library(FactoMineR)#画主成分分析图需要加载这两个包
Warning message:
程辑包‘FactoMineR’是用R版本3.6.1 来建造的 
> library(factoextra)
载入需要的程辑包：ggplot2
Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
> ###进行主成分分析
> dat1.pca <- PCA(dat1[,-ncol(dat1)], graph = FALSE)
> ###画PCA图
> fviz_pca_ind(dat1.pca,
+              geom.ind = "point", # show points only (nbut not "text")
+              col.ind = dat1$group_list_1, # color by groups
+              palette = c("#00AFBB", "#E7B800"),##https://www.114la.com/other/rgb.htm
+              addEllipses = TRUE, # Concentration ellipses
+              legend.title = "Groups"
+ )
Too few points to calculate an ellipse
Too few points to calculate an ellipse
> ?fviz_pca_ind
> ggsave('all_samples_PCA.png')
Saving 9.03 x 5.6 in image
Too few points to calculate an ellipse
Too few points to calculate an ellipse

PCA图

单个基因差异分析

> table(group_list_1)
group_list_1
control  EGF 12 
      3       2 
> exp1[1:4,1:4]
           GSM75272 GSM75273 GSM75274 GSM75275
1367452_at    14.08    14.11    14.28    13.82
1367453_at    13.20    13.10    12.81    13.70
1367454_at    14.19    13.97    14.30    14.70
1367455_at    14.30    14.37    14.31    14.12
> boxplot(exp1[1,]~group_list_1)

image.png

> boxplot(exp1[1,]~group_list_1)
> bp=function(g){         #定义一个函数g，函数为{}里的内容
+   library(ggpubr)
+   df=data.frame(gene=g,group=group_list_1)
+   p <- ggboxplot(df, x = "group", y = "gene",
+                  color = "group", palette = "jco",
+                  add = "jitter")
+   #  Add p-value
+   p + stat_compare_means(label.y = 8)
+ }
> bp(exp1[1,]) ## 调用上面定义好的函数，避免同样的绘图代码重复多次敲。
Loading required package: magrittr

image.png

> #差异分析，用limma包来做
> #需要表达矩阵和group_list，其他都不要动
> library(limma)
> design=model.matrix(~factor(group_list_1))
> fit=lmFit(exp1,design)
> fit=eBayes(fit)
> dim(fit)
[1] 31099     2
> #差异基因排名
> deg=topTable(fit,coef=2,number = Inf)
> head(deg)
              logFC AveExpr      t   P.Value adj.P.Val     B
1370019_at   -5.522  11.145 -34.63 1.736e-06   0.01261 5.299
1376711_at   -5.240  13.389 -32.54 2.274e-06   0.01261 5.170
1387958_at    6.431  10.722  28.22 4.221e-06   0.01261 4.835
1370387_at    5.767   9.032  27.89 4.438e-06   0.01261 4.805
1370982_at   -4.201   6.493 -27.84 4.477e-06   0.01261 4.800
1388054_a_at  4.544  10.726  27.78 4.521e-06   0.01261 4.794
> bp(exp1[rownames(deg)[1],])

image.png

> #为deg数据框添加几列
> #1.加probe_id列，把行名变成一列
> library(dplyr)
> deg <- mutate(deg,probe_id=rownames(deg))
> #tibble::rownames_to_column()
> head(deg)
   logFC AveExpr      t   P.Value adj.P.Val     B     probe_id
1 -5.522  11.145 -34.63 1.736e-06   0.01261 5.299   1370019_at
2 -5.240  13.389 -32.54 2.274e-06   0.01261 5.170   1376711_at
3  6.431  10.722  28.22 4.221e-06   0.01261 4.835   1387958_at
4  5.767   9.032  27.89 4.438e-06   0.01261 4.805   1370387_at
5 -4.201   6.493 -27.84 4.477e-06   0.01261 4.800   1370982_at
6  4.544  10.726  27.78 4.521e-06   0.01261 4.794 1388054_a_at

> #2.加symbol列，火山图要用
#id转换，查找芯片平台对应的包
http://www.bio-info-trainee.com/1399.html
 "GPL1355"###注释包为rat2302.db
[1] "GPL1355"
> library(rat2302.db)
> ls("package:rat2302.db")
 [1] "rat2302"              "rat2302.db"           "rat2302_dbconn"      
 [4] "rat2302_dbfile"       "rat2302_dbInfo"       "rat2302_dbschema"    
 [7] "rat2302ACCNUM"        "rat2302ALIAS2PROBE"   "rat2302CHR"          
[10] "rat2302CHRLENGTHS"    "rat2302CHRLOC"        "rat2302CHRLOCEND"    
[13] "rat2302ENSEMBL"       "rat2302ENSEMBL2PROBE" "rat2302ENTREZID"     
[16] "rat2302ENZYME"        "rat2302ENZYME2PROBE"  "rat2302GENENAME"     
[19] "rat2302GO"            "rat2302GO2ALLPROBES"  "rat2302GO2PROBE"     
[22] "rat2302MAPCOUNTS"     "rat2302ORGANISM"      "rat2302ORGPKG"       
[25] "rat2302PATH"          "rat2302PATH2PROBE"    "rat2302PFAM"         
[28] "rat2302PMID"          "rat2302PMID2PROBE"    "rat2302PROSITE"      
[31] "rat2302REFSEQ"        "rat2302SYMBOL"        "rat2302UNIGENE"      
[34] "rat2302UNIPROT"      
> ids <- toTable(rat2302SYMBOL)
> head(ids)
    probe_id symbol
1 1367452_at  Sumo2
2 1367453_at  Cdc37
3 1367454_at  Copb2
4 1367455_at    Vcp
5 1367457_at  Becn1
6 1367458_at Lypla2
#3.加change列：上调或下调，火山图要用
logFC_t <- 1 #不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
> change=ifelse(deg$P.Value>0.01,'stable', 
+               ifelse( deg$logFC >logFC_t,'up', 
+                       ifelse( deg$logFC < -logFC_t,'down','stable') )
+ )
> deg <- mutate(deg,change)
> table(deg$change)

  down stable     up 
  1254  18732    932

火山图

plot(deg$logFC,-log10(deg$P.Value))

image.png

> library(ggpubr)
> ggscatter(dat1, x = "logFC", y = "v",size=0.5,color = "change")

image.png

> ggscatter(dat1, x = "logFC", y = "v", color = "change",size = 0.5,
+           label = "symbol", repel = T,
+           label.select = dat1$symbol[1:30] ,
+           #label.select = c('CD36','DUSP6'), #挑选一些基因在图中显示出来
+           palette = c("#00AFBB", "#999999", "#FC4E07") )

image.png

绘制热图

> library(pheatmap)
> table(deg$change)

  down stable     up 
  1254  18732    932 
> table(x=='up')

FALSE  TRUE 
19986   932 
> table(x[x=='up'])

 up 
932 
> cg <- names(x[x=='up'])

归一化画热图

> n=t(scale(t(exp1[cg,])))
> thr=2
> n[n>thr]=thr
> n[n< -thr]= -thr
> n[1:4,1:4]
             GSM75272 GSM75273 GSM75274 GSM75275
1387958_at    -0.6468  -0.7954  -0.7451    1.066
1370387_at    -0.6985  -0.8093  -0.6800    1.122
1388054_a_at  -0.6861  -0.7605  -0.7429    1.125
1388142_at    -0.7516  -0.6713  -0.7660    1.124
> pheatmap(n,show_colnames =F,show_rownames = F)
> #显示分组
> ac=data.frame(group=group_list_1)
> rownames(ac)=colnames(n)
> pheatmap(n,show_colnames =T,show_rownames = F,cluster_col = F,annotation_col = ac,color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
#保存
>pdf(file = "heatmap.pdf")