01、02是重点，需着重处理好，后续才能顺利，可形成模版
image.png

安装

  1. options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
  2. if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
  3. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
  5. cran_packages <- c('tidyr',
  6.                    'tibble',
  7.                    'dplyr',
  8.                    'stringr',
  9.                    'ggplot2',
  10.                    'ggpubr',
  11.                    'factoextra',
  12.                    'FactoMineR',
  13.                    'devtools',
  14.                    'cowplot',
  15.                    'patchwork',
  16.                    'basetheme',
  17.                    'paletteer',
  18.                    'AnnoProbe',
  19.                    'ggthemes',
  20.                    'VennDiagram',
  21.                    'tinyarray') 
  22. Biocductor_packages <- c('GEOquery',
  23.                          'hgu133plus2.db',
  24.                          'ggnewscale',
  25.                          "limma",
  26.                          "impute",
  27.                          "GSEABase",
  28.                          "GSVA",
  29.                          "clusterProfiler",
  30.                          "org.Hs.eg.db",
  31.                          "preprocessCore",
  32.                          "enrichplot",
  33.                          "ggplotify")
  35. for (pkg in cran_packages){
  36.   if (! require(pkg,character.only=T) ) {
  37.     install.packages(pkg,ask = F,update = F)
  38.     require(pkg,character.only=T) 
  39.   }
  40. }
  43. for (pkg in Biocductor_packages){
  44.   if (! require(pkg,character.only=T) ) {
  45.     BiocManager::install(pkg,ask = F,update = F)
  46.     require(pkg,character.only=T) 
  47.   }
  48. }
  50. #前面的所有提示和报错都先不要管。主要看这里
  51. for (pkg in c(Biocductor_packages,cran_packages)){
  52.   require(pkg,character.only=T) 
  53. }
  54. #没有任何提示就是成功了，如果有warningxx包不存在，用library检查一下。
  56. #library报错，就单独安装。
  57. library(hgu133plus2.db)

下载及检查数据

  1. #切记先听课再跑代码
  2. #数据下载
  3. rm(list = ls())
  4. library(GEOquery)
  5. #先去网页确定是否是表达芯片数据，不是的话不能用本流程。
  6. #功能参数不需要更改
  7. #getGEO具有判断机制，如果本地有，则不会再下载
  8. gse_number = "GSE56649"
  9. eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
  10. class(eSet)
  11. length(eSet)
  12. eSet = eSet[[1]]
  13. #(1)提取表达矩阵exp
  14. exp <- exprs(eSet)
  15. dim(exp)
  16. exp[1:4,1:4]#列出前4行4列的情况
  17. #检查矩阵是否正常，如果是空的就会报错。对于空的、负值的、异常值的矩阵需要处理原始数据
  18. #如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。
  19. #自行判断是否需要log。一般数值<20，则取过log；若是几百上千，则未取
  20. exp = log2(exp+1)
  21. boxplot(exp)
  22. #(2)提取临床信息
  23. pd <- pData(eSet)
  24. #(3)让exp列名与pd的行名顺序完全一致
  25. p = identical(rownames(pd),colnames(exp));p
  26. if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
  27. #(4)提取芯片平台编号
  28. gpl_number <- eSet@annotation;gpl_number
  29. save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

image.png
image.png
image.png

实验分组&探针注释

  1. # Group(实验分组)和ids(探针注释)
  2. rm(list = ls())  
  3. load(file = "step1output.Rdata")
  4. library(stringr)
  5. # 标准流程代码是二分组，多分组数据的分析后面另讲
  6. # 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
  7. if(F){
  8.   # 1.Group----
  9.   # 第一种方法，有现成的可以用来分组的列
  10.   Group = pd$`disease state:ch1` 
  11. }else if(F){
  12.   # 第二种方法，自己生成
  13.   Group = c(rep("RA",times=13),
  14.             rep("control",times=9))
  15.   Group = rep(c("RA","control"),times = c(13,9))
  16. }else if(T){
  17.   # 第三种方法，使用字符串出理的函数获取分组
  18.   Group=ifelse(str_detect(pd$source_name_ch1,"control"),
  19.                "control",
  20.                "RA")
  21. }
  23. # 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
  24. Group = factor(Group,levels = c("control","RA"))
  25. Group

探针注释：探针与基因的对应关系

