GEO数据挖掘二

生信技能树学徒第二周

一、GEO数据库 芯片数据获取

1.  
   #数据下载
2.  
   rm(list = ls())
3.  
   library(GEOquery)
4.  
   #先去网页确定是否是表达芯片数据，不是的话不能用本流程。以GSE56649为例
5.  
   gse_number = "GSE56649"
6.  
   eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
7.  
   class(eSet)
8.  
   length(eSet)
9.  
   #得出往往是一个长度为一的列表，转换为数据框的格式
10.  
    eSet = eSet[[1]]
11.  
    #(1)提取表达矩阵exp
12.  
    exp <- exprs(eSet)
13.  
    dim(exp)
14.  
    exp[1:4,1:4]
15.  
    #检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。
16.  
    #如果表达矩阵为空，大多数是转录组数据，不能用这个流程。
17.  
    #自行判断是否需要log，一般boxplot 后的纵坐标范围在0-20之间，如果过高就需要取log
18.  
    #exp = log2(exp 1)
19.  
    boxplot(exp)
20.  
    #(2)提取临床信息，会有分组信息，我们根据分组信息来进行命名
21.  
    pd <- pData(eSet)
22.  
    #(3)让exp列名与pd的行名顺序完全一致，
23.  
    p = identical(rownames(pd),colnames(exp));p
24.  
    if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
25.  
    #(4)提取芯片平台编号
26.  
    gpl_number <- eSet@annotation;gpl_number
27.  
    save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

      因此，数据获取并确认的流程，getGEO函数下载数据→转换格式→提取表达矩阵→提取分组信息→确认分组对应

二、确定分组以及探针注释的获取

数据默认分组时冗长的，我们一般用control和treat来进行重新命名分组

1.  
   # 标准流程代码是二分组
2.  
   # 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
3.  
   if(F){
4.  
   # 1.Group----
5.  
   # 第一种方法，有现成的可以用来分组的列
6.  
   Group = pd$`disease state:ch1`
7.  
   }else if(F){
8.  
   # 第二种方法，自己生成
9.  
   Group = c(rep("RA",times=13),
10.  
    rep("control",times=9))
11.  
    Group = rep(c("RA","control"),times = c(13,9))
12.  
    }else if(T){
13.  
    # 第三种方法，使用字符串出理的函数获取分组
14.  
    Group=ifelse(str_detect(pd$source_name_ch1,"control"),
15.  
    "control",
16.  
    "RA")
17.  
    }
18.  
    #结束之后及时进行检查
19.  
    # 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
20.  
    Group = factor(Group,levels = c("control","RA"))#确定control的位置，总是要把它放在第一的位置上！关键词相应
21.  
    Group

探针对应到基因里

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

三、PCA和热图绘制

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
7.  
    
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")
22.  
     
23.  
    # 2.top 1000 sd 热图----
24.  
    cg=names(tail(sort(apply(exp,1,sd)),1000))
25.  
    n=exp[cg,]
26.  
     
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.  
    )
36.  
     
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)#规定取值对应颜色的而范围以及分1000个级别
44.  
    #cluster_rows = TRUE 是否聚类
45.  
    #cluster_cols = TRUE
46.  
    )
47.  
    dev.off()

四、差异分析

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)
10.  
     
11.  
    #为deg数据框添加几列
12.  
    #1.加probe_id列，把行名变成一列
13.  
    library(dplyr)
14.  
    deg <- mutate(deg,probe_id=rownames(deg))
15.  
    #2.加上探针注释，去除一个基因对应多种探针ID，随机去重
16.  
    ids = ids[!duplicated(ids$symbol),]
17.  
    #其他去重方式在zz.去重方式.R
18.  
    deg <- inner_join(deg,ids,by="probe_id")
19.  
    nrow(deg)
20.  
     
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)
28.  
    #4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
29.  
    library(clusterProfiler)
30.  
    library(org.Hs.eg.db)
31.  
    s2e <- bitr(deg$symbol,
32.  
    fromType = "SYMBOL",
33.  
    toType = "ENTREZID",
34.  
    OrgDb = org.Hs.eg.db)#人类
