原始資料格式

我相信做生物分析的都知道很多時候難點往往不在於程式碼怎麼寫，而在於資料格式很多時候與各種包的要求不相符合，因此我先列一下本次分析中用到的資料格式

表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包可以滿足要求¹。其實這個包是在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"))

這樣輸出的矩陣就與熱圖上的位置一致了。

---

1. http://www.bioconductor.org/packages/release/bioc/html/TCC.html ↩︎