image.png
gpl文件解析
image.png
自主注释流程
image.png
可能存在一个探针对应多个基因
image.png

  1. #2.探针注释的获取-----------------
  2. #捷径
  3. library(tinyarray)
  4. find_anno(gpl_number)
  5. ids <- AnnoProbe::idmap('GPL570')
  7. #四种方法，方法1里找不到就从方法2找，以此类推。
  8. #不同的方法，获取的数值可能有些差异，但被允许
  9. #方法1 BioconductorR包(最常用)
  10. gpl_number 
  11. #http://www.bio-info-trainee.com/1399.html
  12. if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
  13. library(hgu133plus2.db)
  14. ls("package:hgu133plus2.db")
  15. ids <- toTable(hgu133plus2SYMBOL)
  16. head(ids)
  18. # 方法2 读取GPL网页的表格文件，按列取子集
  19. ##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
  20. if(F){
  21.   #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  22.   b = read.table("GPL570-55999.txt",header = T,
  23.                  quote = "\"",sep = "\t",check.names = F)
  24.   colnames(b)
  25.   ids2 = b[,c("ID","Gene Symbol")] #按列名取子集
  26.   colnames(ids2) = c("probe_id","symbol")
  27.   ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]
  28. } #去除存在///的探针
  30. # 方法3 官网下载注释文件并读取
  31. ##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
  33. # 方法4 自主注释 
  34. #https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
  35. save(exp,Group,ids,gse_number,file = "step2output.Rdata")

多个探针对应1个基因
image.png
image.png

差异分析图

  1. rm(list = ls()) 
  2. load(file = "step1output.Rdata")
  3. load(file = "step4output.Rdata")
  5. #1.火山图----
  6. library(dplyr)
  7. library(ggplot2)
  8. dat  = deg[!duplicated(deg$symbol),]
  10. p <- ggplot(data = dat, 
  11.             aes(x = logFC, 
  12.                 y = -log10(P.Value))) +
  13.   geom_point(alpha=0.4, size=3.5, 
  14.              aes(color=change)) +
  15.   ylab("-log10(Pvalue)")+
  16.   scale_color_manual(values=c("blue", "grey","red"))+
  17.   geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  18.   geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  19.   theme_bw()
  20. p
  22. for_label <- dat%>% 
  23.   filter(symbol %in% c("HADHA","LRRFIP1")) #添加指定基因
  25. volcano_plot <- p +
  26.   geom_point(size = 3, shape = 1, data = for_label) +
  27.   ggrepel::geom_label_repel(
  28.     aes(label = symbol),
  29.     data = for_label,
  30.     color="black"
  31.   ) #对数据进行筛选，允许对不同图层进行特定设置
  32. volcano_plot
  34. #2.差异基因热图----
  36. load(file = 'step2output.Rdata')
  37. # 表达矩阵行名替换
  38. exp = exp[dat$probe_id,]
  39. rownames(exp) = dat$symbol
  40. if(F){
  41.   #全部差异基因
  42.   cg = dat$symbol[dat$change !="stable"]
  43.   length(cg)
  44. }else{
  45.   #取前10上调和前10下调
  46.   library(dplyr)
  47.   dat2 = dat %>%
  48.     filter(change!="stable") %>%
  49.     arrange(logFC)
  50.   cg = c(head(dat2$symbol,10),
  51.          tail(dat2$symbol,10))
  52.   }
  53. n=exp[cg,]
  54. dim(n)
  55. #差异基因热图
  56. library(pheatmap)
  57. annotation_col=data.frame(group=Group)
  58. rownames(annotation_col)=colnames(n) 
  59. heatmap_plot <- pheatmap(n,show_colnames =F,
  60.                          scale = "row",
  61.                          #cluster_cols = F, 
  62.                          annotation_col=annotation_col,
  63.                          breaks = seq(-3,3,length.out = 100)
  64. ) 
  65. heatmap_plot
  66. #拼图
  67. library(patchwork)
  68. library(ggplotify)
  69. volcano_plot +as.ggplot(heatmap_plot)
  71. # 3.感兴趣基因的相关性----
  72. library(corrplot)
  73. g = deg$symbol[1:10] # 换成自己感兴趣的基因
  74. g
  76. M = cor(t(exp[g,]))
  77. pheatmap(M)
  79. library(paletteer)
  80. my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
  81. my_color = colorRampPalette(my_color)(10)
  82. corrplot(M, type="upper",
  83.          method="pie",
  84.          order="hclust", 
  85.          col=my_color,
  86.          tl.col="black", 
  87.          tl.srt=45)
  88. library(cowplot)
  89. cor_plot <- recordPlot()#抠图
  90. # 拼图
  91. load("pca_plot.Rdata")
  92. plot_grid(pca_plot,cor_plot,
  93.           volcano_plot,heatmap_plot$gtable)
  94. dev.off()
  95. #保存
  96. pdf("deg.pdf", width = 12,height = 12 )#修改图的大小
  97. plot_grid(pca_plot,cor_plot,
  98.           volcano_plot,heatmap_plot$gtable)
  99. dev.off()

