前言
NormqPcR这个包用于对qPCR数据进行均一化(normalization，这里需要注意的是，这里的normalization只是使用内参来校正目标基因的表达的意思，与统计学中的归一化(Normalization)，中心化(center)，标准化(scales)之类的术语没有关系)。这篇笔记主要是翻译了这个包的英文文档，本文档以一个小型案例(vignette)来演示了这个包的主要功能，笔记主要内容如下所示：
- 合并技术重复(technical replicates)（了解即可）；
- 处理低Cq值(undetermined values)（了解即可）；
- 计算用于均一化的最佳管家基因(housekeeping genes)（这个最重要，其余的很多函数其实了解即可，这一部分要了解geNorm和NormFinder的算法原理）；
- 使用管家基因来对数据进行均一化(normalization)（了解）。
注：在这篇笔记中，管家基因(housekeeping genes)与内参(reference)指的是同一个意思。
合并技术重复
我们先使用ReadqPCR包来导入原始数据，这个原始数据中含有技术重复，其后缀都是_TechRep.n，其中n表示总的重复数。读取数据过程如下所示：
|
|
计算结果如下所示：
|
|
上面是读取数据的一些信息，从中我们可以发现，里面有技术重复，例如gene_aj_TechReps.1与gene_aj_TechReps.2就是做了2个技术重复。
在我们计算之前，我们还要计算技术重复Cq的算术平均值(arithmetic mean)，这一步我们使用combineTechReps()函数来计算一下，这个函数会生成一个qPCRBatch对象，这样就计算出了技术重复的算术平均值，如下所示：
|
|
结果如下所示：
|
|
从上面我们可以看到：
- 做了2个技术重复；
- 使用
combineTechReps()函数计算了技术重复的算术平均值； - 又用
apply()函数核对了一下，没问题。
处理未检测值
在qPCR实验中，当Cq值超过一定阈值后，qPCR仪有可能会将这个读数标记为Undetermined，也就是数值太小（Cq值越大，表示信号越弱），在qPCRBatch对象中，这个值就会被标记为NA。不同的人有不同的处理策略，常见的处理策略有两种。
第一种就是将那些大于38的Cq值认为是无法检测值，并将其标记为40，这是Taqman软件中的常见做法。而在我们的这个案例中，我们不会将其替换为40，而是替换为NA值，现在我们读入一个Taqman的qPCR数据，这批数据是用96孔板测的，其中，实验组(mia)4个重复，对照组(no-mia)4个重复，如下所示：
|
|
现在我们只看某个样本的检测值，即Ccl20.Rn00570287_m1这个数据，如下所示：
|
|
结果如下所示：
|
|
现在我们使用replaceAboveCutOff()函数将大于35的数值替换为NA，如下所示：
|
|
结果如下所示：
|
|
但是对于有些同学来说，我不想把大于35的数值替换为NA，而是想把所有的NA都替换为40，那么可以使用replaceNAs()函数，如下所示：
|
|
结果如下所示：
|
|
有的时候我们会遇到这样的一样种情况，在一次实验中，例如我们做了96孔板的qPCR实验（做了好几块96孔板），某个孔的检测器中检测的数值大于设定的阈值，而有的孔的检测器检测到的数值则没有大于设定的阈值。有的同学可能就会认为，即然一些孔的检测器检测到的值设为NAs，那么这个孔的检测器在读取其它数值时，也应该设为NAs。因为否则一个检测器的计数出现了异学，那么会导致某个基因表达水平的计算出现错误估计（这一段不太理解，原话如下）：
In addition, the situation sometimes arises where some readings for a given detector are above a given cycle threshold, but some others are not. The user may decide for example that if a given number of readings are NAs, then all of the readings for this detector should be NAs. This is important because otherwise an unusual reading for one detector might lead to an inaccurate estimate for the expression of a given gene.
对于不同的样本类型这个过程应该分开处理，因为你也许想比较不同样本之间的某个特定的基因。因此有必要使用比较矩阵来指定每个样本类型的重复数。因此，在上面的案例中，我们一共有两个类型，分别是实验组与对照组，每组4个重复，我们可以使用contrastMatrix()函数和sampleMaxMatrix()函数来创建比较矩阵，如下所示：
|
|
结果如下所示：
|
|
有关更多的算法细节可以参考limmma包中的内容，那个包中创建差异基因比较的矩阵也是类似方法。
例如，如果用户决定将那些4个重复中，有3个是NA的值都设为NA的话，可以使用makeAllNewVal()函数，如下所示：
|
|
现在我们来看一下Ccl20.Rn00570287 m1这个样本的检测值，对照组的值都为NA，因为原始的4个技术重复中，有3个是NA，有一个是35，因此全设为NA。不过实验组的数据全部进行了保留，如下所示：
|
|
结果如下所示：
|
|
前面这两部分内容（合并技术重复与处理未检测值）了解一下就行，平时也用不到，最重要的是下面的内容。
选择最稳定内参基因
这一部分涉及选择最稳定的内参基因，其中用的是geNorm算法，这个算法的相关文献信息如下所示：
Jo Vandesompele, Katleen De Preter, Filip Pattyn et al. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biology 2002. 3(7):research0034.1-0034.11. http://genomebiology.com/2002/3/7/research/0034/
geNorm
先加载数据，如下所示：
|
|
数据结构如下所示：
|
|
现在我们对选择的内参基因进行排序。我们在这里使用selectHKs()函数来实现geNorm算法中的逐步处理过程，这里的数据都是文献中材料与方法中提到的数据。
第一步，我们计算所有内参基因稳定检测值M后，排除那些M值最高的内参基因；
第二步，再次计算M值，再排除掉最高M值的内参基因，直到剩下两个内参基因（也就是参数minNrHK=2，如下所示：
|
|
计算过程中，这个函数会给出不断计算的M值，结果如下所示：
|
|
以上只是BM这类组织的计算结果。现在我们把其它的组织都计算出来，如下所示：
|
|
我们就获得了基因的下面排序信息，这个就是文献中的Table 3，文献中的信息如下所示：
我们自己的计算结果：
|
|
结果如下所示：
|
|
现在绘制出每个细胞类型的平均基因表达稳定值M的折线图，这个对应于文献中的Figure 2，如下所示：
代码如下所示：
|
|
绘制每个细胞类型的成对方差(pairwise variation)图，这个对应文献中的Figure 3，如下所示：
代码如下所示：
|
|
注：文献中推荐的成对方差值的截止值为0.15，因此低于此值后，就不需要使用其它的内参基因了。
NormFinder
第二种选择内参的基因的方法是NormFinder，这种方法来源于以下文献：
Andersen, C. L., et al. (2004). “Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets.” Cancer Res 64(15): 5245-5250.
现在我们计算Table 3，如下所示：
第一步，导入数据：：
|
|
结果如下所示：
|
|
继续计算，如下所示：
|
|
结果如下所示：
|
|
再看另外的一个数据集Bladder，如下所示：
|
|
计算结果如下所示：
|
|
我们也可以使用geNorm算法来计算，如下所示：
|
|
计算结果如下所示：
|
|
再来看Baldder的数据，如下所示：
|
|
结果如下所示：
|
|
我们也可以使用NormFinder算法的分步计算方式工来看多个内参基因的情况，也就是文献中附录中的Average control gene部分，selectHKs函数中含有这种算法的参数，如下所示：
|
|
结果如下所示：
|
|
在Bladder数据集中，我们发现，排列前面的2个内参基因是HSPCB与RPS13，也就是Figure 1中的内容，如下所示：
|
|
结果如下所示：
|
|
使用多内参进行均一化
单内参的ΔCq方法
单内参均一化的方法就是，使用某个基因的Cq值减去内参的Cq值。NormqPCR中就有相关的计算函数，如下所示：
|
|
这里我们使用了Actb作为内参，计算结果如下所示：
|
|
生成的结果是一个qPCRBatch对象，很多R包例如heatmap就兼容这类对象。
多内参的ΔCq方法
如果使用前面提到的geNorm或NormFinder算法找到了多个内参，此时想要用这多个内参进行均一化，那么就需要计算这几个内参的均值，构成一个“伪内参“(pseudo-housekeeper)进行均一化。现在我们使用三个内参，分别为GAPDH，beta-2-microglobulin和Beta-actin来进行均一化，如下所示：
|
|
结果如下所示：
|
|
单内参2^-ΔΔCq方法
使用2^-ΔΔCq方法哦可以计算两个样本类型之间的相对表达值。默认情况下，同一个样本的目标基因与其内参基因配对进行计算，此时的标准差就体现在不同基因与不同内参上。但是，如果参数中设定了paired=FALSE，那么标准差的计算公式就是s=sqrt(s1^2+s2^2)，这里的s1是目标基因的标准差，s2是内参基因的标准差，但是，如果当你要比较的基因与内参基因是来源于同一个样本时，不推荐这种方法，例如使用Taqman板的情况就是如此，但是，如果要考虑到完整性(completeness)以及要使用不同的生物学重复中获得的内参基因时，可以考虑使用这种计算方法，例如我们最初设计的内参数据并不好时，当采用事后检验(post hoc)时，就采用这种方法，或者是当NormqPCR用于处理那些单独板孔中的扩增数据时，也是采用这种方法。
现在我们读取原始数据，如下所示：
|
|
deltaDeltaCq()函数在进行计算时需要一个比较矩阵。在这个矩阵中，列是样本名（例如case或control），这与limma包中的数据类似。这个比较矩阵的列中只包括1或0这两个数据，用于指定分类变量，如下所示：
|
|
结果如下所示：
|
|
此时，我们可以指定一个内参基因来进行均一化，然后再看一个某个基因在实验组(case)和对照组(control)的比值，计算结果包括以下信息（按列显示）：
- 每个基因的表达值；
- 实验组的ΔCq值：
- Cq值的sd；
- 对照组的ΔCq值；
- 对照组的ΔCq值的sd；
- 2^-ΔΔCq值，这是对照组与实验组的差异；
- 和8. 分别是1和0时的sd值，如下所示：
|
|
结果如下所示：
|
|
也可以使用每个单独的样本/孔(samples/wells)方法来计算平均值。这里我们单独计算sd，然后将它们汇总起来。因此，同一个样本的内参与检测值的本领信息氷是人干活的，这同时会增加变异，如下所示：
|
|
运行结果如下所示：
|
|
多内参2^-ΔΔCq方法
如果想使用大于1个内参进行均一化，例如，我们想使用NormFinder/geNorm算法来进行均一化。我们就可以使用几何平均数(geometric mean)来构建一个伪内参(pseudo-housekeeper)来实现这一目的。现在我们使用GAPDH，beta-2-microglobulin和Beta-actin这三个内参来进行均一化，如下所示：
|
|
计算结果如下所示：
|
|
我们也可以使用那些在样本间方差相同的内参来进行均一化，这个类似于前面提到的第二个deltaDeltaCq方法，如下所示：
|
|
计算结果如下所示：
|
|
NRQ值计算
导入数据，如下所示：
|
|
合并技术重复，并计算SD，如下所示：
|
|
加载数据集的效应(efficiencies)，并将它们添加到数据集中，如下所示：
|
|
计算数据集的NRQ(normalized relative quantities)，我们只考虑两个特征值作为内参基因，如下所示：
|
|
结果如下所示：
|
|
参考资料
- NormqPCR: Functions for normalisation of RT-qPCR data