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