一行命令将count转为CPM/TPM/FPKM

一行命令将count转为CPM/TPM/FPKM 的软件为rnanorm,是一个基于Python开发的命令行工具。安装可以通过命令安装：

pip install rnanorm

我以featureCounts的输出文件进行举例，用featureCounts对进行基因count计数

featureCounts -T 10 -t exon -g gene_id -Q 10 --primary -p \
    -a gencode.vM25.annotation.gtf -F GTF -o sample.count \
    mapping/*.fil.bam 

得到的gene count在sample.count文件里

# tail -n +2 sample.count 是排除第一行
$ tail -n +2 sample.count | sed '1s/Geneid/FEATURE_ID/g' | cut -f 1,7- > sample.count.tsv
$ head sample.count.tsv
FEATURE_ID      mapping/RNA.fb.e11.5.rep1.fil.bam       mapping/RNA.fb.e11.5.rep2.fil.bam       mapping/RNA.fb.e16.5.rep1.fil.bam       mapping/RNA.fb.e16.5.rep2.fil.bam
ENSMUSG00000102693.1    1       0       0       0
ENSMUSG00000064842.1    1       0       0       0
ENSMUSG00000051951.5    6       11      16      16
ENSMUSG00000102851.1    2       0       0       0
ENSMUSG00000103377.1    1       2       3       0
ENSMUSG00000104017.1    1       1       0       0
ENSMUSG00000103025.1    0       3       2       2
ENSMUSG00000089699.1    0       0       0       0
ENSMUSG00000103201.1    0       1       2       1

用cut将基因ID一列，和count数值的列提取出来。用sed将Geneid换成FEATURE_ID，因为当前版本rnanorm（ 1.5.1）要求第一列的基因ID列名必须为FEATURE_ID

然后就是一行代码将count转为CPM/TPM/FPKM。

rnanorm sample.count.tsv --annotation gencode.vM25.annotation.gtf \
    --cpm-output sample.count.cpm.tsv
    --tpm-output sample.count.tpm.tsv \
    --fpkm-output sample.count.fpkm.tsv \

• 位置参数为基因count文件
• --annotation 基因组注释GTF文件
• --cpm-output CPM输出文件
• --tpm-output TPM输出文件
• --fpkm-output FPKM输出文件

后面的前面

RNA-Seq测序得到的fastq数据比对至参考基因组，得到SAM/BAM文件，然后进行read count计数，得到基因表达信息。然而由于不同基因长度不同，基因长，自然容易获得更read count；不同批次数据测序量不同，数据测序量多，自然read count就会高了。所以不同基因表达量不能直接进行比较，需要进行标准化。

要对基因长度进行矫正，就除以基因长度:

• 1000bp长度基因A count=1000, 矫正后1000/1000=1
• 5000bp长度的基因B count=5000, 矫正后5000/5000=1

要对测序深度进行矫正，就除以总测序量:

• 总read count=100000, 基因A read count 1000, 矫正后1000/100000
• 总read count=200000, 基因B read count 2000, 矫正后2000/200000

CPM

CPM (count per million),对测序深度做了一个矫正。

CPM = (gene_read_count / total_count ) * 10^6

至于为啥乘一个million（10^6），我猜大概是因为gene_read_count / total_count数值太小了，不方便人查看，又英文中计数单位有个，十，百，千，百万，十亿。十亿（billion）较大了， 千(thousand)又太小了， million这个数量级刚好差不多，就乘这么个数值，在百万层级看基因表达量。

RPKM

RPKM(reads per kilobase per million),对测序深度和长度做了一个矫正。用于单端测序数据。

RPKM = (gene_read_count / (total_count*gene_length) * 10^6 * 10^3

基因count除以总read count，除以gene长度。然后避免数值太小，根据测序深度大小乘10^6, 基因长度就乘以一个10^3。

FPKM

RPKM(fragments per kilobase per million) 用于双端数据，一个fragment是一对reads。故FPKM = RPKM / 2

TPM

TPM(transcript per million), 基因FPKM 占总的FPKM的比例， TPM衡量基因在样本中的相对表达量，对同一个基因不同样本，其FPKM可能相同，但相对所有基因而言，其在不同样本中可能所占比例不同。

TPM = gene_fpkm / total_fpkm * 10^6