Skip to content

统计读数

一般来说是统计比对到某个contig，某个基因，某个区域之类的的读数。然后换算为RPKM、FPKM、TPM等值，抑或是直接使用counts数来定量，再进行后面的差异分析。

其中，RPKM是Reads per kilo bases per million mapped reads，计算公式如下：

\[ RPKM = \frac{Reads\_per\_transcript}{\frac{total\_reads}{10^6} × \frac{transcript\_length}{10^3}} \]

其中10^3是为了标准化基因的长度，而10^6是为了标准化测序深度。FPKM是Fragments per kilo bases per million mapped reads，则是针对双端测序，与RPKM类似，计算方法相同。

但是，目前更多的会使用TPM或counts。其中，TPM是transcript per million，与RPKM/FPKM不同之处在于TPM先去除了基因的长度影响，而RPKM/FPKM则先去除测序深度的影响，实际上TPM改进了RPKM/FPKM等跨样本时定量的不准确。

\[ T = \sum^{N}_{i=1}{N_i / L_i} \]
\[ TPM = \frac{N_i / L_i × 10^6}{T} \]

其中\(N_i\)为比对到第i个外显子的reads数，\(L_i\)是第i个外显子的长度。

那么什么时候用counts，什么时候用TPM呢。比较容易得到的结论是在比较同一个基因在不同情况下表达量的差异时，基因长度对分析结果不会有影响，所有可以用counts。而在比较不同基因的表达差异时，用TPM可能会更好。

使用featureCounts或HTseq-count获得counts数同时加入了gtf文件用来注释，后续可以从结果矩阵里计算TPM，这样不管我们是counts还是用TPM来分析都可以。

featureCounts

featureCounts是subread中的一个工具，使用featureCounts进行counts的统计，需要用到gtf基因注释文件。命令中把所有需要的bam文件传入，下面是示例。最好把所有需要的bam文件扔到一个文件夹下。 

1
2
3
4
5
6
featureCounts \
    -T 10 \
    -g gene_id \
    -a annotation/mm39.ncbiRefSeq.gtf \
    -o results/final_featureCounts.txt \
    bam/*.bam

如果是双端的数据，还需要加入-p参数。

结果中的Successfully assigned alignments基本都在60%左右，这个值与gtf文件有关，正常应该在70%~90%左右。先不管，进入下一步差异分析。

HTseq

用HTseq-count来统计hisat2比对的结果。推荐使用上面的featureCounts，因为速度是真滴快太多。 

1
2
3
htseq-count -f bam -r pos -s no \
    -i gene_id bam/*.bam annotation/mm39.ncbiRefSeq.gtf \
    1>results/final_htseq.txt 2>results/final_htseq.summary