主成分分析

image.png
image.png
image.png

  1. rm(list = ls())  
  2. load(file = "step1output.Rdata")
  3. load(file = "step2output.Rdata")
  4. #输入数据：exp和Group
  5. #Principal Component Analysis
  6. #http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
  8. # 1.PCA 图----
  9. dat=as.data.frame(t(exp))
  10. library(FactoMineR)
  11. library(factoextra) 
  12. dat.pca <- PCA(dat, graph = FALSE)
  13. pca_plot <- fviz_pca_ind(dat.pca,
  14.                          geom.ind = "point", # show points only (nbut not "text")
  15.                          col.ind = Group, # color by groups
  16.                          palette = c("#00AFBB", "#E7B800"),
  17.                          addEllipses = TRUE, # Concentration ellipses
  18.                          legend.title = "Groups"
  19. )
  20. pca_plot
  21. save(pca_plot,file = "pca_plot.Rdata")
  23. # 2.top 1000 sd 热图---- 
  24. cg=names(tail(sort(apply(exp,1,sd)),1000))
  25. n=exp[cg,]
  27. # 直接画热图，对比不鲜明
  28. library(pheatmap)
  29. annotation_col=data.frame(group=Group) #指示条
  30. rownames(annotation_col)=colnames(n) 
  31. pheatmap(n,
  32.          show_colnames =F,
  33.          show_rownames = F,
  34.          annotation_col=annotation_col
  35. )
  37. # 按行标准化
  38. pheatmap(n,
  39.          show_colnames =F,
  40.          show_rownames = F,
  41.          annotation_col=annotation_col,
  42.          scale = "row", #实现标准化
  43.          breaks = seq(-3,3,length.out = 100) #设置色带
  44.          ) 
  45. dev.off()
  47. # 关于scale的进一步学习：zz.scale.R

DEG差异分析、可视化

  1. rm(list = ls()) 
  2. load(file = "step2output.Rdata")
  3. #差异分析，用limma包来做
  4. #需要表达矩阵和Group，不需要改
  5. library(limma)
  6. design=model.matrix(~Group)
  7. fit=lmFit(exp,design)
  8. fit=eBayes(fit)
  9. deg=topTable(fit,coef=2,number = Inf)
  11. #为deg数据框添加几列
  12. #1.加probe_id列，把行名变成一列
  13. library(dplyr)
  14. deg <- mutate(deg,probe_id=rownames(deg))
  15. #2.加上探针注释
  16. ids = ids[!duplicated(ids$symbol),] #随机去重
  17. #其他去重方式在zz.去重方式.R
  18. deg <- inner_join(deg,ids,by="probe_id")
  19. nrow(deg)
  21. #3.加change列,标记上下调基因
  22. logFC_t=1
  23. P.Value_t = 0.05
  24. k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
  25. k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
  26. deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
  27. table(deg$change)
  29. #4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
  30. library(clusterProfiler)
  31. library(org.Hs.eg.db)
  32. s2e <- bitr(deg$symbol, 
  33.             fromType = "SYMBOL",
  34.             toType = "ENTREZID",
  35.             OrgDb = org.Hs.eg.db)#人类
  36. #其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
  37. deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
  38. save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

