5. Normalization
Normalization就是为了把两种测序方法得到的数据不一致现象消除，否则如果你通过分析函数知道CpG有差异，你如何知道这个差异是疾病带来的？还是测序方法差异带来的？这一部至关重要，现在已经有很多方法用于完成这个问题，还有科学家在乐此不疲地开发更多算法完成这一步。
在Illumina芯片上，探针有两种不同的设计(称为i型和ii型)，具有不同的杂化化学性质，这意味着来自这两种不同设计的探针将呈现不同的分布。这是一种技术效应，不受i型和ii型探针生物学特性(如CpG密度)差异的影响。i型甲基化分布与ii型甲基化分布最显著的区别在于，ii型甲基化分布呈现出减小的动态范围。
在监督分析中，这可能导致i型探头的选择比ii型探头有偏差。
对于DMR检测，当i型和ii型探头可能落在同一区域时，这个时候，矫正就是很重要的。一般来说，推荐对ii型探针偏差进行调整，我们可以使用champ.norm()函数来执行这种标准化。
在champ.norm()中，我们提供了四种方法来实现这种类型的探针归一化:BMIQ、SWAN、PBC和FunctionalNormliazation。每一种方法都有一些关键的区别，用户可以阅读相关的论文来选择最适合自己分析的方法。
目前，FunctionalNormliazation需要rgSet对象，SWAN需要从IDAT文件生成的mset和rgSet对象(如minfi包描述所述)，FunctionalNormliazation需要rgSet对象。因此，FunctionalNormliazation和SWAN归一化方法仅适用于minfi加载方法的分析。
请注意，BMIQ函数已经更新到1.6版本，它可以为EPIC数组数据提供更好的标准化。
我们注意到，如果样品的甲基化分布明显偏离3-state beta-mixture distribution(这可能发生在甲基化/非甲基化对照中) 或者样品质量差，就可能出现BMIQ不收敛，不能为某些样本提供输出
BMIQ函数的一个新特性是现在可以并行运行，所以如果您的计算机有多个核心，您可以在函数中设置参数“核心”来加速该函数。如果参数PDFplot=TRUE, BMIQ函数的plot将保存在resultsDir中。
进行标准化的代码如下:
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
归一化后，可以使用QC.GUI()函数再次检查结果;记得将QC.GUI()函数中的beta参数设置为myNorm。
QC.GUI(beta=myNorm)
image.png
image.png
规则化之后的结果是显而易见的，type2-BMIQ是规则化之后的分布曲线，和红色的已经一致了，所以测序方法上的差异，不会导致结果出现问题。
SVD
这个分析讲真很简单，但是确实很有用，算法一句话形容就是：把数据进行SVD分解，查看Component有没有与某些重要的Factor相关。SVD大家都知道，主成分分析第一步就是SVD，其性质类似于把重要的Component给找出来了。而与Factor的相关，可以告诉使用者，你的数据中的大部分差异，是你想研究的因素（比如癌症、疾病、年龄等等）导致的，还是由你不想研究的东西，比如样本的种族、数据产出的研究所等等一些无关痛痒的东西决定的。如果是后者，那很不幸你的数据质量也不高……
甲基化数据的奇异值分解方法(SVD)是一种评估和寻找数据集中变异的重要组成部分的数量和性质的一种强大工具。理想情况下，这些变异成分与感兴趣的生物学因素相关，但通常也与技术变异源相关(如批效应)。我们强烈建议用户收集尽可能多的被分析样本的信息(例如杂交日期、采集样本的季节、流行病学信息等)，并在与SVD成分相关时将所有这些因素都包括在内。
如果从.idat文件中加载了样本，当我们在champ.SVD()函数中设置参数RGEffect=TRUE，那么beadchip上的18个内部探针控件(包括亚硫酸盐转换率)将被包括在内, 与ChAMP包的旧版本相比，当前版本的ChAMP . svd()函数将检测所有有效因素进行分析，这意味着该plot包含以下两个条件：
- 协变量与名称、Sample_Name或File_Name没有关联
- 协变量至少包含两个值(例如，要测试作为协变量的“BeadChip ID”，必须包含来自至少两个不同BeadChip的样本)
champ.SVD()函数的作用是:检测pd文件中所有有效的协变量。因此，如果您有一些想要分析的其他的协变量，您可以将它们组合到您的pd文件中，从而允许champ.SVD()测试与数据集中最重要的变化组件之间的关联。
值得注意的是，champ.SVD()只以pd文件的形式使用表型数据。因此，如果您有自己的协变量，希望与当前的pd文件一起分析，请使用cbind()和myLoad$pd合并之后，来一起分析，然后使用这个对象作为champ.SVD()函数的pd参数的输入。
注意，对于不同类型的协变量(分类变量和数字变量)，champ.SVD()使用不同的方法计算协变量和主成分之间的显著性(Krustal.Test and Linear Regression) 。
因此，请确保您的pd对象是一个数据框，数字类型的协变量被转换为数字，分类类型协变量被转换为factor或character类型。
如果你的年龄变量被转成了字符“character”变量，champ.SVD()函数就不会将年龄当成协变量了，必须是数值型，才会被当作协变量分析。
因此，我们必须非常清楚地了解我们的数据集和pd文件。
当我们运行champ.SVD()函数时，所有检测到的协变量都将打印在屏幕上，以便用户检查希望分析的协变量是否正确。生成的结果是与提供的协变量信息最相关的最主要的主成分的热图(保存为SVDsummary.pdf)。在champ.SVD()函数中，我们可以使用isva包中的随机矩阵理论（ Random Matrix Theory）来检测潜在变量的数量。如果我们的方法检测到超过20个组分，我们将只选择前20个。在热图中，较深的颜色代表更显著的p值，表明SVD成分与我们感兴趣的因子之间的相关性更强。
如果SVD分析表明技术因素在变异中占相当大的比例，则有必要实现其他标准化方法(e.g ComBat)这可能有助于消除这种技术上的差异。
ComBat is included in the ChAMP pipeline to remove variation related to the beadchip, position and/or plate, but it may be used to remove other batch effects revealed in the SVD analysis.
您可以使用以下代码来生成SVD图:
champ.SVD(beta=myNorm,pd=myLoad$pd)
image.png
比如上述这批测试数据的结果就非常好：因为数据中的差异，显著与样本的PhenoType互相关联。换言之，你如果做针对Sample_Group（就是Cancer和Normal）的分析，得到的结果就是由于疾病导致的。
额外提一句，如果数据与其他无关的东西相关，而与你想要研究的东西不相关，那么你就需要想很多办法了，比如，如果上述的图片中，红色那个板块（显著相关的意思，越红越相关）是在Array那一行，而不是在Sample_Group那一行，那就意味着你的数据中的差异更多是Array导致的，而Array是芯片中的样本排列方式，换言之就是……你的实验不够好。但有些问题很难避免，比如你的100个病人年龄不同很正常，老人身体上自然会有一些弱化，年轻人可能更好，样本的年龄很自然地会引入一些误差。上述这种情况一旦遇到是非常痛苦的，至今也有很多领域的人在研究如何把这样的错误修正，如何矫正数据等等。ChAMP包集成了一个别人写好的用于修正这一类错误的函数Combat()。个人觉得效果一般，而且时间非常非常非常长（该测试数据只有8个样本，要跑一个半小时！）。但是似乎也没有更好的办法……
如果想要用程序集成的函数，用如下代码修正：
myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Array"))
6. Batch Effect Correction
该函数是为芯片数据开发的ComBat normalization method，sva R包用于实现此功能。ComBat 是明确定义在ChAMP包，以纠正批处理效果。更高级的用户可以使用sva文档实现ComBat，以调整其他批处理效果。您需要设置参数batchname=c("Slide")（根据svd结果来选择校正变量）来在函数中纠正它们。
Note that
champ.runCombat()would extract batch from yourbatchname, so please make sure you understand you dataset perfectly. Another important thing is, previously champ.runCombat() would automatically assume Sample_Group is the covariates need to be analysed, but now, we added a parameter calledvariablenameso that users may assign their own covariates to analysis, thus for some data set whosedisease phenotypeisNOT Sample_Group, please remember to assignvariablenameparameter.
ComBat algorithm采用经验贝叶斯方法对技术变化进行修正。
如果ComBat直接应用于beta值，输出可能不受0到1之间的限制。
因此ComBat adjustment前，ChAMP logit首先转换beta值，然后log转换之后，在进行Combat矫正调整。If the user has chosen to use M-values, please assign logitTrans=FALSE.
在函数champ.runCombat()过程中，程序会自动检测到所有可能被纠正的有效因素，并在屏幕上列出，以便用户检查是否有需要更正的批次。然后，如果用户指定的批处理名称是正确的，函数就会开始进行批处理更正.
在最新版本中，champ.runCombat()还将检查用户的Batch and Variable是否相互混淆。混淆是指你的变量的一个表型的所有样本都来自于你批次的某一个表型，这意味着，如果你批次正确，你的变量包含的信息也消失了，这在SVA包的Combat中是不允许的，这就是为什么我们之前必须要做检查。
你可以使用以下代码进行批次校正：
myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Slide"))
如果您有大量的样本，那么Combat函数可能是一个非常耗时的步骤。
校正后，我们可以使用champ.SVD()函数检查修正后的结果。
7. Differential Methylation Probes（差异甲基化位点）
如果上述步骤都能完成，说明你的数据质量是挺好的，可以分析了！最简单，但同时最重要的一种分析就是差异化CpG，这个步骤的目的就是，找出几百万CpG中的哪些在疾病中发生了变化，而这些变化又是如何导致了基因发生了变化，最终导致了人体生病。而做的方法直观上简单的可怕：
你有100个癌症病人，100个正常人，每个人身体中都有450K个CpG的位点在测序出来，针对其中的每一个CpG，你都有200个数据对不对？如果这一个CpG在100个正常人中和100个癌症病人中的甲基化水平都差不多，你还会继续怀疑它吗，当然不会！但如果，你的100个癌症病人普遍在这一个CpG上的甲基化水平高（不太严谨但是很形象地说，就是DNA那个CpG的序列外层越来越多的部分被甲基附集上去了），而那100个普通人的甲基化水平不高，那这个位点就很有嫌疑了对吧。
为什么需要一定量的样本，比如100个？因为如果你只找两个人，一个德国人一个中国人，那个德国人高，那个中国人矮，你能因此就说德国人在人种上比中国人高吗？当然不能……但如果你找到100个德国人，在找到100个中国人，比较以后，就比较可信了，这涉及t.test()里的power的问题。
差异甲基化位点(DMPs)是最常见的期望输出。
这里我们可以使用champ.DMP()函数计算差异甲基化，也可以使用DMP.GUI()函数检查结果。
champ.DMP() 函数的作用是:使用 limma包，利用线性模型计算差异甲基化的p值。
我们最新的champ.DMP()函数支持年龄等数字变量和包含肿瘤、转移、控制等两种以上表型的分类变量。如果我们的函数检测到年龄等数字变量，我们将对每个CpG位点进行线性回归，以找到与协变量相关的CpGs，例如与年龄相关的CpGs。
当检测到分类变量时，champ.DMP()将对协变量中的两个表型进行对比比较。
Say you have covariate as “tumor”, “metastasis”, “control”, then champ.DMP() would compare “tumor–metastatic”, “tumor-control”, and “metastatis-control”. Each comparison would return one data frame containing significantly differential methylated probes.
champ.DMP()的输出包含两组之间一些p值、t统计量和平均甲基化差异的数据框(这被标记为logFC，因为它类似于基因表达分析中的log fold-change)(仅用于分类共变)。
它还包括每个探针的注释、样本组的平均beta值以及比较中使用的两组的delta值(delta值与logFC相同，我们保留在这里是因为非常老的版本ChAMP使用了它)。
做这一步最讨厌的就是：有时候一把跑出来，几百万位点中，几万甚至十几万都是显著的（做年龄的时候就这样，因为年龄对于人身体的影响太大，可谓是全方位无死角的）
我们可以使用以下代码获取DMPs(这里仍然使用Demo 450K测试数据):
myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)
Remember that in champ.DMP(), each returned result would be stored in an element in a list. Each element in the list is one result of DMPs detected between two phenotypes. If your covariate is numeric variable, the returned data frame name is “NumericVariable”. If your covariate is categorical variable, the returned name would be “PhenoA_to_PhenoB”, where “PhenoA” and “PhenoB” are two name of your phenotypes.
代表的就是PhenoB/PhenoA
在ChAMP的新版本中，我们包含了一个GUI界面，用户可以通过它来检查上述代码生成的myDMP的结果。您只需要提供champ的未修改结果即可。DMP (myDMP)，你用来生成DMP()的原始矩阵和你想要分析的协变。函数的作用是:自动检测协变量是数值型的还是类别型的。如果协变量是数值型的，则DMP.GUI()将为每个CpGs绘制散点图，显示beta值随协变量的变化趋势。如果你的协变量是像癌症/正常这样的分类变量，DMP.GUI()会为显著不同的甲基化CpGs绘制箱线图。您可以使用以下代码打开DMP.GUI()函数。
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)
# myDMP is a list now, each data frame is stored as myDMP[[1]], myDMP[[2]], myDMP[[3]]...
8. Differential Methylation Regions
差异分析是多数研究都要分析的，这里包括三种方法：DMP，DMR，DMB。DMP代表找出Differential Methylation Probe（差异化CpG位点），DMR代表找出Differential Methylation Region（差异化CpG区域，2kb），Block代表Differential Methylation Block（更大范围的差异化region区域，4kb）
简单来说，DMP是找出一个一个的差异甲基化CpG位点，DMR就是一个连续不断都比较长的差异片段，科学家们觉得，这样的连续差异片段，对于基因的影响会更加明显，只找这样的片段，可以使得计算生物学的打击精度更为准确，也可以让最终找出来的结论数据更少，便于实验人员筛选。另外一个类似的东西就是DMB，那个东西出现的原因是，有的科学家觉得，DMR这样的区域还不够显著，DNA上的甲基化出现变化，可能是绵延几千位点的！而且只会在基因以外的区域，但是这些基因以外的区域发生变化，却会导致基因的表达发生变化。你可以想象成，北京周边的河北在大炼钢铁，然后北京也跟着雾霾了，大概就是这意思。
网友评论