转录组GO富集与微生物相关性分析
转录组GO富集与微生物相关性分析
原始数据格式
我相信做生物分析的都知道很多时候难点往往不在于代码怎么写，而在于数据格式很多时候与各种包的要求不相符合，因此我先列一下本次分析中用到的数据格式
表1：转录组数据使用原始矩阵
| A_1 | A_2 | A_3 | B_1 | B_2 | B_3 | C_1 | C_2 | C_3 | D_1 | …… | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Gene0000001 | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… |
| Gene0000002 | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… |
| Gene0000003 | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… |
| Gene0000004 | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… |
| …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… | …… |
表2：分组数据
| sample_ID | group |
|---|---|
| A_1 | A |
| A_2 | A |
| A_3 | A |
| B_1 | B |
| B_2 | B |
| B_2 | B |
| C_1 | C |
| …… | …… |
表3：细菌丰度表
| genus1 | genus2 | genus3 | genus4 | …… | |
|---|---|---|---|---|---|
| A_1 | …… | …… | …… | …… | …… |
| A_2 | …… | …… | …… | …… | …… |
| A_3 | …… | …… | …… | …… | …… |
| B_1 | …… | …… | …… | …… | …… |
| B_2 | …… | …… | …… | …… | …… |
| …… | …… | …… | …… | …… | …… |
表4：基因注释文件
由于我做的是非模式生物，所以进行GO富集时还需要一个基因注释文件，这个文件的后缀是.map（用文本文档按下列格式弄好改后缀就行），gene_id列与GO_number列之间用\t制表符分隔
| gene_id | GO_number |
|---|---|
| Gene0000001 | GO:0005515 |
| Gene0000001 | GO:0051259, GO:0046872, GO:0030176, GO:0005783, GO:0030968, GO:0016874, GO:0016021, GO:0016020, GO:0008270, GO:0000299, GO:0006928, GO:0006512, GO:0006511, GO:0007165, GO:0005515, GO:0030433, GO:0005789, GO:0004842, GO:0000209, GO:0004872 |
| Gene0000003 | GO:0008076, GO:0005242, GO:0030955, GO:0005244 |
| …… | …… |
使用TCC包进行差异基因分析
TCC是什么包？为什么要用TCC包？目前常用的差异基因分析包为DEseq2或者edgeR,但是这两个包在寻找差异基因时只能够找出两组之间的差异基因，即使你的数据里有A、B、C、D四个组，DESeq2分析后也只会返回A vs B，A vs C，A vs D，B vs C……的结果，而不会找出在四个组中有显著性差异的基因。这对于做临床研究、药理研究的同志们来说可能够了，但是对于植物学方面的就不够用了（不同器官、不同生长阶段等）。那么有没有可以多组同时比较的R包呢？经过检索我发现了这个TCC包可以满足要求1。其实这个包是在DESeq，DESeq2和edgeR三个包基础上构建了一个适用于多组分析的pipeline。这里我使用的是基于DESeq2的管道，更多内容可以参考帮助文档。分析代码如下：
#gene_count_matrix，为基因表达矩阵，读入过程略
#group为分组，读入过程略
library(TCC)
#构建TCC类
tcc <- new("TCC", gene_count_matrix, group)
#下面使用TCC包基于DESeq2的管道方案
ssss <- calcNormFactors(tcc, norm.method = "deseq2",
test.method = "deseq2",
iteration = 3)
DE <- estimateDE(ssss,
test.method = "deseq2",
full = ~ group,
reduced = ~ 1)
res <- getResult(DE)
res就是分析后的结果了，根据names可以很清楚这个结果包含哪些信息
使用topGO包进行GO富集分析
这里没什么说的，就是一个典型的topGO分析，有不明白的可以参考其他博主的内容
#准备基因列表，注意，一定要因子化！！！
genelist<-res$estimatedDEG
names(genelist)<-res$gene_id
genelist<-as.factor(genelist)
#读入基因注释文件
gene2GO_map<-readMappings("D:/gene2GO.map")
#构建topGO类，注意MF。CC和BP是分开构建三个。
Godata_BP <- new(
"topGOdata",
ontology = "BP",
allGenes = genelist,
annot = annFUN.gene2GO,
gene2GO = gene2GO_map)
#寻找差异GO
resultFis_BP <- runTest(Godata_BP, algorithm = "classic", statistic = "fisher")
#找出排名前100的GO
sig.tab <- GenTable(
Godata_BP,
Fis = resultFis_BP,
topNodes = 100)
计算转录组与微生物组相关性
这边由于微生物和转录组在两个不同的矩阵中，所以自己写了个函数进行相关性计算，比较懒直接用了循环嵌套，速度有些慢，用apply族写的话速度应该会快很多。需要注意的是，循环数小的应该放在外层，循环数大的应该放在内层
get_matrix_cor<-function(x,y)
{
cormatrix<-as.data.frame(matrix(numeric(0),ncol=length(x),nrow = length(y)))
for (i in (1:length(x))){
for (j in (1:length(y))){
cormatrix[j,i]<-cor(x[,i],y[,j])
}
}
}
计算相关性矩阵
trans_bac_cor_matrix<-get_matrix_cor(bac_genus_matrix，t(gene_count_matrix)) %>% as.data.frame()
根据GO号提取相关性矩阵并绘图
一个转录组数据里面的gene数以万计，画在一张图上就什么都看不见了，因此根据我感兴趣的GO号提取出来画图会好的多。GO号的选取可以根据富集结果的sig.tab选取，这里重点讲述画图和输出结果的过程
提取相关性矩阵
#输入一个GO号
GO="0005985"
#获取感兴趣的GO号相对应的gene名字
a<-genesInTerm(Godata_BP,paste0("GO:",GO))[[1]] %>% as.vector()
#从相关系数矩阵中提取对应基因的行，建议用apply族函数改写，简化代码
b<-c()
for (i in 1:length(a) ){
b[i]<-which(row.names(trans_bac_cor_matrix)==a[i])
}
matrix_BF<-trans_bac_cor_matrix[b,]
这里的matrix_BF就是转录组与微生物组的相关性矩阵了，后期用来进行热图绘制。但是用pheatmap包进行热图绘制前需要先去除NA
#删除NA
matrix_BF<-matrix_BF[-which(is.na(matrix_BF)),]
na_flag <- apply(is.na(matrix_BF), 2, sum)
matrix_BF<-matrix_BF[,na_flag==0]
#注意，我这里删除NA的策略是将某种存在NA的菌直接删去，所以如果有某个基因全是
#NA的话，那matrix_BF就会成为一个空数据框，这个时候需要将这行删去。
绘制热图
#绘制热图
annotation_col = data.frame(
ClassBac = factor(paste0('Cluster',cutree(hclust_BF,10)))
)
rownames(annotation_col) = colnames(matrix_BF)
pdf(paste0("D:/GO_",GO,".pdf"))
heatmap_BF<-pheatmap::pheatmap(matrix_BF,show_colnames = F,annotation_col = annotation_col)
dev.off()
这里展示一张（涉及数据因此基因名盖住了）
输出与热图对应的相关性矩阵
pheatmap会自动聚类，这样绘制完成的热图就与之前提取的相关性矩阵不对应了，因此在这里进行一个处理，使得我们输出的矩阵与热图顺序一致。该部分内容参考了一位博主的方法，但是后来找不到页面了，如有侵权，与我联系
hclust_BF=hclust(dist(t(matrix_BF)))
hclust_BF_row=hclust(dist(matrix_BF))
col_cluster=cutree(hclust_BF,k=10)
row_cluster=cutree(hclust_BF_row,k=2)
newOrder=matrix_BF[hclust_BF_row$order,]
newOrder[,ncol(newOrder)+1]=row_cluster[match(row.names(newOrder),names(row_cluster))]
names(newOrder)[ncol(newOrder)]="Cluster"
names2=names(matrix_BF[,hclust_BF$order])
newOrder[nrow(newOrder)+1,]=col_cluster[match(c(names2,"cluster"),names(col_cluster))]
row.names(newOrder)[nrow(newOrder)]="Cluster"
#输出对应矩阵
write.csv(newOrder,paste0("D:/GO_",GO,".csv"))
这样输出的矩阵就与热图上的位置一致了。