GTEx數據庫-TCGA數據挖掘的好幫手

通常我們在挖掘TCGA數據庫的時候，會發現該項目納入的正常組織測序結果是非常少的，也就是說很多病人都不會有他的正常組織的轉錄組測序結果，比如說乳腺癌吧，1200個左右的轉錄組數據，其中1100左右都是腫瘤組織的測序數據，只有區區100個左右的正常對照。

這個時候我們就需要想辦法加大正常組織測序樣本量，既然TCGA數據庫沒有，我們就從其他數據庫着手。

這里值得大力推薦的是GTEx數據庫 ,Genotype-Tissue Expression (GTEx)

背景知識

一期

2015年，GTEx發布了第一個階段性成果，一次性在Science雜志上發表三篇研究成果，該成果還被選為封面文章。GTEx的研究從175名死者身上采集到了1641個屍檢樣本，這些樣本來自54個不同的身體部位，對幾乎所有轉錄基因的基因表達模式進行了觀察，從而夠確定基因組中影響基因表達的特定區域。另外兩篇文章之一從人所有組織中的基因表達譜進行了描述，證明了組織特異性的某些基因往往決定了組織特異性基因的表達調控；另一篇解釋了截短的蛋白變異體如何影響組織中的基因表達。

二期

在2017年，一次性在nature發表4篇研究成果，GTEx研究聯盟的研究收集並研究了來自449名生前健康的人類捐獻者的7000多份屍檢樣本，涵蓋44個組織（42種不同的組織類型），包括31個實體器官組織、10個腦分區、全血、兩個來自捐獻者血液和皮膚的細胞系，作者利用這些樣本研究基因表達在不同組織和個體中有何差異。題為“Landscape of X chromosome inactivation across human tissues”和“Dynamic landscape and regulation of RNA editing in mammals”的論文，采用GTEx數據探討了與基因表達相關聯的基因變異如何能夠調節RNA編輯和X染色體失活現象。

數據庫內容介紹

通常是直接去 https://gtexportal.org/ 找到可以下載的數據集，如下：

其中，對我們來說最重要的就是表達矩陣，可以下載圖中 gene read counts 這個496M的文件，表達矩陣里面的樣本ID肯定是數據庫組織者自定義的，所以我們還需要找到樣本ID的注釋信息。

更多的是關於這個數據庫的網頁使用介紹，我們生信工程師通常不需要，就不贅述了。

注意一下數據庫的版本信息：

The current release is V7 including 11,688 samples, 53 tissues and 714 donors

首先看數據庫的注釋信息

重點是：

 # SMTS Tissue Type, area from which the tissue sample was taken. 
 # SMTSD Tissue Type, more specific detail of tissue type

可以看到每個樣本屬於哪一種組織，這樣就方便提取他們的信息來輔助自己的研究。

把 gene read counts 這個496M的表達矩陣導入R：

if(F){ options(stringsAsFactors = F) GTEx=read.table('~/Downloads/GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_reads.gct.gz' ,header = T,sep = '\t',skip = 2) GTEx[1:4,1:4] h=head(GTEx) save(h,file = 'GTEx_head.Rdata') }

挑選感興趣的組織的表達矩陣

上面我們詳細了解了不同樣本注釋到的組織，所以代碼很簡單：

 load('~/Desktop/GTEx_all.Rdata') a[1:4,1:4] colnames(a) # SMTS Tissue Type, area from which the tissue sample was taken. # SMTSD Tissue Type, more specific detail of tissue type b=read.table('GTEx_v7_Annotations_SampleAttributesDS.txt', header = T,sep = '\t',quote = '') table(b$SMTS) breat_gtex=a[,gsub('[.]','-',colnames(a)) %in% b[b$SMTS=='Breast',1]] rownames(breat_gtex)=a[,1] dat=breat_gtex

就是把屬於breast這個組織的樣本名挑選出來，在上面的表達矩陣里面取子集即可。

