據(jù)格式到批次效應(yīng)校正)
做TCGA數(shù)據(jù)挖掘很多人的精力都花在差異分析本身等跑完一看結(jié)果要么差異基因數(shù)少得可憐要么多到?jīng)]法解釋回頭排查才發(fā)現(xiàn)是差異分析前處理這步出了問(wèn)題。這篇文章把TCGA數(shù)據(jù)處理的完整前處理流程拆開(kāi)講清楚從數(shù)據(jù)格式辨析、下載渠道選擇、樣本分組、表達(dá)量矩陣構(gòu)建到低表達(dá)基因過(guò)濾、批次效應(yīng)校正、標(biāo)準(zhǔn)化轉(zhuǎn)換每個(gè)環(huán)節(jié)都說(shuō)清做什么、為什么這么做、以及我踩過(guò)的坑。適合正在做TCGA數(shù)據(jù)挖掘、尤其是不太清楚數(shù)據(jù)下載之后該怎么下手的生信新手參考。1. 先搞懂TCGA數(shù)據(jù)的三種形態(tài)不然后面全是坑TCGA全稱(chēng)是The Cancer Genome Atlas癌癥基因組圖譜這是由一個(gè)國(guó)家級(jí)大規(guī)模測(cè)序項(xiàng)目積累下來(lái)的多組學(xué)數(shù)據(jù)庫(kù)覆蓋33種癌癥類(lèi)型、超過(guò)2萬(wàn)例樣本包含了轉(zhuǎn)錄組、拷貝數(shù)變異、甲基化、突變等多層面數(shù)據(jù)。做癌癥相關(guān)生信分析的人基本繞不開(kāi)這個(gè)庫(kù)。但很多人下載完轉(zhuǎn)錄組表達(dá)數(shù)據(jù)后一打開(kāi)文件就懵了因?yàn)橥瑯右粋€(gè)TCGA項(xiàng)目你能下載到好幾種不同格式的表達(dá)量數(shù)據(jù)它們之間的數(shù)值含義完全不同用錯(cuò)了地方會(huì)讓后續(xù)所有分析都失去意義。1.1 Counts、FPKM、TPM到底有什么區(qū)別TCGA的轉(zhuǎn)錄組表達(dá)數(shù)據(jù)在GDC官方門(mén)戶(hù)和第三方平臺(tái)上最常見(jiàn)的就是三種格式HTSeq-Counts、FPKM/FPKM-UQ、TPM。簡(jiǎn)單說(shuō)HTSeq-Counts是基因比對(duì)后統(tǒng)計(jì)到的原始read計(jì)數(shù)它沒(méi)有做任何文庫(kù)大小和基因長(zhǎng)度的校正數(shù)值大小直接受測(cè)序深度影響同一批樣本之間如果測(cè)序量差異大raw counts就不適合直接拿來(lái)比較。FPKMFragments Per Kilobase Million是RNA-seq早期常用的標(biāo)準(zhǔn)化指標(biāo)計(jì)算時(shí)先按基因長(zhǎng)度歸一化再按文庫(kù)大小歸一化。它和TPM看起來(lái)很像但計(jì)算順序不同導(dǎo)致二者在數(shù)值分布上有細(xì)微差別。FPKM在跨樣本比較時(shí)其實(shí)存在一點(diǎn)系統(tǒng)偏差而TPM在實(shí)現(xiàn)上做了修正單位是“每百萬(wàn)條轉(zhuǎn)錄本中來(lái)自該基因的轉(zhuǎn)錄本數(shù)”樣本間可比性更好所以現(xiàn)在很多平臺(tái)默認(rèn)給TPM。那差異分析前處理到底該選哪個(gè)這取決于你后面用什么工具。如果你打算用DESeq2或edgeR跑標(biāo)準(zhǔn)差異分析官方推薦的輸入是raw count矩陣HTSeq-Counts因?yàn)檫@兩個(gè)工具內(nèi)部有自己的一套歸一化邏輯你給它的應(yīng)該是“沒(méi)被過(guò)度處理過(guò)的原始計(jì)數(shù)”。如果你只是想做基于表達(dá)量高低的分組比較、或者直接對(duì)表達(dá)值做t檢驗(yàn)相關(guān)的差異篩選那用TPM然后log2轉(zhuǎn)換是比較常見(jiàn)的做法。需要特別提醒絕對(duì)不要把FPKM或者TPM塞給DESeq2結(jié)果會(huì)非常難看差異基因數(shù)量動(dòng)不動(dòng)就成千上萬(wàn)明顯失真。1.2 數(shù)據(jù)下載GDC官方和UCSC Xena怎么選下載TCGA表達(dá)數(shù)據(jù)主要有兩條主流路徑。第一條是GDC Data Portalportal.gdc.cancer.gov這是官方數(shù)據(jù)倉(cāng)庫(kù)能拿到最原始的HTSeq-Counts文件下載時(shí)可以選擇per-sample的單獨(dú)文件也可以通過(guò)GDC API或者像TCGAbiolinks這樣的R包批量下載適合對(duì)數(shù)據(jù)溯源要求嚴(yán)格的課題。官方渠道的數(shù)據(jù)更新及時(shí)參考基因組版本和注釋信息都標(biāo)得很清楚寫(xiě)方法部分時(shí)引起來(lái)也方便。第二條是UCSC Xenaxenabrowser.net這個(gè)平臺(tái)把TCGA的數(shù)據(jù)做了統(tǒng)一整理可以按需求下載整理好的表達(dá)矩陣提供TPM和FPKM兩種格式還附帶臨床信息。用Xena最舒服的一點(diǎn)是它直接給你一個(gè)基因x樣本的矩陣不需要自己逐樣本合并對(duì)生信基礎(chǔ)不太強(qiáng)的人友好很多。它的數(shù)據(jù)也經(jīng)過(guò)了統(tǒng)一流程處理樣本名、基因名格式都比較規(guī)范。兩套方案怎么選我的建議是如果后續(xù)要用DESeq2/edgeR就走GDC拿raw count如果只做簡(jiǎn)單的表達(dá)量差異篩選Xena的TPM矩陣完全夠用。另外還有一個(gè)常見(jiàn)渠道是cBioPortal但它的表達(dá)數(shù)據(jù)多是處理過(guò)的適合查詢(xún)和可視化不太適合拿來(lái)當(dāng)差異分析的原始輸入。實(shí)際處理TCGA數(shù)據(jù)時(shí)經(jīng)常是下載完才發(fā)現(xiàn)格式不對(duì)所以第一步先把數(shù)據(jù)形態(tài)和后續(xù)分析工具對(duì)應(yīng)好能省下一大堆返工時(shí)間。2. 樣本分組前處理的第一個(gè)分水嶺拿到表達(dá)數(shù)據(jù)后的第一件事不是急著做差異分析而是把樣本分組搞清楚。TCGA的樣本編號(hào)不是隨便起的里面藏著樣本類(lèi)型、組織來(lái)源、是否配對(duì)等信息如果分組出錯(cuò)后面所有差異比較都是白做。2.1 腫瘤和正常樣本到底怎么區(qū)分TCGA樣本編號(hào)是有一套規(guī)則的四段式結(jié)構(gòu)比如TCGA-XX-XXXX-01A-11R-XXXX-XX其中第四段的開(kāi)頭兩位數(shù)字代表樣本類(lèi)型。01開(kāi)頭的是原發(fā)性實(shí)體瘤11開(kāi)頭的是癌旁正常組織solid tissue normal05是原發(fā)性血液腫瘤06是轉(zhuǎn)移性腫瘤還有02復(fù)發(fā)性腫瘤、03原發(fā)性血液腫瘤等少見(jiàn)類(lèi)型。日常分析中最常見(jiàn)的就是01和11兩類(lèi)。這里面最容易被忽略的一點(diǎn)是不是所有癌癥類(lèi)型都同時(shí)有01和11兩類(lèi)樣本。像GBM膠質(zhì)母細(xì)胞瘤、LGG腦低級(jí)別膠質(zhì)瘤這類(lèi)中樞神經(jīng)系統(tǒng)腫瘤以及一些罕見(jiàn)癌種正常對(duì)照樣本數(shù)量極少甚至沒(méi)有。這時(shí)候如果你還想著用TCGA自身樣本做配對(duì)差異分析樣本量就會(huì)非常尷尬。常見(jiàn)的替代方案是去GTEx數(shù)據(jù)庫(kù)拉正常組織數(shù)據(jù)來(lái)補(bǔ)充對(duì)照不過(guò)這又會(huì)引入跨數(shù)據(jù)集的批次效應(yīng)問(wèn)題處理起來(lái)要另外花心思。在寫(xiě)分組文件通常叫sample_info或者phenotype表的時(shí)候有幾個(gè)細(xì)節(jié)需要注意第一列名要簡(jiǎn)潔統(tǒng)一一般就是sample_id和group兩列第二組別命名建議用明確的class標(biāo)簽比如Tumor和Normal而不是01和11這種數(shù)字代號(hào)后面可視化的時(shí)候圖例會(huì)更清晰第三一定記得檢查有沒(méi)有重復(fù)的樣本ID有些樣本是多組織部位取樣如果你下載的時(shí)候沒(méi)嚴(yán)格過(guò)濾會(huì)出現(xiàn)一個(gè)病人對(duì)應(yīng)多條記錄的情況。2.2 表達(dá)量矩陣構(gòu)建和基因ID轉(zhuǎn)換如果你的數(shù)據(jù)是從GDC下載的per-sample文件每個(gè)樣本是一個(gè)獨(dú)立的counts文件那第一步就是把所有文件合并成一個(gè)矩陣。合并的邏輯很簡(jiǎn)單讀入所有文件把里面每行基因的counts值一一對(duì)應(yīng)到矩陣的列上。但實(shí)際操作有幾個(gè)很容易出問(wèn)題的坑。第一個(gè)坑是gene_id的格式。GDC下載的HTSeq-Counts文件里行名通常是Ensembl基因ID的帶版本格式比如ENSG00000242268.11。你要用sub函數(shù)把小數(shù)點(diǎn)和后面的版本號(hào)去掉再去做ID注釋轉(zhuǎn)換。對(duì)應(yīng)的基因symbol轉(zhuǎn)換最常用的是用biomaRt包實(shí)時(shí)查詢(xún)或者用org.Hs.eg.db包做本地注釋。我自己的經(jīng)驗(yàn)是網(wǎng)絡(luò)不穩(wěn)定時(shí)biomaRt很容易超時(shí)預(yù)先把Ensembl ID和symbol的映射關(guān)系保存成一份本地文件離線環(huán)境下也能直接轉(zhuǎn)換效率高很多。第二個(gè)坑是重復(fù)基因名的處理。Ensembl ID轉(zhuǎn)換到symbol之后大概率會(huì)出現(xiàn)多個(gè)ID對(duì)應(yīng)同一個(gè)symbol的情況也就是基因名重復(fù)。這時(shí)候不能直接保留需要按累積表達(dá)量或者最大表達(dá)量合并去重否則后面差異分析會(huì)報(bào)錯(cuò)或者產(chǎn)生重復(fù)行。還有一個(gè)細(xì)節(jié)是染色體上的小RNA、假基因等非蛋白編碼基因要不要過(guò)濾如果是常規(guī)mRNA表達(dá)譜差異分析我一般會(huì)保留protein_coding基因用注釋文件篩掉非編碼轉(zhuǎn)錄本能在源頭減少不少噪音。3. 低表達(dá)基因過(guò)濾這一步做不做結(jié)果差非常多處理完ID轉(zhuǎn)換和矩陣構(gòu)建之后接下來(lái)就是過(guò)濾。別小看這一步它對(duì)后續(xù)差異分析的影響非常大但也是最容易被新手跳過(guò)的一步。3.1 為什么不能把全基因集直接拿去做差異分析TCGA的HTSeq-Counts矩陣動(dòng)輒六萬(wàn)行但里面真正在樣本中穩(wěn)定表達(dá)的基因其實(shí)遠(yuǎn)沒(méi)有那么多。很多基因在絕大多數(shù)樣本里表達(dá)量為0或者只有幾個(gè)read這些基因就是噪音。你如果不做過(guò)濾直接跑DESeq2會(huì)產(chǎn)生兩個(gè)問(wèn)題一是多重檢驗(yàn)校正時(shí)需要比較的基因數(shù)量大幅膨脹padj會(huì)變嚴(yán)格一些真陽(yáng)性可能被壓掉二是大量全零或近零的基因會(huì)干擾后續(xù)離散度估計(jì)影響差異檢驗(yàn)的穩(wěn)健性。過(guò)濾的基本思路是設(shè)定一個(gè)保留標(biāo)準(zhǔn)比如“至少在20%的樣本中counts大于等于10”這類(lèi)規(guī)則。實(shí)際操作里DESeq2官方文檔建議的是一個(gè)簡(jiǎn)單但有據(jù)可依的過(guò)濾方式先對(duì)所有基因計(jì)算該基因在各樣本中的平均表達(dá)量或最大表達(dá)量再設(shè)一個(gè)閾值。edgeR則提供了filterByExpr函數(shù)可以綜合考慮最小計(jì)數(shù)、樣本量、分組信息來(lái)自動(dòng)推薦過(guò)濾條件。如果你用的是TPM矩陣做后續(xù)分析通常會(huì)把閾值設(shè)成TPM 1且在多少比例樣本中滿(mǎn)足條件。3.2 過(guò)濾閾值怎么定三個(gè)場(chǎng)景對(duì)比我整理了幾種常用的過(guò)濾標(biāo)準(zhǔn)你可以根據(jù)數(shù)據(jù)量和平時(shí)的習(xí)慣來(lái)選擇。過(guò)濾方式規(guī)則示例適用場(chǎng)景注意事項(xiàng)總量閾值型行和 10 或 mean counts 5全轉(zhuǎn)錄組初步篩選速度快偏寬松可能留下較多低表達(dá)基因比例閾值型至少在90%樣本中 counts 10關(guān)注優(yōu)勢(shì)表達(dá)基因的課題嚴(yán)格容易誤刪條件特異性表達(dá)基因工具推薦型filterByExpr自動(dòng)判斷常規(guī)差異分析準(zhǔn)備綜合庫(kù)大小和分組信息個(gè)人最常用第三種是edgeR推薦型底層會(huì)綜合所有樣本的測(cè)序深度和庫(kù)大小來(lái)設(shè)置閾值我個(gè)人用得最多因?yàn)樗诒Ａ粽鎸?shí)信號(hào)和去掉噪音之間平衡得比較好。順便說(shuō)一句這里有一個(gè)容易混淆的概念過(guò)濾和標(biāo)準(zhǔn)化誰(shuí)先誰(shuí)后的問(wèn)題。通常先過(guò)濾低表達(dá)基因再做標(biāo)準(zhǔn)化/歸一化邏輯上更順。原因是低表達(dá)基因的存在會(huì)影響某些標(biāo)準(zhǔn)化方法對(duì)文庫(kù)大小的估計(jì)先把確定是噪音的行去掉標(biāo)準(zhǔn)化會(huì)更穩(wěn)定。如果你用的是DESeq2它的median-of-ratios因素估計(jì)本身對(duì)低表達(dá)基因也敏感所以提前過(guò)濾是防患于未然。4. 批次效應(yīng)處理和數(shù)據(jù)標(biāo)準(zhǔn)化批次效應(yīng)是TCGA數(shù)據(jù)處理里最讓人頭疼的問(wèn)題之一但也是差異分析前必須面對(duì)的一關(guān)。TCGA的樣本不是一天之內(nèi)測(cè)完的樣本來(lái)源遍布多個(gè)組織中心、測(cè)序平臺(tái)和批次這些技術(shù)差異如果混進(jìn)你的分析里得到的結(jié)果很可能不是生物學(xué)差異而是技術(shù)噪音。4.1 批次效應(yīng)怎么發(fā)現(xiàn)主成分分析和聚類(lèi)熱圖批次效應(yīng)是指樣本在測(cè)序批次、文庫(kù)制備、芯片或平臺(tái)不同等因素影響下產(chǎn)生的系統(tǒng)性差異它和真實(shí)的生物學(xué)差異混雜在一起輕則讓PCA圖上的Tumor和Normal分不開(kāi)重則直接讓差異分析結(jié)果不可信。很多資料把批次效應(yīng)放在差分析之后才檢查但我建議在處理階段就提前看一遍免得后續(xù)返工。最直觀的方法是PCA。對(duì)log2標(biāo)準(zhǔn)化后的表達(dá)矩陣做PCA然后按樣本的分組信息、測(cè)序平臺(tái)如Illumina GA vs HiSeq、或者來(lái)源組織中心給點(diǎn)著色觀察樣本是否明顯按非生物學(xué)因素聚類(lèi)。如果發(fā)現(xiàn)按批次聚類(lèi)的現(xiàn)象很明顯就該考慮校正了。還有一個(gè)輔助方法是畫(huà)樣本相關(guān)性的熱圖如果同批次的樣本聚成了清晰的模塊而模塊內(nèi)既有Tumor也有Normal那基本可以判斷存在批次效應(yīng)。4.2 ComBat-seq和limma的removeBatchEffect怎么選處理批次效應(yīng)有不少工具最常用的是sva包的ComBat系列。ComBat適合處理微陣列和RNA-seq的表達(dá)矩陣近似連續(xù)數(shù)值ComBat-seq則是針對(duì)RNA-seq raw count專(zhuān)門(mén)開(kāi)發(fā)的它不改變數(shù)據(jù)的整數(shù)特性輸出結(jié)果可以繼續(xù)喂給DESeq2或edgeR做差異分析。如果你已經(jīng)把數(shù)據(jù)log2轉(zhuǎn)換成了連續(xù)值也可以用limma包的removeBatchEffect但它適合在標(biāo)準(zhǔn)化之后、差異分析之前對(duì)表達(dá)矩陣做殘差化處理。我的組合拳實(shí)踐經(jīng)驗(yàn)是先用filterByExpr過(guò)濾掉低表達(dá)基因再用ComBat_seq對(duì)raw counts校正批次校正完再跑DESeq2比較穩(wěn)健。如果你拿的是Xena的TPM矩陣TPM是連續(xù)值就log2(TPM1)之后用removeBatchEffect然后基于殘差矩陣做后續(xù)分析。這里有個(gè)細(xì)節(jié)批次信息最好是樣本的真實(shí)測(cè)序批次plate、seq center、tissue source site不要簡(jiǎn)單用“下載日期”替代因?yàn)橄螺d日期跟生物學(xué)變量完全無(wú)關(guān)反而可能引入新的混淆。4.3 log2轉(zhuǎn)換和標(biāo)準(zhǔn)化方法選擇的經(jīng)驗(yàn)log2轉(zhuǎn)換幾乎是TCGA表達(dá)數(shù)據(jù)可視化和差異篩選的標(biāo)配操作但要注意兩個(gè)容易被忽略的點(diǎn)。第一log2(x1)和log2(CPM1)的區(qū)別前者針對(duì)raw count或TPM后者是先將counts轉(zhuǎn)成CPM再做log2二者數(shù)值分布不同下游算法對(duì)輸入類(lèi)型敏感建議從頭到尾保持一致。第二log2轉(zhuǎn)換只適合方差穩(wěn)定的數(shù)據(jù)場(chǎng)景如果你要跑的是方差依賴(lài)的統(tǒng)計(jì)模型比如DESeq2負(fù)二項(xiàng)模型千萬(wàn)不能自己先log2再喂進(jìn)去應(yīng)該把raw counts原樣交給DESeq2處理。至于TPM和CPM的選擇我個(gè)人觀點(diǎn)是TCGA的TPM矩陣已經(jīng)考慮到基因長(zhǎng)度影響適合做表達(dá)定量但如果你的分析目標(biāo)是比較同一樣本內(nèi)部基因間的表達(dá)水平TPM比CPM合適。實(shí)際處理中從一個(gè)矩陣出發(fā)先明確下游分析類(lèi)型再?zèng)Q定采用哪種數(shù)據(jù)形態(tài)不要事到臨頭才來(lái)回切換切換過(guò)程中數(shù)值分布的變化很可能把你的差異分析結(jié)果帶偏。5. 實(shí)操?gòu)脑嘉募娇芍苯幼霾町惙治龅臄?shù)據(jù)理論說(shuō)了一大堆我放一套自己常用的R腳本流程出來(lái)。這個(gè)流程從GDC下載的per-sample HTSeq-Counts文件入手最終得到可以直接做差異分析的數(shù)據(jù)結(jié)構(gòu)。你只需準(zhǔn)備好兩個(gè)目錄一個(gè)放所有樣本的counts文件每個(gè)文件兩列g(shù)ene_id和count一個(gè)放樣本注釋表格包含sample_id、group和batch等列。5.1 準(zhǔn)備工作與讀取數(shù)據(jù)library(DESeq2) library(edgeR) library(sva) library(biomaRt) library(dplyr) library(tibble) # 讀入所有樣本的counts文件 files - list.files(path counts_dir, pattern *.txt, full.names TRUE) sample_names - gsub(\\.txt$, , basename(files)) count_list - lapply(files, function(f) { df - read.table(f, header TRUE, row.names 1, sep \t) return(df$count) })這里有個(gè)經(jīng)驗(yàn)如果你用GDC的下載方式文件名帶長(zhǎng)串UUID建議先重命名成樣本ID避免后面矩陣列名和分組表對(duì)不上?？梢杂靡粋€(gè)簡(jiǎn)單的批量重命名腳本或者直接在R里用sample_names映射總之要保持文件名和樣本ID的對(duì)應(yīng)關(guān)系清晰。5.2 構(gòu)建表達(dá)矩陣和分組信息expr_raw - do.call(cbind, count_list) colnames(expr_raw) - sample_names # 去掉Ensembl ID的小數(shù)版本號(hào)并注釋成symbol ensembl - gsub(\\..*, , rownames(expr_raw)) rownames(expr_raw) - ensembl # 用biomaRt做注釋也可以預(yù)存本地映射 mart - useMart(ensembl, dataset hsapiens_gene_ensembl) annot - getBM(attributes c(ensembl_gene_id, hgnc_symbol), filters ensembl_gene_id, values ensembl, mart mart) expr_symbol - expr_raw[annot$ensembl_gene_id, ] rownames(expr_symbol) - annot$hgnc_symbol # 處理重復(fù)symbol按行求和保留 expr_symbol - expr_symbol[!is.na(rownames(expr_symbol)), ] expr_symbol - expr_symbol[rownames(expr_symbol) ! , ] expr_symbol - as.data.frame(expr_symbol) %% rownames_to_column(symbol) %% group_by(symbol) %% summarise(across(everything(), sum)) %% column_to_rownames(symbol)需要注意如果biomaRt連接不穩(wěn)定建議一次性把所有Ensembl ID查完后把注釋結(jié)果存成csv后面重跑時(shí)直接read.csv讀取避免反復(fù)等待網(wǎng)絡(luò)。如果你的分析不要求轉(zhuǎn)symbol也可以保留Ensembl ID后續(xù)用注釋文件做功能富集時(shí)再映射各有各的方便。5.3 過(guò)濾、校正、標(biāo)準(zhǔn)化三步走# 假設(shè)sample_info包含sample_id, group, batch三列 # 確保矩陣列的順序與sample_info的sample_id完全一致 expr_raw - expr_raw[, sample_info$sample_id] # 1. 低表達(dá)基因過(guò)濾edgeR推薦方式 dge - DGEList(counts expr_raw, group sample_info$group) keep - filterByExpr(dge, group sample_info$group) dge - dge[keep, , keep.lib.sizes FALSE] # 2. 批次效應(yīng)校正ComBat_seq輸入raw count counts_corrected - ComBat_seq(counts dge$counts, batch sample_info$batch, group sample_info$group) # 3. 標(biāo)準(zhǔn)化轉(zhuǎn)換供可視化和常規(guī)差異篩選用 # 例如轉(zhuǎn)CPM后log2 expr_cpm - cpm(counts_corrected, log TRUE, prior.count 1) # 如果要喂給DESeq2做差異分析則用counts_corrected構(gòu)建DESeqDataSet dds - DESeqDataSetFromMatrix(countData counts_corrected, colData sample_info, design ~ group) dds - DESeq(dds) res - results(dds, contrast c(group, Tumor, Normal))這段代碼里ComBat_seq的group參數(shù)是必填的它在校正批次效應(yīng)的同時(shí)會(huì)盡量保留真實(shí)的組間差異。如果你漏了這個(gè)參數(shù)ComBat_seq會(huì)在無(wú)監(jiān)督模式下運(yùn)行可能會(huì)把真實(shí)的生物學(xué)差異也一并“校正”掉結(jié)果就是差異分析什么都篩不出來(lái)。5.4 驗(yàn)證處理效果處理完之后一定要驗(yàn)證別急著進(jìn)差異分析。常用兩個(gè)檢查第一重新跑一次PCA看Tumor和Normal是否按預(yù)期的分組分開(kāi)了第二繪制處理前后的批次聚類(lèi)熱圖對(duì)比確認(rèn)批次效應(yīng)有所緩解。如果PCA上樣本仍然明顯按批次聚類(lèi)說(shuō)明批次信息可能沒(méi)找對(duì)或者批次效應(yīng)與生物因素高度混雜可能需要更復(fù)雜的模型處理。PCA的可視化可以用基礎(chǔ)R也可以ggplot2畫(huà)比如提取前兩個(gè)主成分按樣本分組著色再用geom_text標(biāo)上樣本ID方便找離群點(diǎn)。這一步花不了五分鐘但能幫你避免跑到差異分析階段才發(fā)現(xiàn)數(shù)據(jù)質(zhì)量有問(wèn)題的大返工。6. 常見(jiàn)問(wèn)題與排查技巧實(shí)錄處理TCGA數(shù)據(jù)的過(guò)程里很多問(wèn)題都是反復(fù)出現(xiàn)的我把一些典型場(chǎng)景整理成速查表遇到問(wèn)題可以直接對(duì)照排查。6.1 常見(jiàn)報(bào)錯(cuò)場(chǎng)景速查表場(chǎng)景典型現(xiàn)象排查方向基因ID轉(zhuǎn)換后全是NAbiomaRt返回大量NA檢查Ensembl版本是否匹配考慮改用org.Hs.eg.db矩陣列名和分組表順序不一致DESeq2報(bào)錯(cuò)樣本不匹配用match()按順序重排列徹底解決順序問(wèn)題大量基因在過(guò)濾后仍然全零過(guò)濾條件太寬松或注釋比例低檢查注釋文件是否只覆蓋了蛋白編碼基因批次校正后組間差異反而變小ComBat_seq參數(shù)不當(dāng)檢查group參數(shù)是否正確指定不能用無(wú)監(jiān)督模式TPM矩陣跑DESeq2結(jié)果高度顯著但基因數(shù)異常多DESeq2不接收TPM改用TPM矩陣做線性差異篩選一個(gè)病人有多個(gè)樣本同一病人在Tumor和Normal組各出現(xiàn)多次按病人ID去重避免偽重復(fù)混淆檢驗(yàn)這些場(chǎng)景我在幫別人看代碼時(shí)幾乎都遇到過(guò)。尤其是矩陣列名順序的問(wèn)題看似小事跑DESeq2時(shí)一旦報(bào)錯(cuò)新手往往摸不著頭腦其實(shí)根源就是列順序不一致。用match函數(shù)把表達(dá)矩陣的列按sample_info的順序重排一下問(wèn)題立刻消失。6.2 避坑經(jīng)驗(yàn)我從這些錯(cuò)誤中學(xué)到的事最容易踩的坑是盲目照搬代碼。網(wǎng)上的教程常常直接用Xena下載好的矩陣但你要跑的是自己從GDC下載的per-sample文件流程就不一樣。我建議每一步都檢查一下中間產(chǎn)物的行數(shù)和列數(shù)至少確認(rèn)表達(dá)矩陣的基因數(shù)在過(guò)濾前后分別有多少樣本數(shù)是否和分組表完全一致。數(shù)據(jù)規(guī)模對(duì)不上后面跑出什么結(jié)果都不要覺(jué)得奇怪。第二個(gè)經(jīng)驗(yàn)是版本記錄。TCGA數(shù)據(jù)本身有版本更新比如GDC上同一個(gè)TCGA項(xiàng)目的表達(dá)數(shù)據(jù)會(huì)隨參考基因組版本更新而重新比對(duì)你下載時(shí)的release版本會(huì)直接影響Ensembl ID的注釋結(jié)果。建議把下載日期、數(shù)據(jù)版本、參考基因組信息都記在一個(gè)README文件里這不僅是可重復(fù)性的要求后面寫(xiě)論文方法部分也會(huì)需要。第三個(gè)經(jīng)驗(yàn)是時(shí)間成本管理。從GDC批量下載幾百個(gè)per-sample文件再合并如果網(wǎng)速不行會(huì)比較痛苦。我曾經(jīng)處理一個(gè)LUSC項(xiàng)目下載和整理就花了大半天后來(lái)改用TCGAbiolinks的GDCquery函數(shù)或直接從UCSC Xena拿TPM矩陣半小時(shí)內(nèi)搞定。根據(jù)自己的分析目標(biāo)選擇合適的數(shù)據(jù)獲取方式省下的時(shí)間足夠你多排查好幾個(gè)報(bào)錯(cuò)了。我做TCGA數(shù)據(jù)處理這幾年最大的體會(huì)就是前處理沒(méi)有想象中那么“機(jī)械”每一步都需要結(jié)合數(shù)據(jù)本身和分析目標(biāo)來(lái)做決定。同樣是差異分析前處理用DESeq2的人和用limma的人在過(guò)濾、標(biāo)準(zhǔn)化、批次校正的選擇上可能完全不一樣但核心邏輯是一致的讓數(shù)據(jù)干凈、可比、可解釋。上面這套流程是我自己反復(fù)用過(guò)的不敢說(shuō)最優(yōu)但至少能幫你少走幾段彎路。如果你在處理過(guò)程中遇到別的坑歡迎交流畢竟生信這條路很多經(jīng)驗(yàn)都是踩坑踩出來(lái)的。 TCGA差異分析前處理全流程：從數(shù)據(jù)格式到批次效應(yīng)校正
做TCGA數(shù)據(jù)挖掘很多人的精力都花在差異分析本身等跑完一看結(jié)果要么差異基因數(shù)少得可憐要么多到?jīng)]法解釋回頭排查才發(fā)現(xiàn)是差異分析前處理這步出了問(wèn)題。這篇文章把TCGA數(shù)據(jù)處理的完整前處理流程拆開(kāi)講清楚從數(shù)據(jù)格式辨析、下載渠道選擇、樣本分組、表達(dá)量矩陣構(gòu)建到低表達(dá)基因過(guò)濾、批次效應(yīng)校正、標(biāo)準(zhǔn)化轉(zhuǎn)換每個(gè)環(huán)節(jié)都說(shuō)清做什么、為什么這么做、以及我踩過(guò)的坑。適合正在做TCGA數(shù)據(jù)挖掘、尤其是不太清楚數(shù)據(jù)下載之后該怎么下手的生信新手參考。1. 先搞懂TCGA數(shù)據(jù)的三種形態(tài)不然后面全是坑TCGA全稱(chēng)是The Cancer Genome Atlas癌癥基因組圖譜這是由一個(gè)國(guó)家級(jí)大規(guī)模測(cè)序項(xiàng)目積累下來(lái)的多組學(xué)數(shù)據(jù)庫(kù)覆蓋33種癌癥類(lèi)型、超過(guò)2萬(wàn)例樣本包含了轉(zhuǎn)錄組、拷貝數(shù)變異、甲基化、突變等多層面數(shù)據(jù)。做癌癥相關(guān)生信分析的人基本繞不開(kāi)這個(gè)庫(kù)。但很多人下載完轉(zhuǎn)錄組表達(dá)數(shù)據(jù)后一打開(kāi)文件就懵了因?yàn)橥瑯右粋€(gè)TCGA項(xiàng)目你能下載到好幾種不同格式的表達(dá)量數(shù)據(jù)它們之間的數(shù)值含義完全不同用錯(cuò)了地方會(huì)讓后續(xù)所有分析都失去意義。1.1 Counts、FPKM、TPM到底有什么區(qū)別TCGA的轉(zhuǎn)錄組表達(dá)數(shù)據(jù)在GDC官方門(mén)戶(hù)和第三方平臺(tái)上最常見(jiàn)的就是三種格式HTSeq-Counts、FPKM/FPKM-UQ、TPM。簡(jiǎn)單說(shuō)HTSeq-Counts是基因比對(duì)后統(tǒng)計(jì)到的原始read計(jì)數(shù)它沒(méi)有做任何文庫(kù)大小和基因長(zhǎng)度的校正數(shù)值大小直接受測(cè)序深度影響同一批樣本之間如果測(cè)序量差異大raw counts就不適合直接拿來(lái)比較。FPKMFragments Per Kilobase Million是RNA-seq早期常用的標(biāo)準(zhǔn)化指標(biāo)計(jì)算時(shí)先按基因長(zhǎng)度歸一化再按文庫(kù)大小歸一化。它和TPM看起來(lái)很像但計(jì)算順序不同導(dǎo)致二者在數(shù)值分布上有細(xì)微差別。FPKM在跨樣本比較時(shí)其實(shí)存在一點(diǎn)系統(tǒng)偏差而TPM在實(shí)現(xiàn)上做了修正單位是“每百萬(wàn)條轉(zhuǎn)錄本中來(lái)自該基因的轉(zhuǎn)錄本數(shù)”樣本間可比性更好所以現(xiàn)在很多平臺(tái)默認(rèn)給TPM。那差異分析前處理到底該選哪個(gè)這取決于你后面用什么工具。如果你打算用DESeq2或edgeR跑標(biāo)準(zhǔn)差異分析官方推薦的輸入是raw count矩陣HTSeq-Counts因?yàn)檫@兩個(gè)工具內(nèi)部有自己的一套歸一化邏輯你給它的應(yīng)該是“沒(méi)被過(guò)度處理過(guò)的原始計(jì)數(shù)”。如果你只是想做基于表達(dá)量高低的分組比較、或者直接對(duì)表達(dá)值做t檢驗(yàn)相關(guān)的差異篩選那用TPM然后log2轉(zhuǎn)換是比較常見(jiàn)的做法。需要特別提醒絕對(duì)不要把FPKM或者TPM塞給DESeq2結(jié)果會(huì)非常難看差異基因數(shù)量動(dòng)不動(dòng)就成千上萬(wàn)明顯失真。1.2 數(shù)據(jù)下載GDC官方和UCSC Xena怎么選下載TCGA表達(dá)數(shù)據(jù)主要有兩條主流路徑。第一條是GDC Data Portalportal.gdc.cancer.gov這是官方數(shù)據(jù)倉(cāng)庫(kù)能拿到最原始的HTSeq-Counts文件下載時(shí)可以選擇per-sample的單獨(dú)文件也可以通過(guò)GDC API或者像TCGAbiolinks這樣的R包批量下載適合對(duì)數(shù)據(jù)溯源要求嚴(yán)格的課題。官方渠道的數(shù)據(jù)更新及時(shí)參考基因組版本和注釋信息都標(biāo)得很清楚寫(xiě)方法部分時(shí)引起來(lái)也方便。第二條是UCSC Xenaxenabrowser.net這個(gè)平臺(tái)把TCGA的數(shù)據(jù)做了統(tǒng)一整理可以按需求下載整理好的表達(dá)矩陣提供TPM和FPKM兩種格式還附帶臨床信息。用Xena最舒服的一點(diǎn)是它直接給你一個(gè)基因x樣本的矩陣不需要自己逐樣本合并對(duì)生信基礎(chǔ)不太強(qiáng)的人友好很多。它的數(shù)據(jù)也經(jīng)過(guò)了統(tǒng)一流程處理樣本名、基因名格式都比較規(guī)范。兩套方案怎么選我的建議是如果后續(xù)要用DESeq2/edgeR就走GDC拿raw count如果只做簡(jiǎn)單的表達(dá)量差異篩選Xena的TPM矩陣完全夠用。另外還有一個(gè)常見(jiàn)渠道是cBioPortal但它的表達(dá)數(shù)據(jù)多是處理過(guò)的適合查詢(xún)和可視化不太適合拿來(lái)當(dāng)差異分析的原始輸入。實(shí)際處理TCGA數(shù)據(jù)時(shí)經(jīng)常是下載完才發(fā)現(xiàn)格式不對(duì)所以第一步先把數(shù)據(jù)形態(tài)和后續(xù)分析工具對(duì)應(yīng)好能省下一大堆返工時(shí)間。2. 樣本分組前處理的第一個(gè)分水嶺拿到表達(dá)數(shù)據(jù)后的第一件事不是急著做差異分析而是把樣本分組搞清楚。TCGA的樣本編號(hào)不是隨便起的里面藏著樣本類(lèi)型、組織來(lái)源、是否配對(duì)等信息如果分組出錯(cuò)后面所有差異比較都是白做。2.1 腫瘤和正常樣本到底怎么區(qū)分TCGA樣本編號(hào)是有一套規(guī)則的四段式結(jié)構(gòu)比如TCGA-XX-XXXX-01A-11R-XXXX-XX其中第四段的開(kāi)頭兩位數(shù)字代表樣本類(lèi)型。01開(kāi)頭的是原發(fā)性實(shí)體瘤11開(kāi)頭的是癌旁正常組織solid tissue normal05是原發(fā)性血液腫瘤06是轉(zhuǎn)移性腫瘤還有02復(fù)發(fā)性腫瘤、03原發(fā)性血液腫瘤等少見(jiàn)類(lèi)型。日常分析中最常見(jiàn)的就是01和11兩類(lèi)。這里面最容易被忽略的一點(diǎn)是不是所有癌癥類(lèi)型都同時(shí)有01和11兩類(lèi)樣本。像GBM膠質(zhì)母細(xì)胞瘤、LGG腦低級(jí)別膠質(zhì)瘤這類(lèi)中樞神經(jīng)系統(tǒng)腫瘤以及一些罕見(jiàn)癌種正常對(duì)照樣本數(shù)量極少甚至沒(méi)有。這時(shí)候如果你還想著用TCGA自身樣本做配對(duì)差異分析樣本量就會(huì)非常尷尬。常見(jiàn)的替代方案是去GTEx數(shù)據(jù)庫(kù)拉正常組織數(shù)據(jù)來(lái)補(bǔ)充對(duì)照不過(guò)這又會(huì)引入跨數(shù)據(jù)集的批次效應(yīng)問(wèn)題處理起來(lái)要另外花心思。在寫(xiě)分組文件通常叫sample_info或者phenotype表的時(shí)候有幾個(gè)細(xì)節(jié)需要注意第一列名要簡(jiǎn)潔統(tǒng)一一般就是sample_id和group兩列第二組別命名建議用明確的class標(biāo)簽比如Tumor和Normal而不是01和11這種數(shù)字代號(hào)后面可視化的時(shí)候圖例會(huì)更清晰第三一定記得檢查有沒(méi)有重復(fù)的樣本ID有些樣本是多組織部位取樣如果你下載的時(shí)候沒(méi)嚴(yán)格過(guò)濾會(huì)出現(xiàn)一個(gè)病人對(duì)應(yīng)多條記錄的情況。2.2 表達(dá)量矩陣構(gòu)建和基因ID轉(zhuǎn)換如果你的數(shù)據(jù)是從GDC下載的per-sample文件每個(gè)樣本是一個(gè)獨(dú)立的counts文件那第一步就是把所有文件合并成一個(gè)矩陣。合并的邏輯很簡(jiǎn)單讀入所有文件把里面每行基因的counts值一一對(duì)應(yīng)到矩陣的列上。但實(shí)際操作有幾個(gè)很容易出問(wèn)題的坑。第一個(gè)坑是gene_id的格式。GDC下載的HTSeq-Counts文件里行名通常是Ensembl基因ID的帶版本格式比如ENSG00000242268.11。你要用sub函數(shù)把小數(shù)點(diǎn)和后面的版本號(hào)去掉再去做ID注釋轉(zhuǎn)換。對(duì)應(yīng)的基因symbol轉(zhuǎn)換最常用的是用biomaRt包實(shí)時(shí)查詢(xún)或者用org.Hs.eg.db包做本地注釋。我自己的經(jīng)驗(yàn)是網(wǎng)絡(luò)不穩(wěn)定時(shí)biomaRt很容易超時(shí)預(yù)先把Ensembl ID和symbol的映射關(guān)系保存成一份本地文件離線環(huán)境下也能直接轉(zhuǎn)換效率高很多。第二個(gè)坑是重復(fù)基因名的處理。Ensembl ID轉(zhuǎn)換到symbol之后大概率會(huì)出現(xiàn)多個(gè)ID對(duì)應(yīng)同一個(gè)symbol的情況也就是基因名重復(fù)。這時(shí)候不能直接保留需要按累積表達(dá)量或者最大表達(dá)量合并去重否則后面差異分析會(huì)報(bào)錯(cuò)或者產(chǎn)生重復(fù)行。還有一個(gè)細(xì)節(jié)是染色體上的小RNA、假基因等非蛋白編碼基因要不要過(guò)濾如果是常規(guī)mRNA表達(dá)譜差異分析我一般會(huì)保留protein_coding基因用注釋文件篩掉非編碼轉(zhuǎn)錄本能在源頭減少不少噪音。3. 低表達(dá)基因過(guò)濾這一步做不做結(jié)果差非常多處理完ID轉(zhuǎn)換和矩陣構(gòu)建之后接下來(lái)就是過(guò)濾。別小看這一步它對(duì)后續(xù)差異分析的影響非常大但也是最容易被新手跳過(guò)的一步。3.1 為什么不能把全基因集直接拿去做差異分析TCGA的HTSeq-Counts矩陣動(dòng)輒六萬(wàn)行但里面真正在樣本中穩(wěn)定表達(dá)的基因其實(shí)遠(yuǎn)沒(méi)有那么多。很多基因在絕大多數(shù)樣本里表達(dá)量為0或者只有幾個(gè)read這些基因就是噪音。你如果不做過(guò)濾直接跑DESeq2會(huì)產(chǎn)生兩個(gè)問(wèn)題一是多重檢驗(yàn)校正時(shí)需要比較的基因數(shù)量大幅膨脹padj會(huì)變嚴(yán)格一些真陽(yáng)性可能被壓掉二是大量全零或近零的基因會(huì)干擾后續(xù)離散度估計(jì)影響差異檢驗(yàn)的穩(wěn)健性。過(guò)濾的基本思路是設(shè)定一個(gè)保留標(biāo)準(zhǔn)比如“至少在20%的樣本中counts大于等于10”這類(lèi)規(guī)則。實(shí)際操作里DESeq2官方文檔建議的是一個(gè)簡(jiǎn)單但有據(jù)可依的過(guò)濾方式先對(duì)所有基因計(jì)算該基因在各樣本中的平均表達(dá)量或最大表達(dá)量再設(shè)一個(gè)閾值。edgeR則提供了filterByExpr函數(shù)可以綜合考慮最小計(jì)數(shù)、樣本量、分組信息來(lái)自動(dòng)推薦過(guò)濾條件。如果你用的是TPM矩陣做后續(xù)分析通常會(huì)把閾值設(shè)成TPM 1且在多少比例樣本中滿(mǎn)足條件。3.2 過(guò)濾閾值怎么定三個(gè)場(chǎng)景對(duì)比我整理了幾種常用的過(guò)濾標(biāo)準(zhǔn)你可以根據(jù)數(shù)據(jù)量和平時(shí)的習(xí)慣來(lái)選擇。過(guò)濾方式規(guī)則示例適用場(chǎng)景注意事項(xiàng)總量閾值型行和 10 或 mean counts 5全轉(zhuǎn)錄組初步篩選速度快偏寬松可能留下較多低表達(dá)基因比例閾值型至少在90%樣本中 counts 10關(guān)注優(yōu)勢(shì)表達(dá)基因的課題嚴(yán)格容易誤刪條件特異性表達(dá)基因工具推薦型filterByExpr自動(dòng)判斷常規(guī)差異分析準(zhǔn)備綜合庫(kù)大小和分組信息個(gè)人最常用第三種是edgeR推薦型底層會(huì)綜合所有樣本的測(cè)序深度和庫(kù)大小來(lái)設(shè)置閾值我個(gè)人用得最多因?yàn)樗诒Ａ粽鎸?shí)信號(hào)和去掉噪音之間平衡得比較好。順便說(shuō)一句這里有一個(gè)容易混淆的概念過(guò)濾和標(biāo)準(zhǔn)化誰(shuí)先誰(shuí)后的問(wèn)題。通常先過(guò)濾低表達(dá)基因再做標(biāo)準(zhǔn)化/歸一化邏輯上更順。原因是低表達(dá)基因的存在會(huì)影響某些標(biāo)準(zhǔn)化方法對(duì)文庫(kù)大小的估計(jì)先把確定是噪音的行去掉標(biāo)準(zhǔn)化會(huì)更穩(wěn)定。如果你用的是DESeq2它的median-of-ratios因素估計(jì)本身對(duì)低表達(dá)基因也敏感所以提前過(guò)濾是防患于未然。4. 批次效應(yīng)處理和數(shù)據(jù)標(biāo)準(zhǔn)化批次效應(yīng)是TCGA數(shù)據(jù)處理里最讓人頭疼的問(wèn)題之一但也是差異分析前必須面對(duì)的一關(guān)。TCGA的樣本不是一天之內(nèi)測(cè)完的樣本來(lái)源遍布多個(gè)組織中心、測(cè)序平臺(tái)和批次這些技術(shù)差異如果混進(jìn)你的分析里得到的結(jié)果很可能不是生物學(xué)差異而是技術(shù)噪音。4.1 批次效應(yīng)怎么發(fā)現(xiàn)主成分分析和聚類(lèi)熱圖批次效應(yīng)是指樣本在測(cè)序批次、文庫(kù)制備、芯片或平臺(tái)不同等因素影響下產(chǎn)生的系統(tǒng)性差異它和真實(shí)的生物學(xué)差異混雜在一起輕則讓PCA圖上的Tumor和Normal分不開(kāi)重則直接讓差異分析結(jié)果不可信。很多資料把批次效應(yīng)放在差分析之后才檢查但我建議在處理階段就提前看一遍免得后續(xù)返工。最直觀的方法是PCA。對(duì)log2標(biāo)準(zhǔn)化后的表達(dá)矩陣做PCA然后按樣本的分組信息、測(cè)序平臺(tái)如Illumina GA vs HiSeq、或者來(lái)源組織中心給點(diǎn)著色觀察樣本是否明顯按非生物學(xué)因素聚類(lèi)。如果發(fā)現(xiàn)按批次聚類(lèi)的現(xiàn)象很明顯就該考慮校正了。還有一個(gè)輔助方法是畫(huà)樣本相關(guān)性的熱圖如果同批次的樣本聚成了清晰的模塊而模塊內(nèi)既有Tumor也有Normal那基本可以判斷存在批次效應(yīng)。4.2 ComBat-seq和limma的removeBatchEffect怎么選處理批次效應(yīng)有不少工具最常用的是sva包的ComBat系列。ComBat適合處理微陣列和RNA-seq的表達(dá)矩陣近似連續(xù)數(shù)值ComBat-seq則是針對(duì)RNA-seq raw count專(zhuān)門(mén)開(kāi)發(fā)的它不改變數(shù)據(jù)的整數(shù)特性輸出結(jié)果可以繼續(xù)喂給DESeq2或edgeR做差異分析。如果你已經(jīng)把數(shù)據(jù)log2轉(zhuǎn)換成了連續(xù)值也可以用limma包的removeBatchEffect但它適合在標(biāo)準(zhǔn)化之后、差異分析之前對(duì)表達(dá)矩陣做殘差化處理。我的組合拳實(shí)踐經(jīng)驗(yàn)是先用filterByExpr過(guò)濾掉低表達(dá)基因再用ComBat_seq對(duì)raw counts校正批次校正完再跑DESeq2比較穩(wěn)健。如果你拿的是Xena的TPM矩陣TPM是連續(xù)值就log2(TPM1)之后用removeBatchEffect然后基于殘差矩陣做后續(xù)分析。這里有個(gè)細(xì)節(jié)批次信息最好是樣本的真實(shí)測(cè)序批次plate、seq center、tissue source site不要簡(jiǎn)單用“下載日期”替代因?yàn)橄螺d日期跟生物學(xué)變量完全無(wú)關(guān)反而可能引入新的混淆。4.3 log2轉(zhuǎn)換和標(biāo)準(zhǔn)化方法選擇的經(jīng)驗(yàn)log2轉(zhuǎn)換幾乎是TCGA表達(dá)數(shù)據(jù)可視化和差異篩選的標(biāo)配操作但要注意兩個(gè)容易被忽略的點(diǎn)。第一log2(x1)和log2(CPM1)的區(qū)別前者針對(duì)raw count或TPM后者是先將counts轉(zhuǎn)成CPM再做log2二者數(shù)值分布不同下游算法對(duì)輸入類(lèi)型敏感建議從頭到尾保持一致。第二log2轉(zhuǎn)換只適合方差穩(wěn)定的數(shù)據(jù)場(chǎng)景如果你要跑的是方差依賴(lài)的統(tǒng)計(jì)模型比如DESeq2負(fù)二項(xiàng)模型千萬(wàn)不能自己先log2再喂進(jìn)去應(yīng)該把raw counts原樣交給DESeq2處理。至于TPM和CPM的選擇我個(gè)人觀點(diǎn)是TCGA的TPM矩陣已經(jīng)考慮到基因長(zhǎng)度影響適合做表達(dá)定量但如果你的分析目標(biāo)是比較同一樣本內(nèi)部基因間的表達(dá)水平TPM比CPM合適。實(shí)際處理中從一個(gè)矩陣出發(fā)先明確下游分析類(lèi)型再?zèng)Q定采用哪種數(shù)據(jù)形態(tài)不要事到臨頭才來(lái)回切換切換過(guò)程中數(shù)值分布的變化很可能把你的差異分析結(jié)果帶偏。5. 實(shí)操?gòu)脑嘉募娇芍苯幼霾町惙治龅臄?shù)據(jù)理論說(shuō)了一大堆我放一套自己常用的R腳本流程出來(lái)。這個(gè)流程從GDC下載的per-sample HTSeq-Counts文件入手最終得到可以直接做差異分析的數(shù)據(jù)結(jié)構(gòu)。你只需準(zhǔn)備好兩個(gè)目錄一個(gè)放所有樣本的counts文件每個(gè)文件兩列g(shù)ene_id和count一個(gè)放樣本注釋表格包含sample_id、group和batch等列。5.1 準(zhǔn)備工作與讀取數(shù)據(jù)library(DESeq2) library(edgeR) library(sva) library(biomaRt) library(dplyr) library(tibble) # 讀入所有樣本的counts文件 files - list.files(path counts_dir, pattern *.txt, full.names TRUE) sample_names - gsub(\\.txt$, , basename(files)) count_list - lapply(files, function(f) { df - read.table(f, header TRUE, row.names 1, sep \t) return(df$count) })這里有個(gè)經(jīng)驗(yàn)如果你用GDC的下載方式文件名帶長(zhǎng)串UUID建議先重命名成樣本ID避免后面矩陣列名和分組表對(duì)不上?？梢杂靡粋€(gè)簡(jiǎn)單的批量重命名腳本或者直接在R里用sample_names映射總之要保持文件名和樣本ID的對(duì)應(yīng)關(guān)系清晰。5.2 構(gòu)建表達(dá)矩陣和分組信息expr_raw - do.call(cbind, count_list) colnames(expr_raw) - sample_names # 去掉Ensembl ID的小數(shù)版本號(hào)并注釋成symbol ensembl - gsub(\\..*, , rownames(expr_raw)) rownames(expr_raw) - ensembl # 用biomaRt做注釋也可以預(yù)存本地映射 mart - useMart(ensembl, dataset hsapiens_gene_ensembl) annot - getBM(attributes c(ensembl_gene_id, hgnc_symbol), filters ensembl_gene_id, values ensembl, mart mart) expr_symbol - expr_raw[annot$ensembl_gene_id, ] rownames(expr_symbol) - annot$hgnc_symbol # 處理重復(fù)symbol按行求和保留 expr_symbol - expr_symbol[!is.na(rownames(expr_symbol)), ] expr_symbol - expr_symbol[rownames(expr_symbol) ! , ] expr_symbol - as.data.frame(expr_symbol) %% rownames_to_column(symbol) %% group_by(symbol) %% summarise(across(everything(), sum)) %% column_to_rownames(symbol)需要注意如果biomaRt連接不穩(wěn)定建議一次性把所有Ensembl ID查完后把注釋結(jié)果存成csv后面重跑時(shí)直接read.csv讀取避免反復(fù)等待網(wǎng)絡(luò)。如果你的分析不要求轉(zhuǎn)symbol也可以保留Ensembl ID后續(xù)用注釋文件做功能富集時(shí)再映射各有各的方便。5.3 過(guò)濾、校正、標(biāo)準(zhǔn)化三步走# 假設(shè)sample_info包含sample_id, group, batch三列 # 確保矩陣列的順序與sample_info的sample_id完全一致 expr_raw - expr_raw[, sample_info$sample_id] # 1. 低表達(dá)基因過(guò)濾edgeR推薦方式 dge - DGEList(counts expr_raw, group sample_info$group) keep - filterByExpr(dge, group sample_info$group) dge - dge[keep, , keep.lib.sizes FALSE] # 2. 批次效應(yīng)校正ComBat_seq輸入raw count counts_corrected - ComBat_seq(counts dge$counts, batch sample_info$batch, group sample_info$group) # 3. 標(biāo)準(zhǔn)化轉(zhuǎn)換供可視化和常規(guī)差異篩選用 # 例如轉(zhuǎn)CPM后log2 expr_cpm - cpm(counts_corrected, log TRUE, prior.count 1) # 如果要喂給DESeq2做差異分析則用counts_corrected構(gòu)建DESeqDataSet dds - DESeqDataSetFromMatrix(countData counts_corrected, colData sample_info, design ~ group) dds - DESeq(dds) res - results(dds, contrast c(group, Tumor, Normal))這段代碼里ComBat_seq的group參數(shù)是必填的它在校正批次效應(yīng)的同時(shí)會(huì)盡量保留真實(shí)的組間差異。如果你漏了這個(gè)參數(shù)ComBat_seq會(huì)在無(wú)監(jiān)督模式下運(yùn)行可能會(huì)把真實(shí)的生物學(xué)差異也一并“校正”掉結(jié)果就是差異分析什么都篩不出來(lái)。5.4 驗(yàn)證處理效果處理完之后一定要驗(yàn)證別急著進(jìn)差異分析。常用兩個(gè)檢查第一重新跑一次PCA看Tumor和Normal是否按預(yù)期的分組分開(kāi)了第二繪制處理前后的批次聚類(lèi)熱圖對(duì)比確認(rèn)批次效應(yīng)有所緩解。如果PCA上樣本仍然明顯按批次聚類(lèi)說(shuō)明批次信息可能沒(méi)找對(duì)或者批次效應(yīng)與生物因素高度混雜可能需要更復(fù)雜的模型處理。PCA的可視化可以用基礎(chǔ)R也可以ggplot2畫(huà)比如提取前兩個(gè)主成分按樣本分組著色再用geom_text標(biāo)上樣本ID方便找離群點(diǎn)。這一步花不了五分鐘但能幫你避免跑到差異分析階段才發(fā)現(xiàn)數(shù)據(jù)質(zhì)量有問(wèn)題的大返工。6. 常見(jiàn)問(wèn)題與排查技巧實(shí)錄處理TCGA數(shù)據(jù)的過(guò)程里很多問(wèn)題都是反復(fù)出現(xiàn)的我把一些典型場(chǎng)景整理成速查表遇到問(wèn)題可以直接對(duì)照排查。6.1 常見(jiàn)報(bào)錯(cuò)場(chǎng)景速查表場(chǎng)景典型現(xiàn)象排查方向基因ID轉(zhuǎn)換后全是NAbiomaRt返回大量NA檢查Ensembl版本是否匹配考慮改用org.Hs.eg.db矩陣列名和分組表順序不一致DESeq2報(bào)錯(cuò)樣本不匹配用match()按順序重排列徹底解決順序問(wèn)題大量基因在過(guò)濾后仍然全零過(guò)濾條件太寬松或注釋比例低檢查注釋文件是否只覆蓋了蛋白編碼基因批次校正后組間差異反而變小ComBat_seq參數(shù)不當(dāng)檢查group參數(shù)是否正確指定不能用無(wú)監(jiān)督模式TPM矩陣跑DESeq2結(jié)果高度顯著但基因數(shù)異常多DESeq2不接收TPM改用TPM矩陣做線性差異篩選一個(gè)病人有多個(gè)樣本同一病人在Tumor和Normal組各出現(xiàn)多次按病人ID去重避免偽重復(fù)混淆檢驗(yàn)這些場(chǎng)景我在幫別人看代碼時(shí)幾乎都遇到過(guò)。尤其是矩陣列名順序的問(wèn)題看似小事跑DESeq2時(shí)一旦報(bào)錯(cuò)新手往往摸不著頭腦其實(shí)根源就是列順序不一致。用match函數(shù)把表達(dá)矩陣的列按sample_info的順序重排一下問(wèn)題立刻消失。6.2 避坑經(jīng)驗(yàn)我從這些錯(cuò)誤中學(xué)到的事最容易踩的坑是盲目照搬代碼。網(wǎng)上的教程常常直接用Xena下載好的矩陣但你要跑的是自己從GDC下載的per-sample文件流程就不一樣。我建議每一步都檢查一下中間產(chǎn)物的行數(shù)和列數(shù)至少確認(rèn)表達(dá)矩陣的基因數(shù)在過(guò)濾前后分別有多少樣本數(shù)是否和分組表完全一致。數(shù)據(jù)規(guī)模對(duì)不上后面跑出什么結(jié)果都不要覺(jué)得奇怪。第二個(gè)經(jīng)驗(yàn)是版本記錄。TCGA數(shù)據(jù)本身有版本更新比如GDC上同一個(gè)TCGA項(xiàng)目的表達(dá)數(shù)據(jù)會(huì)隨參考基因組版本更新而重新比對(duì)你下載時(shí)的release版本會(huì)直接影響Ensembl ID的注釋結(jié)果。建議把下載日期、數(shù)據(jù)版本、參考基因組信息都記在一個(gè)README文件里這不僅是可重復(fù)性的要求后面寫(xiě)論文方法部分也會(huì)需要。第三個(gè)經(jīng)驗(yàn)是時(shí)間成本管理。從GDC批量下載幾百個(gè)per-sample文件再合并如果網(wǎng)速不行會(huì)比較痛苦。我曾經(jīng)處理一個(gè)LUSC項(xiàng)目下載和整理就花了大半天后來(lái)改用TCGAbiolinks的GDCquery函數(shù)或直接從UCSC Xena拿TPM矩陣半小時(shí)內(nèi)搞定。根據(jù)自己的分析目標(biāo)選擇合適的數(shù)據(jù)獲取方式省下的時(shí)間足夠你多排查好幾個(gè)報(bào)錯(cuò)了。我做TCGA數(shù)據(jù)處理這幾年最大的體會(huì)就是前處理沒(méi)有想象中那么“機(jī)械”每一步都需要結(jié)合數(shù)據(jù)本身和分析目標(biāo)來(lái)做決定。同樣是差異分析前處理用DESeq2的人和用limma的人在過(guò)濾、標(biāo)準(zhǔn)化、批次校正的選擇上可能完全不一樣但核心邏輯是一致的讓數(shù)據(jù)干凈、可比、可解釋。上面這套流程是我自己反復(fù)用過(guò)的不敢說(shuō)最優(yōu)但至少能幫你少走幾段彎路。如果你在處理過(guò)程中遇到別的坑歡迎交流畢竟生信這條路很多經(jīng)驗(yàn)都是踩坑踩出來(lái)的。