关于RNA-seq 的那点事Count 数的标准化 （一) RPKM 和FPKM，TPM及C(R）PM

图片来自网络 

我们都知道，在RNA seq 测序的过程中，我们测完序的最终目的是想根据测序的结果，最终分析得到差异基因以及潜在可能的功能分析，那么在进行差异分析以及对表达量进行分析的时候，对基因原始的Count 进行标准化，消除由于测序过程中单个基因自身的长度以及测序深度对数据的影响，是非常关键的一步。

RNAseq 测序，对于一个基因的Count 的计数呢，主要是基于匹配到该基因的外显子上的数目，那么按照这样理解的话，基因越长，比对到该基因(外显子)上的count 数就越多；而影响Count 的另一个因素就是测序深度，也就是该基因在测序的过程中每百万碱基检测到的数目，测序深度越大，那么本次RNA seq 中的所有read count都会增加，因在差异化以及探索表达量的过程中呢，需要对基因长度和测序深度进行标准化，消除这2个因素带来的影响，从而准确的确定基因在样本中是真实的差异表达。

那讲到了这里我们就需要来理解一下在RNA seq 的Count 数进行标准化的常用方法:

常用的方法，包括

1.C（R)PM究竟指什么呢?在常见的分析中，它出现在那里呢?

参考链接(链接：https://cloud.tencent.com/developer/article/1484078  (名称RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同）)

Deseq2 的分析方法 (https://www.jianshu.com/p/bdf2b72b8761  4. edgeR/limma/DESeq2差异基因分析→ggplot2作火山图→biomaRt转换ID并注释)

RPM/CPM: RPM (Reads per million mapped reads)
Calculate Formula:
RPM=Number of reads mapped to a gene *10^6/ Total number of mapped reads from given library

R(C）PM:通过10^6标准化了测序深度的影响，但是没有考虑测序长度的影响。

RPM适合用于产生的read 读数不受基因长度的影响，比如miRNA-seq测序，miRNA的长度一般在20-24个碱基之间。

Deseq2和edgeR 差异分析时，主要涉及CPM数据的归一化

(参考帖子 https://www.jianshu.com/p/2689e9a1d10c DESeq2详细用法)

通常情况下，Deseq2和edgeR进行差异化的分析的时候，都会对数据进行归一化处理，它们的处理方式呢，主要是基于CPM即RPM的分析，然后再通过标准化因子size factor  进行处理，得到一个近似为同方差的值矩阵（沿均值范围具有恒定的方差），这个矩阵可以用于后续的聚类以及PCA 分析以及差异分析。

关于Deseq2

dds <- DESeqDataSet( se, design = ~ cell + dex) 

dds <- DESeqDataSetFromMatrix( countData = countData, colData = colData, design = ~ Group) 

在进行这一步的分析时，其中的design=~batch+condition

设计公式通常格式为~ batch + conditions，batch和conditions都是colData的一列，是因子型变量。为了方便后续计算，最为关注的分组信息放在最后一位。如果记录了样本的批次信息，或者其它需要抹除的信息可以定义在design参数中，在下游回归分析中会根据design formula来估计batch effect的影响，并在下游分析时减去这个影响。这是处理batch effect的推荐方式。在模型中考虑batch effect并没有在数据矩阵中移除bacth effect，如果下游处理时确实有需要，可以使用limma包的removeBatchEffect来处理。

低丰度的数据过滤

dds <- dds[rowSums(counts(dds)) > 1, ] 

在独立筛选（independent filtering）中，DESeq2可以去掉在所有样品中平均表达量CPM不大于min.CPM的基因，以减少假阴性。
EdgeR是保留在2个或更多样品中表达量大于min.CPM的基因。
可以尝试不同的cutoff，以获得最佳效果。

数据的Normolization 的2种方法

Deseq2可以有2种转换数据的方法，其中一种是方差稳定变换，另一种则是正则化对数变换

（1）方差稳定变换，The variance stabilizing transformation

vsd <- vst(object=dds,blind=T) 

样本信息的列名names(colData(vsd))多了1列sizeFactor，colData(vsd)$sizeFactor
基因信息的列名names(rowData(vsd))多了4列
vst函数快速估计离散趋势并应用方差稳定变换。该函数从拟合的离散-均值关系中计算方差稳定变换(VST)，然后变换count data（除以标准化因子），得到一个近似为同方差的值矩阵（沿均值范围具有恒定的方差）。许多常见的多维数据探索性分析方法，例如聚类或PCA，对于同方差的数据表现良好。数据集小于30个样品可以用rlog，数据集大于30个样品用vst，因为rlog速度慢。
 

(2)正则化对数变换，The regularized-logarithm transformation

rld <- rlog(object=dds,blind=F) 

样本信息多了1列sizeFactor，和vsd的sizeFactor相同
基因信息多了7列
rlog函数将count data转换为log2尺度，以最小化有small counts的行的样本间差异，并使library size标准化。rlog在size factors变化很大的情况下更稳健。

blind，转换时是否忽视实验设计。blind=T，不考虑实验设计，用于样品质量保证（sample quality assurance，QA）。blind=F，考虑实验设计，用于downstream analysis。

为什么需要转换?

为什么要转换？为了确保所有基因有大致相同的贡献。
对于RNA-seq raw counts，方差随均值增长，如果直接用size-factor-normalized read counts：counts(dds, normalized=T) 进行主成分分析，结果通常只取决于少数几个表达最高的基因，因为它们显示了样本之间最大的绝对差异。为了避免这种情况，一个策略是采用the logarithm of the normalized count values plus a small pseudocount：log2(counts(dds2, normalized=T) +1)。但是这样，有很低counts的基因将倾向于主导结果。作为一种解决方案，DESeq2为counts数据提供了stabilize the variance across the mean的转换。其中之一是regularized-logarithm transformation or rlog2。对于counts较高的基因，rlog转换可以得到与普通log2转换相似的结果。然而，对于counts较低的基因，所有样本的值都缩小到基因的平均值。用于绘制PCA图或聚类的数据可以有多种：counts、CPM、log2(counts+1)、log2(CPM+1)、vst、rlog等。

关于edgeR中对于数据的归一化(参考帖子  https://cloud.tencent.com/developer/article/1553526  基因芯片数据分析（六）：DESeq2包的基本原理)

edgeR是按照CPM来过滤数据的，通常情况下情况下，阈值设置为1，保留在2个及2个以上样本中大于1 的基因，剩余的基因就被剔除。