image.png

  1. load("TCGA-CHOL_DEG.Rdata")
  2. load("TCGA-CHOL.Rdata")
  3. cg2 = rownames(DEG2)[DEG2$change !="NOT"]
  4. library(tinyarray)
  5. library(edgeR)
  6. dat = log2(cpm(exp)+1)
  7. draw_heatmap(dat[cg2,],Group,cluster_cols = F)
  8. dat2
  9. Group2 = sort(Group)
  10. dat2 = dat[,order(Group)]
  11. draw_heatmap(dat2[cg2,],Group2,cluster_cols = F)

三大R包差异分析

  1. rm(list = ls())
  2. load("DHA.Rdata")
  3. table(Group)
  4. #> Group
  5. #> DMSO  DHA 
  6. #>    3    3
  7. #deseq2----
  8. library(DESeq2)
  9. colData <- data.frame(row.names =colnames(exp), 
  10.                       condition=Group)
  11. if(!file.exists(paste0(proj,"_dd.Rdata"))){
  12.   dds <- DESeqDataSetFromMatrix(
  13.   countData = exp,
  14.   colData = colData,
  15.   design = ~ condition)
  16.   dds <- DESeq(dds)
  17.   save(dds,file = paste0(proj,"_dd.Rdata"))
  18. }
  19. load(file = paste0(proj,"_dd.Rdata"))
  20. class(dds)
  21. #> [1] "DESeqDataSet"
  22. #> attr(,"package")
  23. #> [1] "DESeq2"
  24. res <- results(dds, contrast = c("condition",rev(levels(Group))))
  25. #constrast
  26. c("condition",rev(levels(Group)))
  27. #> [1] "condition" "DHA"       "DMSO"
  28. class(res)
  29. #> [1] "DESeqResults"
  30. #> attr(,"package")
  31. #> [1] "DESeq2"
  32. DEG1 <- as.data.frame(res)
  33. DEG1 <- DEG1[order(DEG1$pvalue),] 
  34. DEG1 = na.omit(DEG1)
  35. head(DEG1)
  36. #>        baseMean log2FoldChange      lfcSE      stat        pvalue          padj
  37. #> CDH6   5006.549      -4.030551 0.10671818 -37.76818  0.000000e+00  0.000000e+00
  38. #> CXCL8 22087.170       4.958524 0.12501648  39.66296  0.000000e+00  0.000000e+00
  39. #> HMOX1 46652.918       2.995161 0.08194416  36.55124 1.704000e-292 1.163094e-288
  40. #> THBD   2832.828       5.025714 0.13858308  36.26499 5.768388e-288 2.952982e-284
  41. #> LIF    4105.206       4.812167 0.14857784  32.38819 4.025761e-230 1.648710e-226
  42. #> DUSP5  2620.024       4.334570 0.13469580  32.18044 3.314657e-227 1.131237e-223
  44. #添加change列标记基因上调下调
  45. logFC_t = 2
  46. pvalue_t = 0.05
  48. k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
  49. #> k1
  50. #> FALSE  TRUE 
  51. #> 19726   751
  52. k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
  53. #> k2
  54. #> FALSE  TRUE 
  55. #> 19858   619
  56. DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  57. table(DEG1$change)
  58. #> 
  59. #>  DOWN   NOT    UP 
  60. #>   751 19107   619
  61. head(DEG1)
  62. #>        baseMean log2FoldChange      lfcSE      stat        pvalue          padj
  63. #> CDH6   5006.549      -4.030551 0.10671818 -37.76818  0.000000e+00  0.000000e+00
  64. #> CXCL8 22087.170       4.958524 0.12501648  39.66296  0.000000e+00  0.000000e+00
  65. #> HMOX1 46652.918       2.995161 0.08194416  36.55124 1.704000e-292 1.163094e-288
  66. #> THBD   2832.828       5.025714 0.13858308  36.26499 5.768388e-288 2.952982e-284
  67. #> LIF    4105.206       4.812167 0.14857784  32.38819 4.025761e-230 1.648710e-226
  68. #> DUSP5  2620.024       4.334570 0.13469580  32.18044 3.314657e-227 1.131237e-223
  69. #>       change
  70. #> CDH6    DOWN
  71. #> CXCL8     UP
  72. #> HMOX1     UP
  73. #> THBD      UP
  74. #> LIF       UP
  75. #> DUSP5     UP
  77. #edgeR----
  78. library(edgeR)
  80. dge <- DGEList(counts=exp,group=Group)
  81. dge$samples$lib.size <- colSums(dge$counts)
  82. dge <- calcNormFactors(dge) 
  84. design <- model.matrix(~Group)
  86. dge <- estimateGLMCommonDisp(dge, design)
  87. dge <- estimateGLMTrendedDisp(dge, design)
  88. dge <- estimateGLMTagwiseDisp(dge, design)
  90. fit <- glmFit(dge, design)
  91. fit <- glmLRT(fit) 
  93. DEG2=topTags(fit, n=Inf)
  94. class(DEG2)
  95. #> [1] "TopTags"
  96. #> attr(,"package")
  97. #> [1] "edgeR"
  98. DEG2=as.data.frame(DEG2)
  99. head(DEG2)
  100. #>             logFC   logCPM       LR        PValue           FDR
  101. #> LIF      4.780301 7.069581 571.9517 2.112394e-126 4.695851e-122
  102. #> PTGS2    5.970830 7.152826 565.5336 5.258853e-125 5.845215e-121
  103. #> THBD     4.992993 6.534219 564.2212 1.014763e-124 7.519395e-121
  104. #> CDH6    -4.062412 7.384277 532.6938 7.326559e-118 4.071735e-114
  105. #> CXCL8    4.926708 9.496259 470.0676 3.097901e-104 1.377327e-100
  106. #> TMEM158  5.147555 6.279693 437.4408  3.903283e-97  1.446166e-93
  108. k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t)
  109. k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t)
  110. DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  112. head(DEG2)
  113. #>             logFC   logCPM       LR        PValue           FDR change
  114. #> LIF      4.780301 7.069581 571.9517 2.112394e-126 4.695851e-122     UP
  115. #> PTGS2    5.970830 7.152826 565.5336 5.258853e-125 5.845215e-121     UP
  116. #> THBD     4.992993 6.534219 564.2212 1.014763e-124 7.519395e-121     UP
  117. #> CDH6    -4.062412 7.384277 532.6938 7.326559e-118 4.071735e-114   DOWN
  118. #> CXCL8    4.926708 9.496259 470.0676 3.097901e-104 1.377327e-100     UP
  119. #> TMEM158  5.147555 6.279693 437.4408  3.903283e-97  1.446166e-93     UP
  120. table(DEG2$change)
  121. #> 
  122. #>  DOWN   NOT    UP 
  123. #>   936 20586   708
  124. ###limma----
  125. library(limma)
  126. dge <- edgeR::DGEList(counts=exp)
  127. dge <- edgeR::calcNormFactors(dge)
  128. design <- model.matrix(~Group)
  129. v <- voom(dge,design, normalize="quantile")
  131. fit <- lmFit(v, design)
  132. fit= eBayes(fit)
  134. DEG3 = topTable(fit, coef=2, n=Inf)
  135. DEG3 = na.omit(DEG3)
  137. k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t)
  138. k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t)
  139. DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
  140. table(DEG3$change)
  141. #> 
  142. #>  DOWN   NOT    UP 
  143. #>   846 20611   773
  144. head(DEG3)
  145. #>            logFC   AveExpr         t      P.Value    adj.P.Val        B change
  146. #> HMOX1   3.023387  9.933877  59.51263 1.613427e-10 1.895736e-06 15.06135     UP
  147. #> CXCL8   5.010797  8.042577  58.11745 1.896200e-10 1.895736e-06 14.56756     UP
  148. #> PLIN2   3.029546  9.754344  54.12555 3.077709e-10 1.895736e-06 14.47996     UP
  149. #> SQSTM1  2.528669 10.930792  51.44786 4.346774e-10 1.895736e-06 14.22328     UP
  150. #> AKR1C1  3.492536  8.480250  50.22935 5.116695e-10 1.895736e-06 13.91592     UP
  151. #> CDH6   -4.028740  6.272248 -52.96068 3.568837e-10 1.895736e-06 13.67540   DOWN
  153. tj = data.frame(deseq2 = as.integer(table(DEG1$change)),
  154.            edgeR = as.integer(table(DEG2$change)),
  155.            limma_voom = as.integer(table(DEG3$change)),
  156.            row.names = c("down","not","up")
  157.           );tj
  158. #>      deseq2 edgeR limma_voom
  159. #> down    751   936        846
  160. #> not   19107 20586      20611
  161. #> up      619   708        773
  162. save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))