值得注意的是這個時候的表達矩陣基因名不是symbol，是需要進行ID轉換的，代碼如下：

dat=breat_gtex
 ids=a[,1:2]
 head(ids)
 colnames(ids)=c('probe_id','symbol')
 dat=dat[ids$probe_id,]
 dat[1:4,1:4] 
 ids$median=apply(dat,1,median)
 ids=ids[order(ids$symbol,ids$median,decreasing = T),]
 ids=ids[!duplicated(ids$symbol),]
 dat=dat[ids$probe_id,]
 rownames(dat)=ids$symbol
 dat[1:4,1:4] 
 breat_gtex=dat
 save(breat_gtex,file = 'breat_gtex_counts.Rdata')

表達矩陣如下所示：

正常乳腺組織樣本表達矩陣可以進行的分析

通常情況下應該是去和腫瘤數據進行分析，那樣的分析就多元化了，這里來個簡單點的，可以進行pam50分類：

if(T){ ddata=t(dat) ddata[1:4,1:4] s=colnames(ddata);head(s) library(org.Hs.eg.db) s2g=toTable(org.Hs.egSYMBOL) g=s2g[match(s,s2g$symbol),1];head(g) # probe Gene.symbol Gene.ID dannot=data.frame(probe=s, "Gene.Symbol" =s, "EntrezGene.ID"=g) ddata=ddata[,!is.na(dannot$EntrezGene.ID)] dannot=dannot[!is.na(dannot$EntrezGene.ID),] head(dannot) library(genefu) # c("scmgene", "scmod1", "scmod2","pam50", "ssp2006", "ssp2003", "intClust", "AIMS","claudinLow") s<-molecular.subtyping(sbt.model = "pam50",data=ddata, annot=dannot,do.mapping=TRUE) table(s$subtype) tmp=as.data.frame(s$subtype) subtypes=as.character(s$subtype) } library(genefu) pam50genes=pam50$centroids.map[c(1,3)] pam50genes[pam50genes$probe=='CDCA1',1]='NUF2' pam50genes[pam50genes$probe=='KNTC2',1]='NDC80' pam50genes[pam50genes$probe=='ORC6L',1]='ORC6' x=dat x=x[pam50genes$probe[pam50genes$probe %in% rownames(x)] ,] tmp=data.frame(subtypes=subtypes) rownames(tmp)=colnames(x) library(pheatmap) pheatmap(x,show_rownames = T,show_colnames = F, annotation_col = tmp, filename = 'ht_by_pam50_raw.png') x=t(scale(t(x))) x[x>1.6]=1.6 x[x< -1.6]= -1.6 pheatmap(x,show_rownames = T,show_colnames = F, annotation_col = tmp, filename = 'ht_by_pam50_scale.png')

單獨取出pam50包含的50個基因的表達矩陣進行熱圖聚類：

由上圖可以看到不同基因的表達量是差異很大的，通常我們不會去比較不同基因的表達量，而只是比較同一個基因在不同樣本的表達量差異的。

所以我們沒有必要去看不同基因的表達量高低，那么就可以進行一定程度的歸一化，重新繪圖如下：

可以很明顯的看到哪怕是對正常組織的轉錄組測序結果走pam50的分類也是可以拿到各種各樣的分類結果的。

但是pam50的分類是在乳腺癌患者的芯片表達矩陣進行訓練的模型，是因為我們用錯了地方，可以看看在METEBRIC里面的分類結果：

上面的分類是pam50算法的結果，下面的分類是臨床信息。

可以看到basal的結果還是很統一的，而且都比較符合TNBC的定義，就是PGR,ESR1,ERBB2都表達量都很低。

去除批次效應

如果真的要把GTEx數據庫的轉錄組表達矩陣和TCGA的進行比較，還需要一定程度的去除批次效應。

我以前在生信技能樹多次講解，這里也不再贅述。