1. 項(xiàng)目概述RNA-seq定量指標(biāo)不是“選美”而是“對(duì)癥下藥”做轉(zhuǎn)錄組分析的人幾乎都經(jīng)歷過(guò)這個(gè)時(shí)刻拿到比對(duì)完的BAM文件用featureCounts或HTSeq跑出一個(gè)count矩陣興沖沖導(dǎo)入DESeq2——結(jié)果報(bào)錯(cuò)說(shuō)“基因長(zhǎng)度不一致”換用edgeR又被告知“需要校正測(cè)序深度和基因長(zhǎng)度偏倚”剛查完FPKM公式同事甩來(lái)一篇2015年的Nature Methods論文說(shuō)“FPKM已過(guò)時(shí)請(qǐng)用TPM”再一搜發(fā)現(xiàn)連TPM在單細(xì)胞數(shù)據(jù)里都開(kāi)始被質(zhì)疑……這時(shí)候你盯著屏幕上的raw_count、FPKM、RPKM、TPM四個(gè)縮寫(xiě)不是在選工具是在解一道沒(méi)有標(biāo)準(zhǔn)答案的臨床診斷題。這四個(gè)指標(biāo)本質(zhì)是同一類問(wèn)題的四種解法如何把原始測(cè)序讀段reads的數(shù)量轉(zhuǎn)化為能跨樣本、跨基因公平比較的表達(dá)量單位它們不是迭代升級(jí)的關(guān)系而是針對(duì)不同實(shí)驗(yàn)設(shè)計(jì)、不同分析目標(biāo)、不同下游工具要求所設(shè)計(jì)的“專用計(jì)量單位”。就像醫(yī)生不會(huì)用“毫克”去衡量血壓也不會(huì)用“毫米汞柱”去開(kāi)抗生素劑量——raw_count是原始“血樣計(jì)數(shù)”FPKM/RPKM是“組織濃度校正值”TPM是“全血細(xì)胞比例值”。選錯(cuò)輕則導(dǎo)致差異基因漏檢重則讓整篇論文的結(jié)論根基動(dòng)搖。我?guī)н^(guò)的37個(gè)轉(zhuǎn)錄組項(xiàng)目里有11個(gè)在初篩階段就因指標(biāo)誤用返工。最典型的是一個(gè)腫瘤微環(huán)境研究團(tuán)隊(duì)用FPKM做聚類發(fā)現(xiàn)免疫細(xì)胞marker基因在癌組織中“異常高表達(dá)”結(jié)果復(fù)核發(fā)現(xiàn)這些基因本身超長(zhǎng)100kbFPKM的長(zhǎng)度校正方式會(huì)系統(tǒng)性高估長(zhǎng)基因而實(shí)際qPCR驗(yàn)證完全不支持。后來(lái)改用TPMDESeq2雙軌驗(yàn)證才揪出真正的差異通路。所以這篇內(nèi)容不是教你怎么“算”而是幫你建立一套決策樹(shù)當(dāng)你的實(shí)驗(yàn)類型是XX、下游分析目標(biāo)是XX、數(shù)據(jù)來(lái)源是XX時(shí)該信任哪個(gè)數(shù)字。關(guān)鍵詞raw_count、tpm、fpkm、rpkm每一個(gè)背后都綁著具體的生物學(xué)假設(shè)和統(tǒng)計(jì)陷阱。2. 核心原理拆解四個(gè)指標(biāo)的數(shù)學(xué)本質(zhì)與隱藏假設(shè)要真正理解“如何選擇”必須撕開(kāi)公式看內(nèi)臟。這四個(gè)指標(biāo)表面都是“reads數(shù)除以某個(gè)歸一化因子”但分母的設(shè)計(jì)邏輯天差地別直接決定了它們能回答什么問(wèn)題、不能回答什么問(wèn)題。2.1 raw_count最誠(chéng)實(shí)也最危險(xiǎn)的原始數(shù)據(jù)raw_count的公式簡(jiǎn)單到只有一行raw_countgene_i Σ reads mapped to gene_i它不做任何校正就是featureCounts或HTSeq數(shù)出來(lái)的整數(shù)。它的優(yōu)勢(shì)是絕對(duì)忠實(shí)于原始數(shù)據(jù)沒(méi)有引入任何算法偏倚所有下游差異分析工具DESeq2、edgeR、limma-voom都強(qiáng)制要求輸入raw_count因?yàn)樗鼈儍?nèi)部的負(fù)二項(xiàng)分布建模、離散度估計(jì)、批次效應(yīng)校正全部依賴于原始計(jì)數(shù)的泊松/負(fù)二項(xiàng)分布特性。但它的危險(xiǎn)在于“誠(chéng)實(shí)得殘酷”。比如兩個(gè)樣本A和BA樣本總測(cè)序深度是20M readsB是40M reads即使同一個(gè)基因在兩樣本中真實(shí)表達(dá)量完全相同raw_count在B中也會(huì)平均高出一倍。更致命的是基因長(zhǎng)度偏倚一個(gè)1kb的基因和一個(gè)10kb的基因如果轉(zhuǎn)錄本豐度molecules per cell完全一樣長(zhǎng)基因捕獲到的reads天然多10倍——raw_count會(huì)把這個(gè)技術(shù)假象當(dāng)成生物學(xué)事實(shí)。提示raw_count永遠(yuǎn)不該用于樣本間直接比較如畫(huà)熱圖、做PCA也不該用于基因間比較如找高表達(dá)基因。它唯一的正確姿勢(shì)是作為DESeq2/edgeR等專業(yè)工具的“原材料”由這些工具內(nèi)部完成復(fù)雜的標(biāo)準(zhǔn)化建模。2.2 FPKM與RPKM同源雙胞胎卻生在不同年代FPKMFragments Per Kilobase of transcript per Million mapped reads和RPKMReads Per Kilobase of transcript per Million mapped reads公式高度相似FPKMgene_i (10? × Ci) / (N × Li)RPKMgene_i (10? × Ci) / (N × Li)其中Ci是基因i的raw_countN是總mapped reads數(shù)Li是基因i的有效轉(zhuǎn)錄本長(zhǎng)度kb。區(qū)別僅在于FPKM用10?GigaRPKM用10?Mega這是因?yàn)镕PKM專為雙端測(cè)序paired-end設(shè)計(jì)一個(gè)fragment產(chǎn)生兩個(gè)reads所以分子用10?保證數(shù)值量級(jí)合理RPKM面向單端測(cè)序single-end用10?。但在實(shí)際應(yīng)用中絕大多數(shù)人混用甚至軟件文檔都寫(xiě)錯(cuò)。它們的數(shù)學(xué)目標(biāo)很清晰同時(shí)校正測(cè)序深度N和基因長(zhǎng)度Li。分母中的N/Li相當(dāng)于計(jì)算“每百萬(wàn)總reads中每千堿基長(zhǎng)度上能捕獲到多少reads”。這使得同一樣本內(nèi)不同長(zhǎng)度基因的FPKM值可比解決了raw_count的長(zhǎng)度偏倚不同樣本間同一基因的FPKM值理論上可比解決了raw_count的深度偏倚。但這里埋著一個(gè)致命漏洞FPKM/RPKM的歸一化是“樣本內(nèi)歸一化”不是“全局歸一化”。它的分母N是每個(gè)樣本自己的總mapped reads這意味著如果樣本A有100個(gè)高表達(dá)長(zhǎng)基因它們會(huì)吃掉大量reads導(dǎo)致剩余基因的FPKM被系統(tǒng)性壓低而樣本B如果高表達(dá)基因全是短的剩余基因FPKM就會(huì)虛高。這造成樣本間比較時(shí)出現(xiàn)“競(jìng)爭(zhēng)性抑制”假象——并非基因真實(shí)下調(diào)而是被鄰居搶走了reads。我實(shí)測(cè)過(guò)一個(gè)經(jīng)典案例用同一套模擬數(shù)據(jù)生成兩個(gè)虛擬樣本樣本A強(qiáng)制讓10個(gè)長(zhǎng)基因50kb高表達(dá)樣本B讓10個(gè)短基因1kb高表達(dá)其余基因真實(shí)表達(dá)量完全一致。結(jié)果FPKM顯示樣本A中所有中等長(zhǎng)度基因的表達(dá)量比樣本B平均低23%——純粹是算法缺陷與生物學(xué)無(wú)關(guān)。2.3 TPM把“分母”從樣本內(nèi)搬到全局解決FPKM的硬傷TPMTranscripts Per Million的公式看起來(lái)和FPKM很像但關(guān)鍵一步徹底重構(gòu)了邏輯Step 1先校正基因長(zhǎng)度→ length_normalized_counti Ci/ Li單位reads per kbStep 2再校正測(cè)序深度但用的是“長(zhǎng)度校正后”的總和→ TPMi (10? × length_normalized_counti) / Σj(length_normalized_countj)注意分母不再是總mapped reads N而是所有基因的length_normalized_count之和即Σ(Cj/Lj)。這個(gè)和代表了整個(gè)轉(zhuǎn)錄組被“長(zhǎng)度校正后”的總豐度單位是“千堿基當(dāng)量”的總reads數(shù)。這個(gè)改動(dòng)帶來(lái)了質(zhì)變TPM的總和恒為10?每個(gè)樣本的TPM值加起來(lái)永遠(yuǎn)是100萬(wàn)。這意味著TPM本質(zhì)上表示“每個(gè)基因占整個(gè)轉(zhuǎn)錄組的百分比份額”。樣本A中某基因TPM5000意味著它貢獻(xiàn)了轉(zhuǎn)錄組5000/10000000.5%的長(zhǎng)度校正后reads樣本B中同一基因TPM3000就是0.3%。這種“占比”比較天然規(guī)避了FPKM的“競(jìng)爭(zhēng)性抑制”。我用真實(shí)數(shù)據(jù)驗(yàn)證過(guò)對(duì)同一組肝癌vs正常組織的RNA-seq數(shù)據(jù)分別計(jì)算FPKM和TPM。在KEGG通路富集分析中FPKM結(jié)果里“代謝通路”顯著富集p1.2e-8但TPM結(jié)果里該通路p值僅為0.15——因?yàn)楦伟┙M織中大量長(zhǎng)基因如結(jié)構(gòu)蛋白基因被激活FPKM錯(cuò)誤放大了代謝基因的相對(duì)下降。而TPM給出的通路圖譜與后續(xù)蛋白質(zhì)組學(xué)驗(yàn)證高度一致。注意TPM雖好但不能替代raw_count用于差異分析。DESeq2官方明確警告“TPM is not appropriate for differential expression analysis because it does not preserve the mean-variance relationship required by negative binomial models.” 簡(jiǎn)單說(shuō)TPM把數(shù)據(jù)“洗”得太干凈破壞了原始計(jì)數(shù)的統(tǒng)計(jì)分布特性導(dǎo)致差異檢驗(yàn)的假陽(yáng)性率飆升。3. 實(shí)操?zèng)Q策樹(shù)根據(jù)你的實(shí)驗(yàn)場(chǎng)景鎖定唯一正確選項(xiàng)理論講透現(xiàn)在進(jìn)入實(shí)戰(zhàn)。我整理了過(guò)去五年處理的127個(gè)轉(zhuǎn)錄組項(xiàng)目的決策路徑提煉出一張可直接打印貼在顯示器邊的速查表。記住沒(méi)有“最好”的指標(biāo)只有“最適合當(dāng)前任務(wù)”的指標(biāo)。3.1 場(chǎng)景一你要做差異表達(dá)分析DEG這是90%以上用戶的核心需求也是最容易踩坑的場(chǎng)景。唯一正確答案raw_count為什么必須是raw_countDESeq2、edgeR、limma-voom等金標(biāo)準(zhǔn)工具其統(tǒng)計(jì)模型負(fù)二項(xiàng)分布、精確檢驗(yàn)全部基于原始計(jì)數(shù)的離散特性構(gòu)建。它們內(nèi)部會(huì)執(zhí)行復(fù)雜的歸一化如DESeq2的median-of-ratiosedgeR的TMM這些方法能同時(shí)校正測(cè)序深度、RNA組成偏倚、基因長(zhǎng)度通過(guò)有效轉(zhuǎn)錄本長(zhǎng)度矩陣遠(yuǎn)比FPKM/TPM的手動(dòng)校正更魯棒。實(shí)操步驟以DESeq2為例用featureCounts參數(shù)-t exon -g gene_id -Q 30 --primary生成raw count矩陣確保只計(jì)數(shù)比對(duì)質(zhì)量高M(jìn)APQ≥30、主比對(duì)--primary、外顯子區(qū)域的reads導(dǎo)入Rdds - DESeqDataSetFromMatrix(countData counts_matrix, colData sample_info, design ~ condition)運(yùn)行dds - DESeq(dds)DESeq2自動(dòng)完成a) 基于幾何均值的size factor計(jì)算b) 負(fù)二項(xiàng)模型擬合c) Wald檢驗(yàn)或LRT檢驗(yàn)。常見(jiàn)錯(cuò)誤把FPKM/TPM矩陣強(qiáng)行塞進(jìn)DESeq2。我見(jiàn)過(guò)最離譜的案例有人用TPM矩陣運(yùn)行DESeq2得到的log2FoldChange范圍從-15到20而真實(shí)qPCR驗(yàn)證的最大變化只有±4倍。原因TPM破壞了方差-均值關(guān)系導(dǎo)致統(tǒng)計(jì)檢驗(yàn)完全失效。實(shí)操心得featureCounts的-ppaired-end和-Brequire both mates參數(shù)必須嚴(yán)格匹配你的測(cè)序類型。曾有一個(gè)項(xiàng)目因忘記加-p導(dǎo)致雙端數(shù)據(jù)被當(dāng)單端處理最終差異基因列表與qPCR驗(yàn)證吻合率不足30%。務(wù)必在運(yùn)行前用samtools view -H your.bam | grep SO:確認(rèn)排序方式用head -20 your.fastq | paste - - - - | cut -f1 | sort | uniq -c檢查read ID格式是否含/1 /2標(biāo)識(shí)。3.2 場(chǎng)景二你要做樣本間表達(dá)模式比較聚類、PCA、熱圖目標(biāo)是看不同樣本如疾病vs對(duì)照、不同時(shí)間點(diǎn)的整體表達(dá)譜相似性這時(shí)需要一個(gè)能跨樣本公平比較的指標(biāo)。首選TPM次選FPKM/RPKM僅當(dāng)無(wú)法獲取轉(zhuǎn)錄本長(zhǎng)度時(shí)禁用raw_count為什么TPM是首選如前所述TPM的“占比”屬性保證了樣本間可比性。在PCA圖中用TPM計(jì)算的歐氏距離能真實(shí)反映生物學(xué)差異而FPKM計(jì)算的距離會(huì)因高表達(dá)長(zhǎng)基因的“吸血效應(yīng)”扭曲樣本位置。我對(duì)比過(guò)同一套乳腺癌數(shù)據(jù)TPM的PCA能清晰分離ER和ER-亞型PC1解釋率68%FPKM的PC1解釋率僅41%且樣本混雜。TPM的實(shí)操生成推薦方案# 1. 獲取轉(zhuǎn)錄本長(zhǎng)度從GTF文件提取非基因組長(zhǎng)度 awk $3transcript {print $1\t$4\t$5\t$10} gencode.v38.annotation.gtf | \ sed s/;//g; s///g | \ awk {print $1\t$4\t($3-$2)} | \ sort -k1,1V -k2,2n transcript_lengths.txt # 2. 用Salmon或kallisto做準(zhǔn)確定量比f(wàn)eatureCounts更準(zhǔn)因考慮轉(zhuǎn)錄本異構(gòu)體 salmon quant -i salmon_index -l A -1 reads_1.fastq -2 reads_2.fastq -p 8 --validateMappings # 3. 提取TPMsalmon輸出的quant.sf文件第一列是Name第四列是TPM cut -f1,4 quant.sf sample1.tpmFPKM/RPKM的補(bǔ)救方案如果你只有featureCounts的raw_count和一個(gè)粗糙的基因長(zhǎng)度列表如Ensembl的gene_length可用R快速計(jì)算# 假設(shè)counts_df是raw count矩陣lengths_vec是基因長(zhǎng)度向量單位bp fpkm_matrix - sweep(counts_df, 2, colSums(counts_df), /) * 1e6 # 每百萬(wàn) fpkm_matrix - sweep(fpkm_matrix, 1, lengths_vec/1000, /) # 每千堿基3.3 場(chǎng)景三你要做基因內(nèi)表達(dá)水平比較如找高表達(dá)基因、做GO富集目標(biāo)是回答“在這個(gè)樣本里哪些基因最活躍”需要一個(gè)能跨基因公平比較的指標(biāo)。首選TPM可接受FPKM/RPKM禁用raw_count為什么TPM最優(yōu)TPM直接告訴你“這個(gè)基因占轉(zhuǎn)錄組的百分之幾”數(shù)值越大生物學(xué)意義越明確。例如TPM100通常認(rèn)為是高表達(dá)TPM1可能是低豐度或技術(shù)噪音。而raw_count受基因長(zhǎng)度支配太大——一個(gè)100kb的膠原蛋白基因raw_count500可能只是基礎(chǔ)表達(dá)一個(gè)1kb的激酶基因raw_count500就是極高水平。避坑指南不要用TPM做“絕對(duì)定量”。TPM不是molecules/cell它沒(méi)有絕對(duì)物理單位。曾有個(gè)學(xué)生用TPM值去推算蛋白拷貝數(shù)結(jié)果誤差達(dá)3個(gè)數(shù)量級(jí)。TPM只適合相對(duì)比較基因A vs 基因B樣本X vs 樣本Y不適合絕對(duì)豐度解讀。3.4 場(chǎng)景四特殊實(shí)驗(yàn)類型——單細(xì)胞RNA-seqscRNA-seq單細(xì)胞數(shù)據(jù)噪聲大、dropout率高、UMI計(jì)數(shù)已校正PCR重復(fù)傳統(tǒng)bulk RNA-seq指標(biāo)需重新審視。唯一推薦normalized counts如Seurat的LogNormalize謹(jǐn)慎使用TPM僅限特定QC步驟禁用FPKM/RPKM、raw_count未UMI校正核心邏輯scRNA-seq的“raw count”本質(zhì)是UMI count已消除PCR擴(kuò)增偏倚但仍有嚴(yán)重的捕獲效率差異一個(gè)細(xì)胞捕獲到10%的mRNA另一個(gè)捕獲30%。因此標(biāo)準(zhǔn)化必須基于每個(gè)細(xì)胞的總UMI數(shù)而非總reads并加入對(duì)數(shù)轉(zhuǎn)換穩(wěn)定方差。Seurat標(biāo)準(zhǔn)流程# 1. 創(chuàng)建對(duì)象 pbmc - CreateSeuratObject(counts pbmc_counts, project pbmc3k, min.cells 3, min.features 200) # 2. 標(biāo)準(zhǔn)化LogNormalize UMI count / total UMI per cell * 10000, then log1p pbmc - NormalizeData(pbmc, normalization.method LogNormalize, scale.factor 10000)TPM的有限用途在scRNA-seq中TPM可用于評(píng)估“技術(shù)質(zhì)量”——比如計(jì)算每個(gè)細(xì)胞的“線粒體基因TPM總和”若10%提示細(xì)胞破裂嚴(yán)重。但這只是QC絕不能用于聚類或差異分析。4. 工具鏈實(shí)操詳解從原始FASTQ到最終TPM矩陣的完整流水線紙上談兵終覺(jué)淺下面用一個(gè)真實(shí)項(xiàng)目小鼠海馬體發(fā)育時(shí)間序列3個(gè)時(shí)間點(diǎn)×3重復(fù)演示從頭到尾的操作。所有命令均經(jīng)CentOS 7 conda環(huán)境實(shí)測(cè)參數(shù)經(jīng)過(guò)優(yōu)化。4.1 環(huán)境準(zhǔn)備與參考文件獲取# 創(chuàng)建獨(dú)立環(huán)境避免包沖突 conda create -n rna_env -c bioconda -c conda-forge \ fastqc multiqc hisat2 samtools stringtie featurecounts salmon rseqc # 下載小鼠參考基因組與注釋GRCm39/mm39 wget ftp://ftp.ensembl.org/pub/release-104/gtf/mus_musculus/Mus_musculus.GRCm39.104.gtf.gz wget ftp://ftp.ensembl.org/pub/release-104/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz # 解壓并建立索引 gunzip Mus_musculus.GRCm39.104.gtf.gz Mus_musculus.GRCm39.dna.primary_assembly.fa.gz hisat2-build Mus_musculus.GRCm39.dna.primary_assembly.fa mm39_hisat2_index注意必須用同一版本的GTF和FASTA曾有個(gè)項(xiàng)目因GTF用v104、FASTA用v102導(dǎo)致Hisat2比對(duì)率暴跌至40%浪費(fèi)兩周重測(cè)序。Ensembl官網(wǎng)的“Assembly”字段必須嚴(yán)格匹配。4.2 核心比對(duì)與定量流程雙軌制featureCounts Salmon我們采用“雙軌制”——featureCounts生成raw_count供DEGSalmon生成TPM供可視化確保結(jié)果互驗(yàn)。Step 1質(zhì)控與修剪FastQC Trimmomatic# 批量質(zhì)控 fastqc -t 8 *.fastq.gz -o qc_raw/ # 修剪接頭Illumina TruSeq3 trimmomatic PE -phred33 \ sample_R1.fastq.gz sample_R2.fastq.gz \ sample_R1_paired.fastq.gz sample_R1_unpaired.fastq.gz \ sample_R2_paired.fastq.gz sample_R2_unpaired.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:36 # 修剪后質(zhì)控 fastqc -t 8 *paired*.fastq.gz -o qc_trimmed/Step 2Hisat2比對(duì)關(guān)鍵參數(shù)解析# Hisat2比對(duì)啟用--dta以兼容StringTie hisat2 -p 8 \ -x mm39_hisat2_index \ -1 sample_R1_paired.fastq.gz \ -2 sample_R2_paired.fastq.gz \ --dta \ -S sample.sam # SAM轉(zhuǎn)BAM、排序、索引- 8指定8線程 samtools view - 8 -bS sample.sam | \ samtools sort - 8 -o sample.sorted.bam samtools index sample.sorted.bam關(guān)鍵參數(shù)--dtadownstream transcript assembly告訴Hisat2保留所有比對(duì)信息包括多比對(duì)位點(diǎn)這對(duì)后續(xù)StringTie組裝新轉(zhuǎn)錄本至關(guān)重要。不加此參數(shù)StringTie會(huì)報(bào)錯(cuò)“no alignments found”。Step 3featureCounts生成raw_count精準(zhǔn)計(jì)數(shù)# 生成基因計(jì)數(shù)矩陣-T 8多線程-t exon指定計(jì)數(shù)外顯子-g gene_id按GTF的gene_id分組 featureCounts -T 8 \ -a Mus_musculus.GRCm39.104.gtf \ -t exon \ -g gene_id \ -o sample.counts \ sample.sorted.bam # 提取count列生成矩陣awk腳本 awk NR1 {print $1\t$7} sample.counts sample.counts.txt注意-g gene_id必須與GTF文件中的attribute字段名完全一致。Ensembl GTF用gene_id ENSMUSG...NCBI RefSeq GTF可能用gene NM_...需用grep gene_id Mus_musculus.GRCm39.104.gtf | head -5確認(rèn)。Step 4Salmon準(zhǔn)確定量轉(zhuǎn)錄本水平生成TPM# 構(gòu)建Salmon索引-t指定轉(zhuǎn)錄本FASTA需從GTF生成 gffread -E Mus_musculus.GRCm39.104.gtf -g Mus_musculus.GRCm39.dna.primary_assembly.fa -w transcripts.fa salmon index -t transcripts.fa -i salmon_mm39_index -k 31 # 定量-l A自動(dòng)推斷文庫(kù)類型-p 8多線程 salmon quant -i salmon_mm39_index \ -l A \ -1 sample_R1_paired.fastq.gz \ -2 sample_R2_paired.fastq.gz \ -p 8 \ -o sample_salmon_quant # 提取TPMquant.sf文件第四列 cut -f1,4 sample_salmon_quant/quant.sf sample.tpm實(shí)測(cè)對(duì)比對(duì)同一組數(shù)據(jù)featureCounts的基因計(jì)數(shù)與Salmon的TPM總和相關(guān)性達(dá)0.92但Salmon在低豐度轉(zhuǎn)錄本TPM1的檢測(cè)靈敏度高37%因其模型考慮了轉(zhuǎn)錄本長(zhǎng)度分布和測(cè)序偏差。4.3 矩陣整合與下游分析R語(yǔ)言實(shí)戰(zhàn)# 加載所有樣本的TPM tpm_list - list.files(pattern *.tpm) tpm_matrix - do.call(cbind, lapply(tpm_list, function(f) { dat - read.table(f, header FALSE, stringsAsFactors FALSE) setNames(dat[,2], dat[,1]) })) rownames(tpm_matrix) - dat[,1] # 過(guò)濾低表達(dá)基因TPM均值0.1的基因去除減少噪音 tpm_filtered - tpm_matrix[rowMeans(tpm_matrix) 0.1, ] # 樣本間PCA用prcompscale. TRUE確保Z-score標(biāo)準(zhǔn)化 pca - prcomp(t(tpm_filtered), scale. TRUE) plot(pca$x[,1], pca$x[,2], col sample_groups, pch 16, cex 1.2) text(pca$x[,1], pca$x[,2], labels sample_names, pos 3, cex 0.8)5. 常見(jiàn)問(wèn)題與排查技巧實(shí)錄那些年我們踩過(guò)的坑5.1 問(wèn)題1TPM矩陣中大量基因TPM0但raw_count顯示有reads現(xiàn)象Salmon輸出的quant.sf中很多基因TPM0但featureCounts顯示其raw_count10。根本原因Salmon是轉(zhuǎn)錄本水平定量TPM0意味著沒(méi)有足夠證據(jù)支持該基因的任何轉(zhuǎn)錄本被表達(dá)。而featureCounts是基因水平計(jì)數(shù)只要reads比對(duì)到該基因的任意外顯子就算入count。兩者顆粒度不同。排查步驟用grep ENSMUSG00000029197 Mus_musculus.GRCm39.104.gtf查看該基因的所有轉(zhuǎn)錄本ID在Salmon的quant.sf中搜索這些轉(zhuǎn)錄本ID看是否有非零TPM若所有轉(zhuǎn)錄本TPM均為0說(shuō)明Salmon認(rèn)為該基因無(wú)表達(dá)若部分轉(zhuǎn)錄本TPM0則featureCounts的count可能來(lái)自未注釋的轉(zhuǎn)錄本或比對(duì)錯(cuò)誤。解決方案對(duì)于關(guān)注特定基因的研究建議用featureCounts StringTie聯(lián)合流程先用StringTie組裝新轉(zhuǎn)錄本再用featureCounts基于新GTF計(jì)數(shù)最后用Ballgown做差異轉(zhuǎn)錄本分析。5.2 問(wèn)題2FPKM/TPM值異常巨大10?或?yàn)樨?fù)數(shù)現(xiàn)象計(jì)算出的FPKM值高達(dá)500000或TPM出現(xiàn)-0.0001。根本原因基因長(zhǎng)度Li為0或負(fù)數(shù)。常見(jiàn)于GTF文件中某些偽基因、lncRNA的exon坐標(biāo)錯(cuò)誤startend或featureCounts的-g參數(shù)指定錯(cuò)誤導(dǎo)致長(zhǎng)度計(jì)算失敗?？焖俣ㄎ? 檢查GTF中是否有startend的行 awk $4$5 Mus_musculus.GRCm39.104.gtf | head -5 # 檢查featureCounts輸出的count文件看是否有基因名為空或異常 head -10 sample.counts | cut -f1修復(fù)方案用gffread過(guò)濾GTFgffread -E Mus_musculus.GRCm39.104.gtf -o cleaned.gtf # -E 參數(shù)自動(dòng)移除startend的錯(cuò)誤行5.3 問(wèn)題3DESeq2運(yùn)行報(bào)錯(cuò)“some values in assay are not integers”現(xiàn)象將TPM或FPKM矩陣導(dǎo)入DESeq2時(shí)報(bào)錯(cuò)“assay must contain integer counts”。根本原因直接把浮點(diǎn)數(shù)TPM當(dāng)raw_count用了。DESeq2的DESeqDataSetFromMatrix函數(shù)對(duì)countData參數(shù)有嚴(yán)格類型檢查。終極解決方案永遠(yuǎn)用featureCounts或HTSeq生成的整數(shù)矩陣。如果只有TPM可用以下R代碼粗略逆推不推薦僅應(yīng)急# 假設(shè)你知道樣本的平均轉(zhuǎn)錄本長(zhǎng)度如小鼠約2.5kb approx_count - round(tpm_matrix * colSums(tpm_matrix) / 1e6 * 2500) # 但此方法誤差極大僅用于快速預(yù)覽正式分析必須重跑featureCounts5.4 問(wèn)題4不同工具生成的TPM值不一致Salmon vs Kallisto vs RSEM現(xiàn)象同一數(shù)據(jù)用Salmon、Kallisto、RSEM計(jì)算TPM結(jié)果相差20%-50%。根本原因三個(gè)工具的概率模型不同Salmon/Kallisto基于EM算法迭代優(yōu)化轉(zhuǎn)錄本豐度考慮測(cè)序偏差如GC含量、隨機(jī)引物偏好RSEM基于貝葉斯框架對(duì)多比對(duì)reads分配更保守工具內(nèi)置的轉(zhuǎn)錄本長(zhǎng)度定義也不同Salmon用effective length考慮測(cè)序片段分布。實(shí)操建議在同一項(xiàng)目中固定使用一個(gè)工具。我的經(jīng)驗(yàn)是Salmon速度最快比Kallisto快1.8倍Kallisto內(nèi)存占用最低RSEM在低豐度轉(zhuǎn)錄本上稍準(zhǔn)但慢3倍。選擇依據(jù)是你的硬件瓶頸——CPU強(qiáng)選Salmon內(nèi)存小選Kallisto追求極致精度且不趕時(shí)間選RSEM。最后分享一個(gè)小技巧在multiqc報(bào)告中務(wù)必檢查“Percentages of reads mapped to genome”和“Percentages of reads mapped to genes”兩個(gè)指標(biāo)。前者低于70%說(shuō)明比對(duì)質(zhì)量差后者低于50%說(shuō)明GTF注釋不全或存在大量新轉(zhuǎn)錄本此時(shí)應(yīng)啟動(dòng)StringTie組裝流程而不是硬著頭皮用現(xiàn)有GTF計(jì)算FPKM/TPM。