---limma分析差异基因---
在经过了前两期中的数据下载，数据基本处理之后，解决了一个探针对应多个基因数的以及多个探针对应一个基因求平均值，在此基础上运用limma包分析差异基因；
除此以外，包括绘制火山图，热图，PCA等，都在本文中解决
数据载入
1if(T){
2Sys.setlocale(\"LC_ALL\",\"C\")
3library(dplyr)
4##
5if(T){
6load("expma.Rdata")
7load("probe.Rdata")
8}
9expma[1:5,1:5]
10boxplot(expma)##看下表达情况
11metdata[1:5,1:5]
12head(probe)
13
14##查看GeneSymbol是否有重复
15table(duplicated(probe$`GeneSymbol`))##12549FALSE
16
17##整合注释信息到表达矩阵
18ID<-rownames(expma)
19expma<-apply(expma,1,function(x){log2(x+1)})
20expma<-t(expma)
21eset<-expma[ID%in%probe$ID,]%>%cbind(probe)
22eset[1:5,1:5]
23colnames(eset)
24}
image.png1##[1]"GSM188013""GSM188014""GSM188016""GSM188018"
2##[5]"GSM188020""GSM188022""ID""GeneSymbol"
3##[9]"ENTREZ_GENE_ID"
多个探针求平均值
1test1<-aggregate(x=eset[,1:6],by=list(eset$`GeneSymbol`),FUN=mean,na.rm=T)
2test1[1:5,1:5]##与去重结果相吻合
3##Group.1GSM188013GSM188014GSM188016GSM188018
4##18.4388468.3685137.3224427.813573
5##2A1CF10.97902510.6169269.94077310.413311
6##3A2M6.5652766.4221128.1421945.652593
7##4A4GALT7.7286287.8189668.6798857.048563
8##5A4GNT10.24338810.1823829.3919918.779887
9dim(test1)##
10##[1]125497
11colnames(test1)
12##[1]"Group.1""GSM188013""GSM188014""GSM188016""GSM188018""GSM188020"
13##[7]"GSM188022"
14rownames(test1)<-test1$Group.1
15test1<-test1[,-1]
16eset_dat<-test1
PCA plot
Principal Component Analysis (PCA)分析使用的是基于R语言的 prcomp() and princomp()函数.
完成PCA分析一般有两种方法：princomp()使用的是一种的spectral decomposition方法
prcomp() and PCA()[FactoMineR]使用的是SVD法
1data<-eset_dat
2data[1:5,1:5]##表达矩阵
3##GSM188013GSM188014GSM188016GSM188018GSM188020
4##8.4388468.3685137.3224427.8135737.244615
5##A1CF10.97902510.6169269.94077310.4133119.743305
6##A2M6.5652766.4221128.1421945.6525935.550033
7##A4GALT7.7286287.8189668.6798857.0485635.929258
8##A4GNT10.24338810.1823829.3919918.7798879.431585
9metdata[1:5,1:5]
10##title
11##GSM188013DMSOtreatedMCF7breastcancercells[HG-U133A]Exp1
12##GSM188014DioxintreatedMCF7breastcancercells[HG-U133A]Exp1
13##GSM188016DMSOtreatedMCF7breastcancercells[HG-U133A]Exp2
14##GSM188018DioxintreatedMCF7breastcancercells[HG-U133A]Exp2
15##GSM188020DMSOtreatedMCF7breastcancercells[HG-U133A]Exp3
16##geo_accessionstatussubmission_date
17##GSM188013GSM188013PubliconMay12May08
18##GSM188014GSM188014PubliconMay12May08
19##GSM188016GSM188016PubliconMay12May08
20##GSM188018GSM188018PubliconMay12May08
21##GSM188020GSM188020PubliconMay12May08
22##last_update_date
23##GSM188013May11
24##GSM188014May11
25##GSM188016May11
26##GSM188018May11
27##GSM188020May11
28##构建group_list
29group_list<-rep(c("Treat","Control"),3)
30colnames(data)<-group_list
31library(factoextra)
32##Warning:package\"factoextra\"wasbuiltunderRversion3.5.3
33##Loadingrequiredpackage:ggplot2
34##Welcome!RelatedBooks:`PracticalGuideToClusterAnalysisinR`athttps://goo.gl/13EFCZ
35##计算PCA
36data<-t(data)##转换数据至行为sample,列为gene
37data<-as.data.frame(data)##注意数据要转换为数据框
38res.pca<-prcomp(data,scale=TRUE)
39##展示主成分对差异的贡献
40fviz_eig(res.pca)
Fig21##可视化结果
2fviz_pca_ind(res.pca,
3col.ind=group_list,#颜色对应group信息
4palette=c("#00AFBB","#FC4E07"),
5addEllipses=TRUE,#Concentrationellipses
6ellipse.type="confidence",
7legend.,##Legend名称
8repel=TRUE
9)
Fig3层次聚类图-聚类结果也与PCA相似，结果并不好
聚类分析的结果也同样可以进一步美化，但这里不做
计算距离时同样需进行转置，但在前一步PCA分析中的data已经经过转置，故未重复
1dd<-dist(data,method="euclidean")##data是经过行列转换的
2hc<-hclust(dd,method="ward.D2")
3plot(hc)
Fig41##对结果进行美化
2#Converthclustintoadendrogramandplot
3hcd<-as.dendrogram(hc)
4#DefinenodePar
5nodePar<-list(lab.cex=0.6,pch=c(NA,19),
6cex=0.7,col="blue")
7#Customizedplot;removelabels
8plot(hcd,ylab="Height",nodePar=nodePar,leaflab="none")
Fig5limma分析差异基因
limma包的具体用法参考 limma Users Guide
构建分组信息，构建好比较矩阵是关键
注意这里的表达矩阵信息 eset_dat是经过处理后的，为转置，行为gene，列为sample
1library(limma)
2library(dplyr)
3group_list
4##[1]"Treat""Control""Treat""Control""Treat""Control"
5design<-model.matrix(~0+factor(group_list))
6colnames(design)=levels(factor(group_list))
7design
8##ControlTreat
9##101
10##210
11##301
12##410
13##501
14##610
15##attr(,"assign")
16##[1]11
17##attr(,"contrasts")
18##attr(,"contrasts")$`factor(group_list)`
19##[1]"contr.treatment"
20##比较信息
21contrast.matrix<-makeContrasts("Treat-Control",
22levels=design)
23contrast.matrix##查看比较矩阵的信息，这里我们设置的是TreatVs.control
24##Contrasts
25##LevelsTreat-Control
26##Control-1
27##Treat1
28##拟合模型
29fit<-lmFit(eset_dat,design)
30fit2<-contrasts.fit(fit,contrast.matrix)
31fit2<-eBayes(fit2)
32DEG<-topTable(fit2,coef=1,n=Inf)%>%na.omit()##coef比较分组n基因数
33head(DEG)
34##logFCAveExprtP.Valueadj.P.ValB
35##ALDH3A1-3.22726310.302323-10.7103064.482850e-050.3134585-4.048355
36##CYP1B1-3.03368413.287607-10.5058884.995753e-050.3134585-4.049713
37##CYP1A1-9.00335311.481268-8.3714761.762905e-040.7374232-4.069681
38##HHLA2-1.5505876.595658-7.4434313.337672e-040.9308066-4.083411
39##SLC7A5-2.47033313.628775-7.2988683.708688e-040.9308066-4.085966
40##TIPARP-1.58127412.764218-7.0242524.552834e-040.9522252-4.091192
41dim(DEG)
42##[1]125496
43save(DEG,file="DEG_all.Rdata")
绘制火山图
火山图其实仅仅是一种可视化的方式，能够从整体上让我们对整体的差异分析情况有个了解
筛选到差异基因后，可以直接绘制出火山图
火山图的横坐标为logFC, 纵坐标为-log10(pvalue)，因此其实理论上讲plot即可完成火山图绘制
最简单原始的画法
1colnames(DEG)
2##[1]"logFC""AveExpr""t""P.Value""adj.P.Val""B"
3plot(DEG$logFC,-log10(DEG$P.Value))
Fig6高级的画法
借助于有人开发的更高级的包，用于完成某些特殊的功能，或者更美观
1require(EnhancedVolcano)
2EnhancedVolcano(DEG,
3
4lab=rownames(DEG),
5
6x="logFC",
7
8y="P.Value",
9
10selectLab=rownames(DEG)[1:5],
11
12xlab=bquote(~Log[2]~"foldchange"),
13
14ylab=bquote(~-Log[10]~italic(P)),
15
16pCutoff=0.05,##pvalue阈值
17
18FCcutoff=1,##FCcutoff
19
20xlim=c(-5,5),
21
22transcriptPointSize=1.8,
23
24transcriptLabSize=5.0,
25
26colAlpha=1,
27
28legend=c("NS","Log2FC","p-value",
29"p-value&Log2FC"),
30
31legendPosition="bottom",
32
33legendLabSize=10,
34
35legendIconSize=3.0)
Fig7绘制热图
热图的使用比较频繁，得到差异基因后可以直接绘制热图
相对简单好用的要属pheatmap包了
管道中的常规提取需要加上特殊的占位符.
1##首先提取出想要画的数据
2head(DEG)
3##提取FC前50
4up_50<-DEG%>%as_tibble()%>%
5mutate(genename=rownames(DEG))%>%
6dplyr::arrange(desc(logFC))%>%
7.$genename%>%.[1:50]##管道符中的提取
8##FC低前50
9down_50<-DEG%>%as_tibble()%>%
10mutate(genename=rownames(DEG))%>%
11dplyr::arrange(logFC)%>%
12.$genename%>%.[1:50]##管道符中的提取
13index<-c(up_50,down_50)
14
15##开始绘图-最简单的图
16library(pheatmap)
17pheatmap(eset_dat[index,],show_colnames=F,show_rownames=F)
Fig8稍微调整细节
1index_matrix<-t(scale(t(eset_dat[index,])))##归一化
2index_matrix[index_matrix>1]=1
3index_matrix[index_matrix<-1]=-1
4head(index_matrix)
5
6##添加注释
7anno=data.frame(group=group_list)
8rownames(anno)=colnames(index_matrix)
9anno##注释信息的数据框
10pheatmap(index_matrix,
11show_colnames=F,
12show_rownames=F,
13cluster_cols=T,
14annotation_col=anno)
15
Fig9
本期内容就到这里，我是白介素2，下期再见
作者：白介素2
如果觉得《GEO实操|limma分析差异基因》对你有帮助，请点赞、收藏，并留下你的观点哦！