35.  
    #其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
36.  
    deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
37.  
    save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

      差异分析主要是用limma包来做的，通过lmFIT和ebays(用来拟合和校正t检验之中的方差部分)，之后添加探针注释，标记上下调基因，加RNTREZID列来用于差异分析

五、差异分析后画图部分

1.  
   rm(list = ls())
2.  
   load(file = "step1output.Rdata")
3.  
   load(file = "step4output.Rdata")
4.  
   #1.火山图----
5.  
   library(dplyr)
6.  
   library(ggplot2)
7.  
   #这次去重复是去除多个表达量对应一个基因的情况，随机去重
8.  
   dat = deg[!duplicated(deg$symbol),]
9.  
   ## 规定横纵轴
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
21.  
    #优秀的管道符号，将一个函数的运算结果传输到下一个函数中
22.  
    for_label <- dat%>%
23.  
    可以更换自己关注的基因名称
24.  
    filter(symbol %in% c("HADHA","LRRFIP1"))##
25.  
    # 火山图
26.  
    volcano_plot <- p
27.  
    geom_point(size = 3, shape = 1, data = for_label)
28.  
    ggrepel::geom_label_repel(
29.  
    aes(label = symbol),
30.  
    data = for_label,
31.  
    color="black"
32.  
    )
33.  
    volcano_plot
34.  
     
35.  
    #2.差异基因热图----
36.  
     
37.  
    load(file = 'step2output.Rdata')
38.  
    # 表达矩阵行名替换
39.  
    exp = exp[dat$probe_id,]
40.  
    rownames(exp) = dat$symbol
41.  
    if(F){
42.  
    #全部差异基因
43.  
    cg = dat$symbol[dat$change !="stable"]
44.  
    length(cg)
45.  
    }else{
46.  
    #取前10上调和前10下调
47.  
    library(dplyr)
48.  
    dat2 = dat %>%
49.  
    filter(change!="stable") %>%
50.  
    arrange(logFC)
51.  
    cg = c(head(dat2$symbol,10),
52.  
    tail(dat2$symbol,10))
53.  
    }
54.  
    n=exp[cg,]
55.  
    dim(n)
56.  
     
57.  
    #差异基因热图
58.  
    library(pheatmap)
59.  
    annotation_col=data.frame(group=Group)
60.  
    rownames(annotation_col)=colnames(n)
61.  
    heatmap_plot <- pheatmap(n,show_colnames =F,
62.  
    scale = "row",
63.  
    #cluster_cols = F,
64.  
    annotation_col=annotation_col,
65.  
    breaks = seq(-3,3,length.out = 100)
66.  
    )
67.  
    heatmap_plot
68.  
    #拼图
69.  
    library(patchwork)
70.  
    library(ggplotify)
71.  
    volcano_plot as.ggplot(heatmap_plot)
72.  
     
73.  
    # 3.感兴趣基因的相关性----
74.  
    library(corrplot)
75.  
    g = deg$symbol[1:10] # 换成自己感兴趣的基因
76.  
    g
77.  
     
78.  
    M = cor(t(exp[g,]))
79.  
    pheatmap(M)
80.  
     
81.  
    library(paletteer)
82.  
    my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
83.  
    my_color = colorRampPalette(my_color)(10)
84.  
    corrplot(M, type="upper",
85.  
    method="pie",
86.  
    order="hclust",
87.  
    col=my_color,
88.  
    tl.col="black",
89.  
    tl.srt=45)
90.  
    library(cowplot)
91.  
    cor_plot <- recordPlot()
92.  
     
93.  
    # 拼图
94.  
    load("pca_plot.Rdata")
95.  
    plot_grid(pca_plot,cor_plot,
96.  
    volcano_plot,heatmap_plot$gtable)
97.  
    dev.off()
98.  
    #保存
99.  
    pdf("deg.pdf")
100.  
     plot_grid(pca_plot,cor_plot,
101.  
     volcano_plot,heatmap_plot$gtable)
102.  
     dev.off()

这篇好文章是转载于：学新通技术网