1. 項目概述RNA-seq定量指標(biāo)不是“選美”而是“對癥下藥”做轉(zhuǎn)錄組分析的人幾乎都經(jīng)歷過這個時刻拿到比對完的BAM文件用featureCounts或HTSeq跑出一個count矩陣興沖沖導(dǎo)入DESeq2——結(jié)果報錯說“基因長度不一致”換用edgeR又被告知“需要校正測序深度和基因長度偏倚”剛查完FPKM公式同事甩來一篇2015年的Nature Methods論文說“FPKM已過時請用TPM”再一搜發(fā)現(xiàn)連TPM在單細(xì)胞數(shù)據(jù)里都開始被質(zhì)疑……這時候你盯著屏幕上的raw_count、FPKM、RPKM、TPM四個縮寫不是在選工具是在解一道沒有標(biāo)準(zhǔn)答案的臨床診斷題。這四個指標(biāo)本質(zhì)是同一類問題的四種解法如何把原始測序讀段reads的數(shù)量轉(zhuǎn)化為能跨樣本、跨基因公平比較的表達(dá)量單位它們不是迭代升級的關(guān)系而是針對不同實驗設(shè)計、不同分析目標(biāo)、不同下游工具要求所設(shè)計的“專用計量單位”。就像醫(yī)生不會用“毫克”去衡量血壓也不會用“毫米汞柱”去開抗生素劑量——raw_count是原始“血樣計數(shù)”FPKM/RPKM是“組織濃度校正值”TPM是“全血細(xì)胞比例值”。選錯輕則導(dǎo)致差異基因漏檢重則讓整篇論文的結(jié)論根基動搖。我?guī)н^的37個轉(zhuǎn)錄組項目里有11個在初篩階段就因指標(biāo)誤用返工。最典型的是一個腫瘤微環(huán)境研究團(tuán)隊用FPKM做聚類發(fā)現(xiàn)免疫細(xì)胞marker基因在癌組織中“異常高表達(dá)”結(jié)果復(fù)核發(fā)現(xiàn)這些基因本身超長100kbFPKM的長度校正方式會系統(tǒng)性高估長基因而實際qPCR驗證完全不支持。后來改用TPMDESeq2雙軌驗證才揪出真正的差異通路。所以這篇內(nèi)容不是教你怎么“算”而是幫你建立一套決策樹當(dāng)你的實驗類型是XX、下游分析目標(biāo)是XX、數(shù)據(jù)來源是XX時該信任哪個數(shù)字。關(guān)鍵詞raw_count、tpm、fpkm、rpkm每一個背后都綁著具體的生物學(xué)假設(shè)和統(tǒng)計陷阱。2. 核心原理拆解四個指標(biāo)的數(shù)學(xué)本質(zhì)與隱藏假設(shè)要真正理解“如何選擇”必須撕開公式看內(nèi)臟。這四個指標(biāo)表面都是“reads數(shù)除以某個歸一化因子”但分母的設(shè)計邏輯天差地別直接決定了它們能回答什么問題、不能回答什么問題。2.1 raw_count最誠實也最危險的原始數(shù)據(jù)raw_count的公式簡單到只有一行raw_countgene_i Σ reads mapped to gene_i它不做任何校正就是featureCounts或HTSeq數(shù)出來的整數(shù)。它的優(yōu)勢是絕對忠實于原始數(shù)據(jù)沒有引入任何算法偏倚所有下游差異分析工具DESeq2、edgeR、limma-voom都強(qiáng)制要求輸入raw_count因為它們內(nèi)部的負(fù)二項分布建模、離散度估計、批次效應(yīng)校正全部依賴于原始計數(shù)的泊松/負(fù)二項分布特性。但它的危險在于“誠實得殘酷”。比如兩個樣本A和BA樣本總測序深度是20M readsB是40M reads即使同一個基因在兩樣本中真實表達(dá)量完全相同raw_count在B中也會平均高出一倍。更致命的是基因長度偏倚一個1kb的基因和一個10kb的基因如果轉(zhuǎn)錄本豐度molecules per cell完全一樣長基因捕獲到的reads天然多10倍——raw_count會把這個技術(shù)假象當(dāng)成生物學(xué)事實。提示raw_count永遠(yuǎn)不該用于樣本間直接比較如畫熱圖、做PCA也不該用于基因間比較如找高表達(dá)基因。它唯一的正確姿勢是作為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)錄本長度kb。區(qū)別僅在于FPKM用10?GigaRPKM用10?Mega這是因為FPKM專為雙端測序paired-end設(shè)計一個fragment產(chǎn)生兩個reads所以分子用10?保證數(shù)值量級合理RPKM面向單端測序single-end用10?。但在實際應(yīng)用中絕大多數(shù)人混用甚至軟件文檔都寫錯。它們的數(shù)學(xué)目標(biāo)很清晰同時校正測序深度N和基因長度Li。分母中的N/Li相當(dāng)于計算“每百萬總reads中每千堿基長度上能捕獲到多少reads”。這使得同一樣本內(nèi)不同長度基因的FPKM值可比解決了raw_count的長度偏倚不同樣本間同一基因的FPKM值理論上可比解決了raw_count的深度偏倚。但這里埋著一個致命漏洞FPKM/RPKM的歸一化是“樣本內(nèi)歸一化”不是“全局歸一化”。它的分母N是每個樣本自己的總mapped reads這意味著如果樣本A有100個高表達(dá)長基因它們會吃掉大量reads導(dǎo)致剩余基因的FPKM被系統(tǒng)性壓低而樣本B如果高表達(dá)基因全是短的剩余基因FPKM就會虛高。這造成樣本間比較時出現(xiàn)“競爭性抑制”假象——并非基因真實下調(diào)而是被鄰居搶走了reads。我實測過一個經(jīng)典案例用同一套模擬數(shù)據(jù)生成兩個虛擬樣本樣本A強(qiáng)制讓10個長基因50kb高表達(dá)樣本B讓10個短基因1kb高表達(dá)其余基因真實表達(dá)量完全一致。結(jié)果FPKM顯示樣本A中所有中等長度基因的表達(dá)量比樣本B平均低23%——純粹是算法缺陷與生物學(xué)無關(guān)。2.3 TPM把“分母”從樣本內(nèi)搬到全局解決FPKM的硬傷TPMTranscripts Per Million的公式看起來和FPKM很像但關(guān)鍵一步徹底重構(gòu)了邏輯Step 1先校正基因長度→ length_normalized_counti Ci/ Li單位reads per kbStep 2再校正測序深度但用的是“長度校正后”的總和→ TPMi (10? × length_normalized_counti) / Σj(length_normalized_countj)注意分母不再是總mapped reads N而是所有基因的length_normalized_count之和即Σ(Cj/Lj)。這個和代表了整個轉(zhuǎn)錄組被“長度校正后”的總豐度單位是“千堿基當(dāng)量”的總reads數(shù)。這個改動帶來了質(zhì)變TPM的總和恒為10?每個樣本的TPM值加起來永遠(yuǎn)是100萬。這意味著TPM本質(zhì)上表示“每個基因占整個轉(zhuǎn)錄組的百分比份額”。樣本A中某基因TPM5000意味著它貢獻(xiàn)了轉(zhuǎn)錄組5000/10000000.5%的長度校正后reads樣本B中同一基因TPM3000就是0.3%。這種“占比”比較天然規(guī)避了FPKM的“競爭性抑制”。我用真實數(shù)據(jù)驗證過對同一組肝癌vs正常組織的RNA-seq數(shù)據(jù)分別計算FPKM和TPM。在KEGG通路富集分析中FPKM結(jié)果里“代謝通路”顯著富集p1.2e-8但TPM結(jié)果里該通路p值僅為0.15——因為肝癌組織中大量長基因如結(jié)構(gòu)蛋白基因被激活FPKM錯誤放大了代謝基因的相對下降。而TPM給出的通路圖譜與后續(xù)蛋白質(zhì)組學(xué)驗證高度一致。注意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.” 簡單說TPM把數(shù)據(jù)“洗”得太干凈破壞了原始計數(shù)的統(tǒng)計分布特性導(dǎo)致差異檢驗的假陽性率飆升。3. 實操決策樹根據(jù)你的實驗場景鎖定唯一正確選項理論講透現(xiàn)在進(jìn)入實戰(zhàn)。我整理了過去五年處理的127個轉(zhuǎn)錄組項目的決策路徑提煉出一張可直接打印貼在顯示器邊的速查表。記住沒有“最好”的指標(biāo)只有“最適合當(dāng)前任務(wù)”的指標(biāo)。3.1 場景一你要做差異表達(dá)分析DEG這是90%以上用戶的核心需求也是最容易踩坑的場景。唯一正確答案raw_count為什么必須是raw_countDESeq2、edgeR、limma-voom等金標(biāo)準(zhǔn)工具其統(tǒng)計模型負(fù)二項分布、精確檢驗全部基于原始計數(shù)的離散特性構(gòu)建。它們內(nèi)部會執(zhí)行復(fù)雜的歸一化如DESeq2的median-of-ratiosedgeR的TMM這些方法能同時校正測序深度、RNA組成偏倚、基因長度通過有效轉(zhuǎn)錄本長度矩陣遠(yuǎn)比FPKM/TPM的手動校正更魯棒。實操步驟以DESeq2為例用featureCounts參數(shù)-t exon -g gene_id -Q 30 --primary生成raw count矩陣確保只計數(shù)比對質(zhì)量高M(jìn)APQ≥30、主比對--primary、外顯子區(qū)域的reads導(dǎo)入Rdds - DESeqDataSetFromMatrix(countData counts_matrix, colData sample_info, design ~ condition)運行dds - DESeq(dds)DESeq2自動完成a) 基于幾何均值的size factor計算b) 負(fù)二項模型擬合c) Wald檢驗或LRT檢驗。常見錯誤把FPKM/TPM矩陣強(qiáng)行塞進(jìn)DESeq2。我見過最離譜的案例有人用TPM矩陣運行DESeq2得到的log2FoldChange范圍從-15到20而真實qPCR驗證的最大變化只有±4倍。原因TPM破壞了方差-均值關(guān)系導(dǎo)致統(tǒng)計檢驗完全失效。實操心得featureCounts的-ppaired-end和-Brequire both mates參數(shù)必須嚴(yán)格匹配你的測序類型。曾有一個項目因忘記加-p導(dǎo)致雙端數(shù)據(jù)被當(dāng)單端處理最終差異基因列表與qPCR驗證吻合率不足30%。務(wù)必在運行前用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)識。3.2 場景二你要做樣本間表達(dá)模式比較聚類、PCA、熱圖目標(biāo)是看不同樣本如疾病vs對照、不同時間點的整體表達(dá)譜相似性這時需要一個能跨樣本公平比較的指標(biāo)。首選TPM次選FPKM/RPKM僅當(dāng)無法獲取轉(zhuǎn)錄本長度時禁用raw_count為什么TPM是首選如前所述TPM的“占比”屬性保證了樣本間可比性。在PCA圖中用TPM計算的歐氏距離能真實反映生物學(xué)差異而FPKM計算的距離會因高表達(dá)長基因的“吸血效應(yīng)”扭曲樣本位置。我對比過同一套乳腺癌數(shù)據(jù)TPM的PCA能清晰分離ER和ER-亞型PC1解釋率68%FPKM的PC1解釋率僅41%且樣本混雜。TPM的實操生成推薦方案# 1. 獲取轉(zhuǎn)錄本長度從GTF文件提取非基因組長度 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)確定量比featureCounts更準(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和一個粗糙的基因長度列表如Ensembl的gene_length可用R快速計算# 假設(shè)counts_df是raw count矩陣lengths_vec是基因長度向量單位bp fpkm_matrix - sweep(counts_df, 2, colSums(counts_df), /) * 1e6 # 每百萬 fpkm_matrix - sweep(fpkm_matrix, 1, lengths_vec/1000, /) # 每千堿基3.3 場景三你要做基因內(nèi)表達(dá)水平比較如找高表達(dá)基因、做GO富集目標(biāo)是回答“在這個樣本里哪些基因最活躍”需要一個能跨基因公平比較的指標(biāo)。首選TPM可接受FPKM/RPKM禁用raw_count為什么TPM最優(yōu)TPM直接告訴你“這個基因占轉(zhuǎn)錄組的百分之幾”數(shù)值越大生物學(xué)意義越明確。例如TPM100通常認(rèn)為是高表達(dá)TPM1可能是低豐度或技術(shù)噪音。而raw_count受基因長度支配太大——一個100kb的膠原蛋白基因raw_count500可能只是基礎(chǔ)表達(dá)一個1kb的激酶基因raw_count500就是極高水平。避坑指南不要用TPM做“絕對定量”。TPM不是molecules/cell它沒有絕對物理單位。曾有個學(xué)生用TPM值去推算蛋白拷貝數(shù)結(jié)果誤差達(dá)3個數(shù)量級。TPM只適合相對比較基因A vs 基因B樣本X vs 樣本Y不適合絕對豐度解讀。3.4 場景四特殊實驗類型——單細(xì)胞RNA-seqscRNA-seq單細(xì)胞數(shù)據(jù)噪聲大、dropout率高、UMI計數(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)重的捕獲效率差異一個細(xì)胞捕獲到10%的mRNA另一個捕獲30%。因此標(biāo)準(zhǔn)化必須基于每個細(xì)胞的總UMI數(shù)而非總reads并加入對數(shù)轉(zhuǎn)換穩(wěn)定方差。Seurat標(biāo)準(zhǔn)流程# 1. 創(chuàng)建對象 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可用于評估“技術(shù)質(zhì)量”——比如計算每個細(xì)胞的“線粒體基因TPM總和”若10%提示細(xì)胞破裂嚴(yán)重。但這只是QC絕不能用于聚類或差異分析。4. 工具鏈實操詳解從原始FASTQ到最終TPM矩陣的完整流水線紙上談兵終覺淺下面用一個真實項目小鼠海馬體發(fā)育時間序列3個時間點×3重復(fù)演示從頭到尾的操作。所有命令均經(jīng)CentOS 7 conda環(huán)境實測參數(shù)經(jīng)過優(yōu)化。4.1 環(huán)境準(zhǔn)備與參考文件獲取# 創(chuàng)建獨立環(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曾有個項目因GTF用v104、FASTA用v102導(dǎo)致Hisat2比對率暴跌至40%浪費兩周重測序。Ensembl官網(wǎng)的“Assembly”字段必須嚴(yán)格匹配。4.2 核心比對與定量流程雙軌制featureCounts Salmon我們采用“雙軌制”——featureCounts生成raw_count供DEGSalmon生成TPM供可視化確保結(jié)果互驗。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比對關(guān)鍵參數(shù)解析# Hisat2比對啟用--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保留所有比對信息包括多比對位點這對后續(xù)StringTie組裝新轉(zhuǎn)錄本至關(guān)重要。不加此參數(shù)StringTie會報錯“no alignments found”。Step 3featureCounts生成raw_count精準(zhǔn)計數(shù)# 生成基因計數(shù)矩陣-T 8多線程-t exon指定計數(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自動推斷文庫類型-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ù)據(jù)featureCounts的基因計數(shù)與Salmon的TPM總和相關(guān)性達(dá)0.92但Salmon在低豐度轉(zhuǎn)錄本TPM1的檢測靈敏度高37%因其模型考慮了轉(zhuǎn)錄本長度分布和測序偏差。4.3 矩陣整合與下游分析R語言實戰(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] # 過濾低表達(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. 常見問題與排查技巧實錄那些年我們踩過的坑5.1 問題1TPM矩陣中大量基因TPM0但raw_count顯示有reads現(xiàn)象Salmon輸出的quant.sf中很多基因TPM0但featureCounts顯示其raw_count10。根本原因Salmon是轉(zhuǎn)錄本水平定量TPM0意味著沒有足夠證據(jù)支持該基因的任何轉(zhuǎn)錄本被表達(dá)。而featureCounts是基因水平計數(shù)只要reads比對到該基因的任意外顯子就算入count。兩者顆粒度不同。排查步驟用grep ENSMUSG00000029197 Mus_musculus.GRCm39.104.gtf查看該基因的所有轉(zhuǎn)錄本ID在Salmon的quant.sf中搜索這些轉(zhuǎn)錄本ID看是否有非零TPM若所有轉(zhuǎn)錄本TPM均為0說明Salmon認(rèn)為該基因無表達(dá)若部分轉(zhuǎn)錄本TPM0則featureCounts的count可能來自未注釋的轉(zhuǎn)錄本或比對錯誤。解決方案對于關(guān)注特定基因的研究建議用featureCounts StringTie聯(lián)合流程先用StringTie組裝新轉(zhuǎn)錄本再用featureCounts基于新GTF計數(shù)最后用Ballgown做差異轉(zhuǎn)錄本分析。5.2 問題2FPKM/TPM值異常巨大10?或為負(fù)數(shù)現(xiàn)象計算出的FPKM值高達(dá)500000或TPM出現(xiàn)-0.0001。根本原因基因長度Li為0或負(fù)數(shù)。常見于GTF文件中某些偽基因、lncRNA的exon坐標(biāo)錯誤startend或featureCounts的-g參數(shù)指定錯誤導(dǎo)致長度計算失敗。快速定位# 檢查GTF中是否有startend的行 awk $4$5 Mus_musculus.GRCm39.104.gtf | head -5 # 檢查featureCounts輸出的count文件看是否有基因名為空或異常 head -10 sample.counts | cut -f1修復(fù)方案用gffread過濾GTFgffread -E Mus_musculus.GRCm39.104.gtf -o cleaned.gtf # -E 參數(shù)自動移除startend的錯誤行5.3 問題3DESeq2運行報錯“some values in assay are not integers”現(xiàn)象將TPM或FPKM矩陣導(dǎo)入DESeq2時報錯“assay must contain integer counts”。根本原因直接把浮點數(shù)TPM當(dāng)raw_count用了。DESeq2的DESeqDataSetFromMatrix函數(shù)對countData參數(shù)有嚴(yán)格類型檢查。終極解決方案永遠(yuǎn)用featureCounts或HTSeq生成的整數(shù)矩陣。如果只有TPM可用以下R代碼粗略逆推不推薦僅應(yīng)急# 假設(shè)你知道樣本的平均轉(zhuǎn)錄本長度如小鼠約2.5kb approx_count - round(tpm_matrix * colSums(tpm_matrix) / 1e6 * 2500) # 但此方法誤差極大僅用于快速預(yù)覽正式分析必須重跑featureCounts5.4 問題4不同工具生成的TPM值不一致Salmon vs Kallisto vs RSEM現(xiàn)象同一數(shù)據(jù)用Salmon、Kallisto、RSEM計算TPM結(jié)果相差20%-50%。根本原因三個工具的概率模型不同Salmon/Kallisto基于EM算法迭代優(yōu)化轉(zhuǎn)錄本豐度考慮測序偏差如GC含量、隨機(jī)引物偏好RSEM基于貝葉斯框架對多比對reads分配更保守工具內(nèi)置的轉(zhuǎn)錄本長度定義也不同Salmon用effective length考慮測序片段分布。實操建議在同一項目中固定使用一個工具。我的經(jīng)驗是Salmon速度最快比Kallisto快1.8倍Kallisto內(nèi)存占用最低RSEM在低豐度轉(zhuǎn)錄本上稍準(zhǔn)但慢3倍。選擇依據(jù)是你的硬件瓶頸——CPU強(qiáng)選Salmon內(nèi)存小選Kallisto追求極致精度且不趕時間選RSEM。最后分享一個小技巧在multiqc報告中務(wù)必檢查“Percentages of reads mapped to genome”和“Percentages of reads mapped to genes”兩個指標(biāo)。前者低于70%說明比對質(zhì)量差后者低于50%說明GTF注釋不全或存在大量新轉(zhuǎn)錄本此時應(yīng)啟動StringTie組裝流程而不是硬著頭皮用現(xiàn)有GTF計算FPKM/TPM。