reads计数原理及featureCounts统计counts后

摘录：生信技能树论坛 Kai statquest

要求：实现这个功能的软件也很多，还是烦请大家先自己搜索几个教程，入门请统一用htseq-count，对每个样本都会输出一个表达量文件。需要用脚本合并所有的样本为表达矩阵。参考：生信编程直播第四题：多个同样的行列式文件合并起来对这个表达矩阵可以自己简单在excel或者R里面摸索，求平均值，方差。
看看一些生物学意义特殊的基因表现如何，比如GAPDH,β-ACTIN等等

基因（转录本定量）

定量分为三个水平：基因水平（gene-level）转录本水平（transcript-level）外显子使用水平（exon-usage-level）

基因水平

在基因水平上，常用的软件为HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等，主要用来解决的问题是根据read和基因位置的overlap判断这个read到底是谁家的孩子，不同的工具对multimapping reads处理方式也是不同的，例如 HTSeq-count直接当他们不存在。而 Qualimap则是一人一份，平均分配。

image

对每个基因计数之后得到的count matrix 再后续的分析中，要注意标准化的问题。

比较同一样本（within-sample）不同基因之间的表达情况，需要考虑到转录本长度

比较不同样本（across sample）同一基因表达情况，不必在意转录本长度，但是要考虑到测序深度

除此之外，还必须要考虑到GC%所导致的偏差，以及测序仪器的系统偏差。read count标准化的算法有RPKM（SE）, FPKM（PE），TPM, TMM。
RPKM和FPKM, 前者是以每个reads作为一个单位，在单端测序中应用较多；而后者是以fragment作为一个单位，主要应用在双端测序后的分析。
有一些下游分析的软件会要求是输入的count matrix 是原始数据，未经过标准化，比如说DESeq2，这个时候需要注意上一步使用软件会不会进行标准化。
RPKM：Reads Per Kilobase per Million mapped reads的缩写，RPKM是将map到基因的read数除以map到基因组上的所有read数（以million为单位）与RNA的长度（以KB为单位）。

image.png

RNA-seq是二代测序技术中用来表示基因表达量或丰度的方法。在衡量基因表达量时，若是单纯以map到的read数来计算基因的表达量，在统计上是不合理的。因为在随机抽样的情况下，序列较长的基因被抽到的机率本来就会比序列短的基因较高，如此一来，序列长的基因永远会被认为表达量较高，而错估基因真正的表现量 。

FPKM：Fragments Per Kilobase of transcript per Million fragments mapped，每1百万个map上的reads中map到map到外显子的每1k个碱基上的fragments个数。如果一对paired reads都比对上了，那么这对reads当做一个fragments；如果paired reads中一个比对上，另外一个没比对上，那么比对上的那个reads当做一个fragments。

image.png
FPKM和RPKM的对比

1. FPKM和RPKM计算方法基本一致，主要一个采用reads一个采用Fragments，片段比读段含义更广。因此FPKM包含含义更广。
2. 由于FPKM与RPKM的唯一差别在于前者在reads map上的情况下只计数1，而后者会计数2；所以两者的公式其实是一样的。
3. RPKM不仅对测序深度做了归一化，也对基因长度做了归一化，使得不同长度的基因在不同测序深度下得到的基因表达水平估计具有了可比性，是目前最为常用的基因估计方法。

TPM: Transcripts per million 每百万读段中来自某基因的读段数，考虑了测序深度对读段计数的影响，但没有基因组长度的影响，主要用于miRNA测序。

Snipaste_2020-02-28_10-51-54.png
FPKM、RPKM、TPM的使用优劣
Snipaste_2020-02-28_10-53-06.png

详细应学习https://statquest.org/2015/07/09/rpkm-fpkm-and-tpm-clearly-explained/上的相关解释。

适用于featureCounts数counts后的cpm和tpm计算

参考：落寞的橙子

setwd("~/Desktop")#设置工作目录

countdata<-read.table("counts.txt",skip = 1,sep="\t",header = T,row.names = 1)

metadata <- countdata[,1:5]#提取基因信息count数据前的几列

countdata <- countdata[,6:ncol(countdata)]#提取counts数，counts数据主题部分

prefix<-"couts"#设置输出文件前缀名

cpm <- t(t(countdata)/colSums(countdata) * 1000000)#参考cpm定义

avg_cpm <- data.frame(avg_cpm=rowMeans(cpm))

#-----TPM Calculation------

kb <- metadata$Length / 1000

rpk <- countdata / kb

tpm <- t(t(rpk)/colSums(rpk) * 1000000)

avg_tpm <- data.frame(avg_tpm=rowMeans(tpm))

write.csv(avg_tpm,paste0(prefix,"_avg_tpm.csv"))

write.csv(avg_cpm,paste0(prefix,"_avg_cpm.csv"))

write.csv(tpm,paste0(prefix,"_tpm.csv"))

write.csv(cpm,paste0(prefix,"_cpm.csv"))

Jimmy老师的TPM版本

TPM公式，就是RPKM的百分比的百万倍扩大值

rpkm <- function(counts, lengths) {
rate <- counts / lengths
rate / sum(counts) * 1e9
}

exprSet_rpkm=rpkm(exprSet,len)
exprSet_tpm=1e6*exprSet_rpkm/colSums(exprSet_rpkm)

PS: 这个colSums(exprSet_tpm) 是接近一百万，而不是精确的一百万。