前面通過兩期視頻
給大家詳細(xì)介紹了如何從TCGA數(shù)據(jù)庫下載某一個腫瘤的RNAseq數(shù)據(jù)和miRNA-seq數(shù)據(jù)。如果看過視頻的小伙伴應(yīng)該都知道，TCGA里面不論是RNAseq數(shù)據(jù)還是miRNA-seq數(shù)據(jù)，每一個樣本的數(shù)據(jù)都是一個單獨(dú)的文件
并不像GEO里面下載得到的表達(dá)譜數(shù)據(jù)一般都是一個完整的表達(dá)譜矩陣，每一行是一個基因，每一列是一個樣本。
我們一般下游的分析，不論是差異表達(dá)分析，聚類分析，主成分分析，還是WGCNA分析，都是基于表達(dá)矩陣進(jìn)行的。所以要想利用TCGA的表達(dá)譜進(jìn)行下游數(shù)據(jù)分析，第一步就必須將這些單獨(dú)的文件都合并成一個類似于GEO那樣的表達(dá)矩陣。
下面先給大家捋一捋思路，我們以TCGA-CHOL（膽管癌）這套數(shù)據(jù)為例
- 合并RNAseq數(shù)據(jù)得到表達(dá)矩陣，RNAseq相關(guān)數(shù)據(jù)下載參考下面視頻
?如何從TCGA數(shù)據(jù)庫下載RNAseq數(shù)據(jù)以及臨床信息(一)?
1.1 通過讀取RNAseq_sample_sheet.csv文件，獲取每一個樣本表達(dá)譜數(shù)據(jù)的文件路徑。
Sample sheet文件中第一列是包裹RNA counts文件的文件夾，而第二列就是包含每個樣本RNA counts數(shù)的htseq_counts.gz文件了。需要將這兩列合并起來就可以得到counts文件的絕對路徑了，這樣代碼才能找到這個文件。
1.2 循環(huán)來讀取每一個文件里的內(nèi)容，這里每一個文件都是一個gz壓縮的文件，可以不用事先解壓，我們可以用R代碼來直接讀取壓縮文件里面的內(nèi)容。這個文件有兩列，第一列是gene的ensembl gene ID，第二列是比對到每一個gene上的counts數(shù)。
TCGA-CHOL這套數(shù)據(jù)一共有45個樣本，所以一共有45個counts文件。
1.3 將45個文件的第二列合并起來，第一列的ensembl gene id作為行名。Sample sheet文件中的sample ID作為列名。
1.4 去除ensembl gene id后面的點(diǎn)和數(shù)字。舉個例子ENSG00000000003.13變成ENSG00000000003。
1.5 去掉非基因的行。Counts文件拖到最后，你會發(fā)現(xiàn)有些行并不是基因的counts數(shù)，我們需要將這些行從最后的表達(dá)矩陣?yán)锩鎰h除。
最后我們就可以得到RNAseq的表達(dá)矩陣了
2. 合并mirRNA seq數(shù)據(jù)的到表達(dá)矩陣
整體思路跟合并RNAseq數(shù)據(jù)大同小異
miRNA seq數(shù)據(jù)下載參考下面視頻
?如何從TCGA數(shù)據(jù)庫下載miRNA數(shù)據(jù)(二)?
2.1 通過讀取RNAseq_sample_sheet.csv文件，獲取每一個樣本表達(dá)譜數(shù)據(jù)的文件路徑。
Sample sheet文件中第一列是包裹mirbase21.isoforms.quantification.txt文件的文件夾，而第二列就是包含每個樣本mir counts數(shù)的quantification文件了。需要將這兩列合并起來就可以得到mir counts文件的絕對路徑了，這樣代碼就可以找到這個文件。
2.2 循環(huán)來讀取每一個文件里的內(nèi)容，這里每一個文件就是普通的txt文本文件。這個文件跟RNA counts文件不太一樣。
我們知道有時(shí)候一個miRNA的前體，可以產(chǎn)生兩個miRNA成熟體，3‘端產(chǎn)生一個，5‘端產(chǎn)生一個。上圖就是一個例子，前體都是has-let-7a-1，但是成熟體有兩個MIMAT0000062（hsa-let-7a-5p）和MIMAT0004481（hsa-let-7a-3p）。通過在mirbase（http://www.mirbase.org/）數(shù)據(jù)庫里面查詢這個miRNA可以得到兩個成熟體的信息。
我們這里不應(yīng)該將所有的read_count都算到has-let-7a-1這個前體上，而是應(yīng)該分別算到兩個成熟體上。因?yàn)檫@兩個成熟體的靶基因一般是不一樣的，應(yīng)為他們的種子序列一般是不一樣的（種子序列概念參考：miRNA 靶向預(yù)測軟件targetscan）。TCGA-CHOL這套數(shù)據(jù)一共有45個樣本，所以一共有45個mirbase21.isoforms.quantification.txt文件。
2.3 將45個文件的所有miRNA成熟體的counts和并起來
2.4 我們把miRNA成熟體的ID（eg. MIMAT0000062）轉(zhuǎn)換成miRNA的名字來作為表達(dá)譜的行名。Sample sheet文件中的sample ID作為列名。做ID轉(zhuǎn)換的表可以從mirbase（http://www.mirbase.org/）下載得到，
如何處理可以參考下面的視頻
?如何獲取人的所有miRNA ID?
2.5 將樣本中不存在的miRNA的counts數(shù)設(shè)置為0
最后我們就可以得到RNAseq的表達(dá)矩陣了
大家可以根據(jù)這個思路自己去寫R代碼。當(dāng)然小編也已經(jīng)為大家寫好了R代碼, ?點(diǎn)擊獲取本文中用到的R代碼?。
下面我們演示一下怎么用這段代碼來合并TCGA-CHOL的RNAseq數(shù)據(jù)和miRNA數(shù)據(jù)，得到mRNA和miRNA的表達(dá)矩陣。
?R代碼合并TCGA表達(dá)譜數(shù)據(jù)?
歡迎關(guān)注公眾號 “生信交流平臺”，我們創(chuàng)建了交流群，供大家交流討論生信問題。
參考資料：