可视化

  1. library(ggplot2)
  2. library(tinyarray)
  3. exp[1:4,1:4]
  4. #>         DMSO1 DMSO2 DMSO3 S25uMDNHA1
  5. #> MADD     1097  1586  1055       1527
  6. #> MAST2    2765  2951  2670       5872
  7. #> TRBJ2-5     2     1     1          1
  8. #> TRBJ2-4     1     1     1          1
  9. # cpm 去除文库大小的影响
  10. dat = log2(cpm(exp)+1)
  11. pca.plot = draw_pca(dat,Group);pca.plot
  1. save(pca.plot,file = paste0(proj,"_pcaplot.Rdata"))
  3. cg1 = rownames(DEG1)[DEG1$change !="NOT"]
  4. cg2 = rownames(DEG2)[DEG2$change !="NOT"]
  5. cg3 = rownames(DEG3)[DEG3$change !="NOT"]
  7. h1 = draw_heatmap(dat[cg1,],Group)
  8. h2 = draw_heatmap(dat[cg2,],Group)
  9. h3 = draw_heatmap(dat[cg3,],Group)
  11. v1 = draw_volcano(DEG1,pkg = 1,logFC_cutoff = logFC_t)
  12. v2 = draw_volcano(DEG2,pkg = 2,logFC_cutoff = logFC_t)
  13. v3 = draw_volcano(DEG3,pkg = 3,logFC_cutoff = logFC_t)
  15. library(patchwork)
  16. (h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect') &theme(legend.position = "none")
  1. ggsave(paste0(proj,"_heat_vo.png"),width = 15,height = 10)

三大R包差异基因对比

  1. UP=function(df){
  2.   rownames(df)[df$change=="UP"]
  3. }
  4. DOWN=function(df){
  5.   rownames(df)[df$change=="DOWN"]
  6. }
  8. up = intersect(intersect(UP(DEG1),UP(DEG2)),UP(DEG3))
  9. down = intersect(intersect(DOWN(DEG1),DOWN(DEG2)),DOWN(DEG3))
  10. dat = log2(cpm(exp)+1)
  11. hp = draw_heatmap(dat[c(up,down),],Group)
  13. #上调、下调基因分别画维恩图
  14. up_genes = list(Deseq2 = UP(DEG1),
  15.           edgeR = UP(DEG2),
  16.           limma = UP(DEG3))
  18. down_genes = list(Deseq2 = DOWN(DEG1),
  19.           edgeR = DOWN(DEG2),
  20.           limma = DOWN(DEG3))
  22. up.plot <- draw_venn(up_genes,"UPgene")
  23. down.plot <- draw_venn(down_genes,"DOWNgene")
  25. #维恩图拼图，终于搞定
  27. library(patchwork)
  28. #up.plot + down.plot
  29. # 拼图
  30. pca.plot + hp+up.plot +down.plot+ plot_layout(guides = "collect")
  31. ggsave(paste0(proj,"_heat_ve_pca.png"),width = 15,height = 10)

分组聚类的热图

  1. library(ComplexHeatmap)
  2. library(circlize)
  3. col_fun = colorRamp2(c(-2, 0, 2), c("#2fa1dd", "white", "#f87669"))
  4. top_annotation = HeatmapAnnotation(
  5.   cluster = anno_block(gp = gpar(fill = c("#f87669","#2fa1dd")),
  6.                        labels = levels(Group),
  7.                        labels_gp = gpar(col = "white", fontsize = 12)))
  9. m = Heatmap(t(scale(t(exp[c(up,down),]))),name = " ",
  10.             col = col_fun,
  11.         top_annotation = top_annotation,
  12.         column_split = Group,
  13.         show_heatmap_legend = T,
  14.         border = F,
  15.         show_column_names = F,
  16.         show_row_names = F,
  17.         use_raster = F,
  18.         cluster_column_slices = F,
  19.         column_title = NULL)
  20. m