富集分析(GSEA)

  1. rm(list = ls())  
  2. load(file = 'step4output.Rdata')
  3. library(clusterProfiler)
  4. library(ggthemes)
  5. library(org.Hs.eg.db)
  6. library(dplyr)
  7. library(ggplot2)
  8. library(stringr)
  9. library(enrichplot)
  11. # 1.GO 富集分析----
  13. #(1)输入数据
  14. gene_up = deg$ENTREZID[deg$change == 'up'] 
  15. gene_down = deg$ENTREZID[deg$change == 'down'] 
  16. gene_diff = c(gene_up,gene_down)
  18. #(2)富集
  19. #以下步骤耗时很长，设置了存在即跳过
  20. if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  21.   ego <- enrichGO(gene = gene_diff,
  22.                   OrgDb= org.Hs.eg.db,
  23.                   ont = "ALL",
  24.                   readable = TRUE)
  25.   ego_BP <- enrichGO(gene = gene_diff,
  26.                   OrgDb= org.Hs.eg.db,
  27.                   ont = "BP",
  28.                   readable = TRUE)
  29.   #ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  30.   save(ego,ego_BP,file = paste0(gse_number,"_GO.Rdata"))
  31. }
  32. load(paste0(gse_number,"_GO.Rdata"))
  34. #(3)可视化
  35. #条带图
  36. barplot(ego)
  37. barplot(ego, split = "ONTOLOGY", font.size = 10, 
  38.         showCategory = 5) + 
  39.   facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
  41. #气泡图
  42. dotplot(ego)
  44. dotplot(ego, split = "ONTOLOGY", font.size = 10, 
  45.         showCategory = 5) + 
  46.   facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 
  48. #(3)展示top通路的共同基因，要放大看。
  50. #gl 用于设置下图的颜色
  51. gl = deg$logFC
  52. names(gl)=deg$ENTREZID
  53. #Gene-Concept Network
  54. cnetplot(ego,categorySize="pvalue", foldChange=gl,colorEdge = TRUE,showCategory = 3)
  55. cnetplot(ego, showCategory = 3,foldChange=gl, circular = TRUE, colorEdge = TRUE)
  57. # 2.KEGG pathway analysis----
  58. #上调、下调、差异、所有基因
  59. #（1）输入数据
  60. gene_up = deg[deg$change == 'up','ENTREZID'] 
  61. gene_down = deg[deg$change == 'down','ENTREZID'] 
  62. gene_diff = c(gene_up,gene_down)
  64. #（2）对上调/下调/所有差异基因进行富集分析
  66. if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
  67.   kk.up <- enrichKEGG(gene         = gene_up,
  68.                       organism     = 'hsa')
  69.   kk.down <- enrichKEGG(gene         =  gene_down,
  70.                         organism     = 'hsa')
  71.   kk.diff <- enrichKEGG(gene         = gene_diff,
  72.                         organism     = 'hsa')
  73.   save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
  74. }
  75. load(paste0(gse_number,"_KEGG.Rdata"))
  77. #(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
  78. table(kk.diff@result$p.adjust<0.05)
  79. table(kk.up@result$p.adjust<0.05)
  80. table(kk.down@result$p.adjust<0.05)
  82. #(4)双向图
  83. # 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值，在文中说明即可
  84. down_kegg <- kk.down@result %>%
  85.   filter(pvalue<0.05) %>% #筛选行
  86.   mutate(group=-1) #新增列
  88. up_kegg <- kk.up@result %>%
  89.   filter(pvalue<0.05) %>%
  90.   mutate(group=1)
  92. source("kegg_plot_function.R")
  93. g_kegg <- kegg_plot(up_kegg,down_kegg)
  94. g_kegg
  95. #g_kegg +scale_y_continuous(labels = c(2,0,2,4,6))
  97. # 3.GSEA作kegg和GO富集分析----
  98. #https://www.yuque.com/xiaojiewanglezenmofenshen/dbwkg1/ytawgg
  100. #(1)查看示例数据
  101. data(geneList, package="DOSE")
  102. #(2)将我们的数据转换成示例数据的格式
  103. geneList=deg$logFC
  104. names(geneList)=deg$ENTREZID
  105. geneList=sort(geneList,decreasing = T)
  106. #(3)GSEA富集
  107. kk_gse <- gseKEGG(geneList     = geneList,
  108.                   organism     = 'hsa',
  109.                   verbose      = FALSE)
  110. down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
  111. up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
  112. #(4)可视化
  113. g2 = kegg_plot(up_kegg,down_kegg)
  114. g2
  115. # 4.能看懂的资料越来越多----
  116. # GSEA学习更多：https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
  117. # 富集分析学习更多：http://yulab-smu.top/clusterProfiler-book/index.html
  118. # 弦图：https://www.yuque.com/xiaojiewanglezenmofenshen/dbwkg1/dgs65p
  119. # GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
  120. # 网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~
  1. if(F){
  2.   a = 1 #限速步骤
  3.   save(a,file = "a.Rdata")
  4. }
  5. load("a.Rdata")
  7. if(!file.exists("a.Rdata")){
  8.   a = 1 #限速步骤
  9.   save(a,file = "a.Rdata")
  10. }
  11. load("a.Rdata")

问题和解答

image.png

image.png