🎣®[R]甲基化②ChAMP帮助文档译[下]-甲基化分析流程(The Chip Analysis Methylation Pipeline)

1.甲基化①ChAMP帮助文档译[上]-甲基化分析流程(The Chip Analysis Methylation Pipeline)

6.ChAMP流程说明

如前所述,用户可以选择单独运行每个ChAMP函数,允许与其他流程集成或在每个步骤中保存结果。下面我们进一步详细介绍每个函数。

6.1加载数据

加载数据始终是分析的第一步。ChAMP提供的加载函数是从.idat文件获取数据(其中包含pd(表型)文件(Sample_Sheet.csv))。旧的champ.load()函数利用minfi包从.idat文件中加载数据,但现在我们提供了一种新的“ChAMP”方法来读取数据。两个加载方法都能有效地返回相同的对象,新版本还将返回rgSet和mset 对象,他们可以用作在Champ.norm()中的SWAN和FunctionalNormliazation方法的输入。默认情况下,加载的是当前工作目录中的数据,当前目录中您应该含有所有.idat文件和sample sheet。sample sheet 需要是一个CSV文件。champ.load() 用于加载idat文件中的数据并自动过滤掉可能影响分析结果的SNP。SNP列产生于Zhou’s Nucleic Acids Research paper in 2016 (http://zwdzwd.github.io/InfiniumAnnotation),专为450K和EPIC设计,它也考虑到人口差异。因此,对于450k珠阵列数据,在过滤低质量探针之前,数据集将包括485,512个探针。而对于EPIC珠阵列数据,在过滤探测数据集之前,将包括867,531个探针。

目前,加载的默认方法是“ChAMP”。您可以在champ.load()中设置“method”来使用经典的minfi方式。请注意,champ.load()中的“ChAMP”方法仅仅是champ.import()和champ.filter()的组合。如果您对champ.load()返回的结果不满意,想获取beadcount,detection P matrix, Meth UnMeth matrix,也可以使用champ.import()来生成它们。

用户可以调用以下代码来加载数据集:

myLoad <- champ.load(testDir)

或者,用户也可以从ChAMPdata包中的数据集对象testDataSet加载数据:

data(testDataSet)

检查pd文件(Sample_Sheet.csv)很重要的一件事。大多数pd文件格式相似,但是它们最有价值的信息并不总是存储在同一列中,例如,对于大多数pd文件,Sample_Group表示要研究的表型组,但是在某些pd文件中,可能存储类似的信息在名为“Diagnose”或“CancerType”的列中,因此用户必须很好地了解其pd文件以实现成功的分析。

myLoad$pd
##    Sample_Name Sample_Plate Sample_Group Pool_ID Project Sample_Well
## C1          C1           NA            C      NA      NA         E09
## C2          C2           NA            C      NA      NA         G09
## C3          C3           NA            C      NA      NA         E02
## C4          C4           NA            C      NA      NA         F02
## T1          T1           NA            T      NA      NA         B09
## T2          T2           NA            T      NA      NA         C09
## T3          T3           NA            T      NA      NA         E08
## T4          T4           NA            T      NA      NA         C09
##     Array      Slide                   Basename                  filenames
## C1 R03C02 7990895118 Test450K/7990895118_R03C02 Test450K/7990895118_R03C02
## C2 R05C02 7990895118 Test450K/7990895118_R05C02 Test450K/7990895118_R05C02
## C3 R01C01 9247377086 Test450K/9247377086_R01C01 Test450K/9247377086_R01C01
## C4 R02C01 9247377086 Test450K/9247377086_R02C01 Test450K/9247377086_R02C01
## T1 R06C01 7766130112 Test450K/7766130112_R06C01 Test450K/7766130112_R06C01
## T2 R01C02 7766130112 Test450K/7766130112_R01C02 Test450K/7766130112_R01C02
## T3 R01C01 7990895118 Test450K/7990895118_R01C01 Test450K/7990895118_R01C01
## T4 R01C02 7990895118 Test450K/7990895118_R01C02 Test450K/7990895118_R01C02

正如我们所见,Sample_Group包含两个表型:C和T。在该数据集中,C表示“Control”,T表示“Tumor”。

6.2过滤数据

ChAMP提供了一个全面的过滤函数champ.filter(),可以将任何数据矩阵(beta,M,Meth,UnMeth,intensity)作为输入,并对其进行过滤。 champ.filter()不需要任何特定的对象或类结果,也不需要一个IDAT文件,只要你有一个的beta矩阵,就可以完成过滤。在champ.filter()中,由于某些探针和样本可能由于质量较差而被去除,所以在过滤数据中可能仍然存在NA。因此,如果仍然存在NAs ,我们在函数中提供了一个参数autoimpute来对探针进行补缺失值。实际上,补缺失由champ.impute()函数完成,读者可以学习该函数了解更多信息。

对于某些过滤方法,必须提供附加的数据,例如,如果要使用检测P值执行过滤,则需要提供检测P值矩阵。另外,如果要根据珠数(beadcounts)进行筛选,则需要提供珠粒(beadcount)数据。如果您使用champ.load()函数的默认方法“ChAMP”读取IDAT文件,则返回结果中应包含珠粒(beadcount)信息。另外,champ.filter()被用为对同一数据集的多种类型的数据框(data frames)进行过滤,例如,您有同一.idat文件的beta矩阵,M矩阵,intensity 矩阵,您可以使用champ.filter()同时进行过滤。

  • 首先,您必须输入至少一个数据矩阵。

  • 第二,如果你提供尽可能多的输入矩阵,它们必须有完全相同rownames和colnames。例如,如果要同时对intensity 数据和M值进行过滤,则必须确保它们具有相同的名称,否则ChAMP将假定它们是来自不同来源的数据框。那么,champ.filter()就会终止过滤。

  • 由于在过滤过程中可能会删除低质量的样本(例如具有许多未检测到的探针的样本),因此确保您的pd文件中的Sample_Name标识符与数据矩阵的列名完全相同是很重要的

  • 第四,如果要进行插补,则必须1.提供检测P矩阵,β或M矩阵,2.并且参数ProbeCutoff不能为0,ProbeCutoff代表NA比率的阈值,3.探针的NA比率超过阈值则删除该探针。如果这三个条件中的任何不满足,那么champ.filter()会将autoimpute参数重置为FALSE,那么插补过程将不会被完成。

  • 第五,如果你想通过珠粒(beadcount)过滤,你必须提供珠粒(beadcount),champ.import()函数将返回beads信息。或者,您也可以使用wateRmelon包中的beatcount函数来提取“minfi”方法的rgSet对象。

champ.filter()函数是champ.import()的后续过程,它是一个由ChAMP团队编写的新的加载函数。你可以像下方所示使用它:

myImport <- champ.import(testDir)
myLoad <- champ.filter()

上面的代码中的champ.filter()的结果应该和旧的minfi方法产生的结果相同。

实际上在加载过程中,champ.load()函数也会做一些过滤。以下是在champ.load()和champ.filter()中执行的过滤步骤。

  • 第一个过滤器用于过滤p值(<0.01)的探针。这利用了存储在.idat文件中的检测p值。champ.import()将读取这些信息并将其转换成数据框(frame)。对于每个探针,如果p值小于0.01,则将其视为失效探针。之后将会将failedSamples结果打印到屏幕上:显示每个样本失效的探针的比例。如果这些值中的任何一个值高(> 0.1),您可以考虑从分析中删除该样本。以前,我们发现在许多情况下,大数据集中只有一个或两个样本不能用于分析,它们可能具有70%甚至80%的失效探针。因此,如果我们仅对探针进行过滤,则大约80%的探针将被掩蔽,因此我们开发了一种新的过滤系统,既可以对样品也对探针质量进行过滤。如果某个样本的探针失效率高于特定阈值(默认= 0.1),那么该样品将被丢弃,并对剩余的样品进行探针的过滤。样本和探针的阈值可以通过参数SampleCutoff和ProbeCutoff分别控制。

  • 第二,ChAMP将在每个探针的至少5%的样本中过滤出带有< 3颗珠子的探针

    可以使用filterBeads参数更改此默认值,也可以使用beadCutoff参数调整频率。

  • 第三,ChAMP将默认筛选出数据集中包含的所有非CpG探针。

  • 第四,默认情况下,ChAMP将过滤所有与SNP相关的探针。SNP列表来自 Zhou’s Nucleic Acids Research paper in 2016(http://zwdzwd.github.io/InfiniumAnnotation) 。请注意,如果您知道您的数据来自哪个群体,您可以选择特定群体进行筛选。否则,ChAMP将使用Zhou提供的一般推荐探针进行过滤。您只需要分配“population”参数即可实现此目的。

  • 第五,默认情况下,ChAMP将过滤所有multi-hit探针。multi-hit探针列表来自于Nordlund的2013年“基因组生物学报” https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-9-r105

  • 第六,ChAMP将过滤出染色体X和Y中的所有探针。这也是默认设置,但用户可以使用filterXY参数进行更改。

所有上述过滤由champ.load()和champ.filter()函数中的参数控制,用户可以根据需要进行调整。请注意,尽管大多数ChAMP功能都支持单独的beta矩阵分析,即不依赖于从头开始处理.idat文件,champ.load()不能单独完成beta矩阵的过滤。对于没有.IDAT数据但却拥有beta矩阵和Sample_Sheet.csv的用户,您可以使用champ.filter()函数执行过滤,然后使用后续的函数进行分析。

除了过滤功能之外,我们还提供了一个函数CpG.GUI(),供用户检查CPGs其在染色体,CpG岛,TSS重复上的分布。例如,该函数可以被用于任何CPGs list,比如您从DMR中获得有意义的CpGs 列表,您就可以使用以下函数来检查CpG向量的分布情况:

CpG.GUI(CpG=rownames(myLoad$beta),arraytype="450K")
image.png

这是一个用于展示您的CpG列表的分布的有用函数。如果您获得DMP列表,您就可以使用此函数来检查DMP在染色体上的分布情况。

6.3进一步的质量控制和探索性分析

质量控制是确保数据集适合下游分析的重要过程。champ.QC()函数和QC.GUI()函数将为用户绘制一些图来轻松的检查其数据质量。通常情况下,将生成三个图:

QC.GUI(beta=myLoad$beta,arraytype="450K")
image.png

image.png

image.png
  • mdsPlot(多维尺度分析图):该图展示在所有样本中top 1000个最大可能突变的探针对应的值来可视化样本间的相似性。样本在该图中根据样品组着色。

  • densityPlot:每个样本的beta分布; 用户可以使用该图来查找与其他样本显着偏差的样品,这样的样品可能不具有良好的质量(例如不完全亚硫酸氢盐转化)。注意:如果对照样品已包含在研究中,densityPlot也被用于鉴定和确认甲基化或非甲基化对照样品)。

  • 树形图:是对所有样本的聚类图,您可以选择不同的方法来生成该图。在champ.QC()函数中有一个参数Feature.sel="None"。而“无”表示样本之间的距离将由所有探针直接计算(如果您直接在大数据集上执行此操作,则您的服务器可能会崩溃),“SVD”表示champ.QC()函数将使用SVD方法在beta矩阵上进行解卷积,然后用“isva”包中EstDimRMT()方法挑选重要组件。然后根据这些组件计算样本之间的距离。

或者您可以使用如下的QC.GUI()函数,但这可能是挺占内存的,所以当您在ChAMP中运行此GUI功能时,请确保您具有良好的服务器或计算机。

QC.GUI(beta=myLoad$beta,arraytype="450K")

RR

与champ.QC()相反,QC.GUI()将提供五个交互图。图片如下:mdsPlot,I型和II型密度图,样本β分布图,树状图和前1000个最可变CpG的热图。用户可以轻松地与这些图进行交互,以查看样品是否有足够质量进行进一步分析,并检查聚类结果和top CpGs。

  • Probe-Type-I&II图:您的数据集中的Type-I和Type-II探针的beta分布,可以帮助您检查数据集的归一化状态。

  • Top variable CpGs’ heatmap:该热图将选择最可能变化的CpG,并根据这些位点的甲基化值创建热图。

6.4标准化

在Illumina beadarrays上,探针有两种不同的设计(称为I型和II型),具有不同的杂交化学特性,这意味着来自这两种不同设计的探针的值将表现出不同的分布。这种差异是由I型和II型探针的技术因素引起的,与生物学特性(如CpG密度)所引起的变化无关。I型和II型甲基化分布之间最显着的差异在于II型分布表现出较小的动态范围。在监督分析中,这可能导致I型与II型探针的选择有偏差。对于DMR检测,I型和II型探针可能落在同一区域(显然对此进行调整至关重要)。一般来说,建议对II型探头偏置进行调整,您可以使用champ.norm() 来实现这一调整。

在champ.norm()中,我们提供了四种方法来进行这种II型探针归一化:BMIQ ,SWAN ,PBC 和Functional Normliazation 。每种方法都有一些关键的区别,用户可以阅读相关文章,选择最适合的分析方法(BMIQ ,SWAN ,PBC 和Functional Normliazation的参考文献分别是帮助文档参考文献的第 3 、1、2、4 篇)。

目前,FunctionalNormalization需要rgSet,SWAN需要从IDAT文件生成的mset和rgSet对象(如minfi包描述中所述)。因此,FunctionalNormliazation和SWAN归一化方法仅适用于“minfi”加载方法的分析。

请注意,BMIQ功能已更新到版本1.6,这将为EPIC阵列数据提供更好的标准化。我们注意到,BMIQ可能不会为某些样品提供输出(如果样品的甲基化分布明显偏离3-state beta-mixture 分布,或者因为样品质量差)。BMIQ函数现在可以并行运行,因此如果您的计算机有多个内核,您可以在函数中设置参数“cores”来加速函数的运行。如果设置参数PDFplot=TRUE,BMIQ函数的图将保存在输出路径上。

执行标准化的代码如下:

myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)

标准化后,您可以使用QC.GUI()函数再次检查结果; 记住要将QC.GUI()函数中的beta参数设置为“myNorm”。

6.5 SVD图

由Teschendorff (参考文献第七篇) 实现的奇异值分解方法(SVD)是评估甲基化数据中重要变异组成部分的数量和性质的有力工具。这些变异组分可能与感兴趣的生物学因素相关,但通常也与技术因素(例如批量效应)相关。我们强烈建议用户尽可能收集所分析样品的信息(例如杂交日期,收集样本的季节,流行病学信息等),并在做SVD相关分析时包括所有这些因素。要是你的样本已从.idat文件加载,在设置参数RGEffect = TRUE时,champ.SVD()函数会对beadchip上的18个内部探针进行controls(包括亚硫酸氢盐转换效率)。

与旧版本的ChAMP包相比,当前版本的ChAMP . svd()将检测所有有效的因素进行分析,这意味着该plot包含两个以下条件

  • 协变量不与name、Sample_Name或File_Name相关联

  • 协变量至少包含两个值(for example, for ‘BeadChip ID’ to be tested as a covariate, samples from at least two different BeadChips must be included in the study)

champ.SVD()函数检测您的pd(表型)文件中的所有有效协变量。如果您有一些需要分析的独特协变量,您可以将它们整合到您的pd文件(表型文件)中,让champ . svd()来测试与数据集中最显著的变化组件的关联。请注意,champ.SVD()仅使用pd文件形式的表型数据。所以,如果你有你自己的协变量,你想用你当前的pd文件进行分析,请将cbind()你的协变量和myLoad $ pd一起使用,然后使用这个对象作为champ.SVD()函数的pd参数的输入。

注意,对于不同类型的协变量(分类和数字),champ.SVD()使用不同的方法来计算协变量和组分之间的significance (Krustal.Test和线性回归)。因此,请确保您的pd对象是数据框(data frame),数字协变量被转换为“numeric”类型,而分类协变量被转换为“factor”或“character”类型。如果您的年龄协变量被指定为“character”,则Champ.SVD功能将不会检测(因为它应该要被指定为numeric分析)。因此,我们建议用户非常清楚他们的数据集和pd文件。

当champ.SVD()正在运行时,所有检测到的协变量将被打印在屏幕上,以便用户可以检查您想要分析的协变量是否正确。结果是与协变量信息相关的主成分的热图(保存为SVDsummary.pdf)。在champ.SVD()中,我们使用来自于isva包的随机矩阵理论(Random Matrix Theory)来检测潜在变量的数量。如果我们的方法检测到20个以上的组件,我们将仅选择前20个。在热图中,较暗的颜色代表更显著的p值,表明SVD分量与感兴趣因素之间的相关性更强。如果SVD分析显示技术因素占变异的很大一部分,那么有必要使用其他归一化方法(如ComBat),可以帮助消除这种技术变异。

您可以使用以下代码生成SVD图:

champ.SVD(beta=myNorm,pd=myLoad$pd)
image.png

上图是在我们的示例的 450K阵列数据集上生成的。由于该数据集中只有8个样本,所以最终的headmap看起来并不复杂。下面我们使用GSE40279 生成的另一个图,这个数据包含656个样本和一个复杂的pd文件,包括年龄和种族等协变量。下面我们将展示使用champ.SVD()函数生成的SVD图的结果。

image.png

在上图中,Age是一个数字变量,而其他因素是分类因素。较暗的块意味着去卷积分量和协变量之间的相关性更强。

6.6批效应校正

这个函数实现了ComBat归一化方法。sva 这个R包用于实现此功能。ComBat在ChAMP包中被明确定义用以矫正批次效应。进阶用户可以查阅sva帮助文档实现ComBat对其他批次效应的矫正。您需要在函数中设定参数batchname=c("Slide")以调整它们。

请注意,champ.runCombat()将从批次名称中提取批次,因此请确保您完全了解你的数据集。先前版本的champ.runCombat()会自动地设定Sample_Group 为需要被分析的协变量,但是现在我们添加了一个参数variablename,以便用户可以将自己的协变量用于分析,因此,对于一些数据集,其疾病表型是非 Sample_Group,请记住设定variablename参数。

ComBat算法使用经验Bayes方法来校正技术变量。如果ComBat直接应用于beta值,输出可能不在0和1之间。由于这个原因,ChAMP 在做ComBat调整之前先logit变换beta值,然后在ComBat调整之后计算反向logit转换。如果用户选择使用M值,请分配logitTrans=FALSE。

在函数champ.runCombat()中,程序将自动检测所有将被矫正的有效因素并在屏幕上列出它们,以便用户可以检查是否存在要纠正的批次。然后如果用户分配的批处理名称正确,则该函数将开始批量修改。

在最新版本中,champ.runCombat()也会检查用户的批次和变量是否相互联系。 Confounding意思是 你的变量的一个表型的所有样本都来自于你批次里的特定一个表型 ,这意味着,如果您校正批次,您的变量包含的一些信息也将失去,这种情况在SVA包的Combat函数中是不允许出现的,这就是我们事先要做检查的原因。

您可以使用以下代码进行校正批次:

myCombat <- champ.runCombat(beta=myNorm,pd=myLoad$pd,batchname=c("Slide"))

如果您样本量大,执行ComBat函数可能是一个非常耗时的过程。校正后,您可以使用champ.SVD()函数来检查你的校正结果。

6.7差示甲基化探针

差异化甲基化探针(DMP)是甲基化分析中最常见的输出。用户可以使用champ.DMP()函数来计算差异甲基化,并且可以使用DMP.GUI()来检查结果。

champ.DMP()函数使用limma包计算差异甲基化的p值(用到线性模型)。我们最新的champ.DMP()函数将支持数字变量,如年龄,和包含两个以上表型的分类变量,如“tumor”,“metastasis”,“control”。如果我们的函数检测到像年龄这样的数字变量,则会对每个CpG位点进行线性回归,找到您的协变量相关CpGs,将之看做年龄有关的CpGs。当检测到分类变量时,champ.DMP()将在您的协变量中做两两表型之间的对比。假设你有“tumor”,“metastasis”,“control”这些协变量,那么champ.DMP()将比较“tumor–metastatic”, “tumor-control”, and “metastatis-control”。每个比较将返回一个包含显着差异甲基化探针的数据框(data frame)。

champ.DMP()的输出包括两组之间的P值、t统计量和平均甲基化差异(被标记为logFC,(因为他和基因表达分析中的log fold-change很相似 ))的数据框。它还包括每个探针的注释,样本组的平均β值和比较中使用到的两组的delta beta值(delta值与logFC类似,我们在这里保留它是因为旧版本ChAMP有使用到它)。

在我们最新版本的champ.DMP()中,如果您有年龄或时间序列数据的协变量,numeric variables are fully support not。然而,更进阶的用户可以使用limma的其他设置运行函数,并使用DMR hunter的输出(包括所有探针及其p值)。

用户可以使用以下代码获取DMP(我们仍然使用的是示例的450K测试数据):

myDMP <- champ.DMP(beta = myNorm,pheno=myLoad$pd$Sample_Group)

我们可以检查champ.DMP()的结果(如下)。记住,在champ.DMP()中,每个返回的结果将存储在列表中的元素中。列表中的每个元素是在两种表型之间检测到的DMPs的一个结果。如果您的协变量是数字变量,则返回的数据框名称为“NumericVariable”。如果您的协变量是分类变量,返回的名称将为“PhenoA_to_PhenoB”,其中“PhenoA”和“PhenoB”是您的两种表型名称。

head(myDMP[[1]])

在新版本的ChAMP中,我们提供了一个GUI界面,用于检查从上述代码生成的myDMP的结果。您只需要提供champ.DMP(myDMP)的“unmodified”结果、您用于生成DMP的原始beta矩阵和您要分析的协变量。DMP.GUI()函数将自动检测您的协变量是数字还是分类。如果您的协变量是数字,DMP.GUI()将绘制每个CpG的散点图,显示beta值的趋势以及您的协变量。如果您的协变量是癌症/正常的分类变量,则DMP.GUI()将绘制显着差异甲基化CpG的boxplot。您可以使用以下代码打开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]]...
image.png

左侧面板包含了卡CPGS显著性的阈值参数。目前仅支持两个参数:P值和abs(logFC)。用户可以设置这两个值,然后按“Submit”在右侧生成显显著的CpG表。右侧窗口是显著CpGs表,所有信息都包含在内,表格可以按列排列。

image.png

第二个选项卡提供了您通过控制左侧面板的cutoff值生成的所有显著CpGs的热图。请注意,此热图的行和列已经进行了层次聚类。

image.png

第三个面板提供了每个CpG feature 占比的barplot(例如TSS200,TSS1500,body,opensea,shore),以及hyper和hypo CpGs的cgi摘要信息。您可以点击图例中的标记来关闭或打开某些bar。

第四个选项卡左侧面板上由Gene这个参数控制。您可以选择一个基因名称(默认为NFIX),然后按“submit”。然后这个基因的plot图将被绘制在上方的图中。请注意,这里的x轴不是每个CpG的实际MAPINFO。每个CpG之间的距离是相同的。每种表型有两种线条,实线代表平均值图,虚线表示Loess line,用户可以通过点击图例上的标记来关闭或打开它们。以下是每个CpG的功能和cgi信息。在该标签的底部是从维基提取的基因信息,也是由左侧的基因参数控制。

image.png

第五个选项卡包含两个图。顶部是通过左侧cutoff参数选择的显著CpGs的前70个基因的富集图。每个bar中的数字表示在该基因中包含多少高或低差异的甲基化CpG。用户可以使用这个图来找到感兴趣的基因,然后使用第4个选项卡绘制这些基因。在此选项卡的底部,是由左侧面板上的参数控制的某些CpG的boxplot。

请注意,这里的所有图可以放大或缩小。并可直接下载。您甚至可以更改浏览器的形状以调整绘图大小,然后以不同的分辨率下载。

6.7.1 DMP.GUI对数值变量的影响

以上DMP.GUI图是由我们的示例数据450K阵列数据生成的。由于这是一个分类协变量数据集,下面我们将展示数值协变量数据集的部分结果:GSE40279。在这里,我们将做完整的分析,不过只展示了GSE40279数据集中的champ.DMP()检测到的与Age相关的CpGs的DMP.GUI()网页。请注意,在加载,过滤,normalization之后,我们得到以下结果,然后运行champ.DMP()和DMP.GUI()。

对于数字变量,为了显示beta值随变量变化的趋势,DMP.GUI()会自动将数字变量逐个切换成几个组(默认为4组),然后根据每组绘制线和点。如果你的变量是年龄,DMP.GUI()会自动将您的年龄变量剪切成几个组,如30-,30-50,50-70,70+,然后绘制四条lines和几组点。这样,您就可以看到旧样本组的线与新样本组的线的显著不同。

首先显示DMP.GUI()网页的第三个选项卡(基因相关的选项卡):

image.png

在上面的plot中,我们可以看到,DMP.GUI()自动将年龄协变量分为4组。然后对于每个组,绘制蚂蚁线(dots ant lines)。我们可以看到,在基因ELOVL2中,总共有8个CpG显示年龄相关的趋势。线条显示,老年组(80.5〜101)倾向于比较年轻的群体高甲基化。

接下来我们展示DMP.GUI()网页的第四个选项卡(CpG相关选项卡):

image.png

在上图中,CpGs将被绘制到散点图中,这样我们可以看到CpG随着年龄的增长趋势。

6.7.2羟甲基化分析(Hydroxymethylation Analysis)

有些用户可能想要分析羟甲基化数据,这时可以用champ.DMP()函数来完成,这里我们在下方提供了一些演示代码,供用户检查。

myDMP <- champ.DMP(beta=myNorm, pheno=myLoad$pd$Sample_Group, compare.group=c("oxBS", "BS"))
# In above code, you can set compare.group() as "oxBS" and "BS" to do DMP detection between hydroxymethylatio and normal methylation.

hmc <- myDMP[[1]][myDMP[[1]]$deltaBeta>0,]
# Then you can use above code to extract hydroxymethylation CpGs.

6.8差异甲基化区域

Differentially Methylated Regions (DMRs) are extended segments of the genome,它展示了两组之间DNA甲基化水平的定量变化。champ.DMR()是ChAMP提供的一个函数,用于计算并返回将探针绑定到离散的甲基化区域的数据框,并带有P值。

ChAMP内有三种DMR算法:Bumphunter,ProbeLasso和DMRcate。Bumphunter算法首先将所有探针分组成小簇(或区域),然后应用随机置换方法来估计候选DMRs。这种方法非常用户友好,不依赖于之前函数的任何输出。Bumphunter中的置换步骤可能是计算密集型的,不过用户可以分配更多的内核来加速它。Bumphunter算法的结果是包含所有检测到的DMR的数据框,包含其长度,簇,被注释的CpG数量。

对于ProbeLasso方法,最终的数据框是从探针的边缘输出(limma output)及其关联统计数据中得出的。在此之前,必须将champ.DMP()的结果输入到champ.DMR()函数中才能运行ProbeLasso,但现在我们已经升级了该功能,而且ProbeLasso不需要来自champ.DMP()的结果。对于ProbeLasso功能的原理和机制,用户可以在帮助文档的第11篇参考文献中查询。

最近DMRcate方法被加入到ChAMP中。这是一种数据驱动的方法,该方法对所有的注解都是不可知的(除了对空间的,尤其是染色体坐标)。(This is a data-driven approach that is agnostic to all annotations except for spatial ones, specifically chromosomal coordinates)。DMRcate 能够鲁棒的评估个体化CPG位点(从limma中获得)的差异甲基化。DMRcate通过每个450K探针计算出的主观t统计量的平方来实现DMR-finding功能。 Then DMRcate would apply a Gaussian kernel to smooth this metric within a given window, and also derive an expected value of the smoothed estimate (in other words, one with no experimental effect) from the varying density of CpGs incu by reduced representation and irregular spacing.有关DMRcate功能的更多细节,用户可以阅读帮助文档中的参考文献第13篇和DMRcate R pacakge。

请注意,在这个版本的ChAMP中,我们在champ.DMR()中进行了一些严格的过程检查,现在,champ.DMR()仅采用刚好具有两种表型的分类协变量。如果您的协变量包含多个表型,称为“肿瘤”/“转移”/“控制”,请手动选择其中任意两个,并输入相应的β 矩阵和pheno信息。

用户可以使用以下代码来生成DMR:

myDMR <- champ.DMR(beta=myNorm,pheno=myLoad$pd$Sample_Group,method="Bumphunter")

Bumphunter,DMRcate和ProbeLasso函数的输出都是数据框,但可能包含不同的列。在这里我们可以检查DMRcate的结果。

head(myDMR$DMRcateDMR)
##       seqnames     start       end width strand no.cpgs        minfdr
## DMR_1     chr2 177021702 177030228  8527      *      48  3.243739e-69
## DMR_2     chr2 177012117 177017797  5681      *      33  1.137990e-59
## DMR_3    chr12 115134148 115136308  2161      *      51 3.571093e-108
## DMR_4    chr11  31820441  31828715  8275      *      40  3.087712e-43
## DMR_5    chr17  46655289  46658170  2882      *      27  1.565729e-85
## DMR_6     chr6  32115964  32118811  2848      *      50 5.055382e-123
##           Stouffer maxbetafc meanbetafc
## DMR_1 4.011352e-24 0.6962753  0.4383726
## DMR_2 6.555658e-20 0.6138975  0.4676655
## DMR_3 5.711951e-19 0.5372083  0.3601741
## DMR_4 1.206741e-18 0.5673115  0.4112588
## DMR_5 2.087453e-17 0.6859260  0.4812521
## DMR_6 7.216567e-16 0.5772168  0.2971849
##                                                                                              overlapping.promoters
## DMR_1                                                                       HOXD3-001, HOXD3-004, RP11-387A1.5-001
## DMR_2                                                                             HOXD4-001, MIR10B-201, HOXD3-005
## DMR_3                                                                                                         <NA>
## DMR_4 PAX6-016, PAX6-006, PAX6-014, PAX6-015, PAX6-007, PAX6-028, PAX6-025, PAX6-027, PAX6-029, PAX6-024, PAX6-026
## DMR_5                                                      HOXB4-001, MIR10A-201, HOXB3-009, HOXB3-005, MIR10A-001
## DMR_6                                                        PRRT1-002, PRRT1-201, PRRT1-007, PRRT1-006, PRRT1-003

就像DMP结果一样,我们提供DMR.GUI()DMR结果的接口。所有来自champ.DMR()的三种不同方法的结果都可以在GUI函数中被接受。DMR.GUI()函数会自动检测到DMR结果是什么样的输入。因此,不建议修改myDMR结果的结构。

DMR.GUI(DMR=myDMR)
# It might be a little bit slow to open DMR.GUI() because function need to extract annotation for CpGs from DMR. Might take 30 seconds.

R
image.png

就像上面的DMP.GUI一样,左边有一个参数面板,用户可以设置参数来选择 significant DMRs 。我们目前只提供两个参数:P值和最小数目探针。上面的标题将显示您用来获取DMR结果的方法(这里是DMRcate)。注意,对于不同的方法,P值可以指示不同的列。Bumphunter是pvaluearea,而ProbeLasso是dmrP,并且没有DMRcate结果的cutoff值(这在计算过程中已经完成)。 参数选择完后,按“Submit”,您会得到右侧的significant DMR表格。对于不同的方法,表内容可能不同。因此,用户应该考虑您想要使用哪种方法来查找DMR。

image.png

第二个选项卡提供了在这些DMRs中包含的所有CpG的表,可以将其视为CpG和DMR之间的映射列表。另外如果您不能在第一个选项卡中获取基因信息,你可以在这里获得。在该表中,对于每个CpG,其相应的基因和DMR被标记。

image.png

第三个图是每个DMR的图,就和DMP.GUI()函数中的gene plot图类似。这里x轴不是每个CpG的真实MAPINFO,每个CpGs之间的距离是相等的。Also I have feature and CpGs marked under the plot。在此选项卡的顶部,是该DMR中所有CpG的表。您可以使用左侧面板上的DMR index参数来选择要检查的DMR。

image.png

第四选项卡中提供了基因富集信息和热图。请注意,在这个plot中可能有许多基因和CpGs。所以如果你选择的DMR太多,那么这个plot可能太大了,不能打开。但是你不用担心看不清楚,因为你可以放大看细节。

目前,在一个协变量中,champ.DMR()不支持多个表型,理论上我们可以做配对比较,但是需要很长时间。另外,我们不支持数字变量,比如说寻找“年龄相关的甲基化区域”,它要求很多繁重的编码工作,但我们可能会在将来的版本中支持。

6.9差异甲基化块

识别差异甲基化块(DMB)也可能有潜在的意义(potential interest).这里我们提供一个推断DMBs的函数。在我们的Block-finder函数中,champ.Block()函数将首先基于它们在基因组上的位置来计算全基因组上的小簇(区域)。然后,对于每个簇(区域),将计算其平均值和平均位置,因此您可以假设每个区域将被折叠(collapsed )成单个单元。When we finding DMB, only single unit from open sea would be used to do clustering。这里将使用Bumphunter算法在这些区域(被折叠后的单个单位)上找到“DMR”。在我们以前的论文(帮助文档第23篇)和其他研究工作中(帮助文档第24篇),我们证明差异甲基化可能在各种癌症中显示出普遍的特征

近年来,有关差异甲基化块的研究较多。这里我们提供一个计算甲基化区块的函数。在我们的Block-finder函数中,champ.Block()将首先基于它们在基因组上的位置来计算全基因组上的小簇(区域)。然后对于每个簇(区域),其平均值和平均位置将被计算,因此您可以假设每个区域将被折叠成一个点。然后,Bumphunter算法将用于在这些区域(dots after collapse)中找到“DMR”。请注意,Note that only open sea regions will be used to calculate Blocks.。在我们之前的论文中(帮助文档第23篇),我们证实了差异甲基化块可以显示各种癌症的普遍特征。(译者注:感觉这段和上段有点重叠)

用户可以使用以下代码生成甲基化区块:

myBlock <- champ.Block(beta=myNorm,pheno=myLoad$pd$Sample_Group,arraytype="450K")

champ.Block()函数返回很多对象,但最重要的是第一个列表:Block。您可以查看如下:

head(myBlock$Block)
##         chr     start       end      value     area cluster indexStart
## Block_1  12 125405786 133245207 -0.1398238 61.66232     649      27996
## Block_2   6  28839951  30308328 -0.1725335 40.20032     356      81359
## Block_3  15  24778439  27360551 -0.1752729 38.73531     732      35251
## Block_4  14 100851839 102144080 -0.2166707 35.75066     724      34607
## Block_5  11 130933979 133954227 -0.1833729 34.29073     616      22569
## Block_6   7  40026390  42107777 -0.2616146 24.85339     400      88773
##         indexEnd   L clusterL      p.value  fwer  p.valueArea fwerArea
## Block_1    28436 441     2176 1.698689e-05 0.028 1.698689e-05    0.028
## Block_2    81591 233     2396 1.698689e-05 0.028 1.577354e-04    0.136
## Block_3    35471 221      261 1.698689e-05 0.028 2.196162e-04    0.136
## Block_4    34771 165      708 1.698689e-05 0.028 3.858451e-04    0.188
## Block_5    22755 187     1601 1.698689e-05 0.028 4.622861e-04    0.188
## Block_6    88867  95     1497 1.698689e-05 0.028 1.571287e-03    0.420

并且我们提供一个GUI功能,供用户查看上面的myBlock的结果。

Block.GUI(Block=myBlock,beta=myNorm,pheno=myLoad$pd$Sample_Group,runDMP=TRUE,compare.group=NULL,arraytype="450K")

请注意,这里有两个特殊参数:runDMP 意思是Block.GUI()函数是否为每个CpG计算DMP,if not,GUI()函数将只返回所有涉及的CpG的注释。如果要在运行Block.GUI()时计算DMP,则可以对compare.group赋值,即使我们的协变量中有多个表型,它还是可运行的。

image.png

就像DMP.GUI()和DMR.GUI()一样,P值cutoff和最小集群(类似于DMR.GUI()中的最小探针)控制您选择的显著块的数目。然后显著块的信息将在右侧列出。

image.png

第二个选项卡提供所有块中涉及的所有CpG的信息。对于甲基化块(methylation blocks),它实际上有三种单元集成在一起:块,簇(区域)和CpG。块由一些簇形成,并且簇由一些CpGs形成。因此,第二个选项卡实际上提供了从CpG到块的映射方案。用户可以使用此选项卡找到感兴趣的的基因和CpGs。然后使用DMP.GUI()来整合块(blocks)和DMPs之间的研究。

image.png
第三个选项卡提供每个块的图。请注意,每个Block可能包含许多簇。因此,计算和绘图就会很慢。请确保您有一台好的电脑或服务器来运行这个功能。请注意,这个图与DMR.GUI()和DMP.GUI()略有不同,X轴实际上是每个集群的REAL MAPINFO,因为我还需要在线上方映射CpG岛(那些深红色的点)。

6.10基因集富集分析

基因集富集分析是很多生物信息学研究工作的重要一步。在之前的步骤中,您可能已经获得了一些重要的DMP或DMR,因此您可能想知道这些重要的DMP或DMR中涉及的基因是否富集到特定的biological terms 或pathways。为了实现这一分析,您可以使用champ.GSEA()来进行GSEA分析。

champ.GSEA()将自动提取基因myDMPmyDMR(无论您用什么方法生成DMR)。此外,如果您有自己的显著的CpG列表或基因列表,而且来自于THE SAME甲基化阵列(450K或EPIC)通过其他未在ChAMP中的方法计算。您也可以将其格式化成list并输入函数中以进行GSEA(请将您自己的CpG列表中的每个元素或基因列表分配一个名称,否则可能会触发错误)。champ.GSEA()将自动提取基因信息,将CpG信息转换成基因信息,然后在每个列表上进行GSEA。在CpG与基因之间的映射过程中,如果将多个CpG映射到一个基因,那么该基因只会被计数一次。

现在有两种途径实现GSEA。在以前的版本中,ChAMP使用从MSIMDB下载的通路信息。然后,用Fisher精确检验计算每个通路的富集状态。基因富集分析后,champ.GSEA()函数将自动返回P值小于临界值(cutoff)的通路。

然而,正如Geeleher所指出的,由于不同基因中含有不同数量的CpG,存在的两种情况是一个基因内含50个CpG,但只有一个显示显著甲基化,一个基因内含2个显著CpG甲基化,这两种情况不能同等对待。解决方案是使用基因包含的CpG个数来校正显著的基因。如在gmeth函数(missMethyl包里的函数)中实现的。在gometh函数中,它使用每个基因的CpG数量替换长度作为偏置数据,以纠正此问题。gometh的思路适用于与GSEA相关的基因的CpG数量的曲线(The idea of gometh is fitting a curve for numbers of CpGs across genes related with GSEA),然后使用概率加权函数来校正GO的p值。

注意:如果用户要纠正基因数与CpG数之间的偏差,则可以将champ.GSEA()函数中的方法参数设置为“goseq”,使用goseq方法进行GSEA,或者用户将其设置为“fisher“做正常的基因集富集分析。

用户可以使用以下命令来执行GSEA:

myGSEA <- champ.GSEA(beta=myNorm,DMP=myDMP[[1]], DMR=myDMR, arraytype="450K",adjPval=0.05, method="fisher")
# myDMP and myDMR could (not must) be used directly.

如果您不想在DMP或DMR上执行GSEA,则可以将其设置为NULL。

在上面的代码中,我们还使用了包含在我们的软件包中的450K示例数组来展示champ.GSEA()的效果。计算后,对于每个列表(DMP,DMR或其他列表),将返回一个GSEA结果。用户可以检查结果(如下):

head(myGSEA$DMP)
## NULL
# Above is the GSEA result for differential methylation probes.
head(myGSEA$DMR)
##                                                         Gene_List nList
## BENPORATH_EED_TARGETS                       BENPORATH_EED_TARGETS  1062
## BENPORATH_ES_WITH_H3K27ME3             BENPORATH_ES_WITH_H3K27ME3  1118
## BENPORATH_SUZ12_TARGETS                   BENPORATH_SUZ12_TARGETS  1038
## HOX_GENES                                               HOX_GENES    65
## BENPORATH_PRC2_TARGETS                     BENPORATH_PRC2_TARGETS   652
## MEISSNER_BRAIN_HCP_WITH_H3K27ME3 MEISSNER_BRAIN_HCP_WITH_H3K27ME3   269
##                                  nRep      fRep nOVLAP        OR
## BENPORATH_EED_TARGETS             966 0.9096045    118  4.816181
## BENPORATH_ES_WITH_H3K27ME3       1031 0.9221825    118  4.457708
## BENPORATH_SUZ12_TARGETS           909 0.8757225    102  4.257157
## HOX_GENES                          60 0.9230769     29 28.960565
## BENPORATH_PRC2_TARGETS            594 0.9110429     79  5.037991
## MEISSNER_BRAIN_HCP_WITH_H3K27ME3  259 0.9628253     48  7.200159
##                                       P.value      adjPval
## BENPORATH_EED_TARGETS            1.002738e-36 8.364836e-33
## BENPORATH_ES_WITH_H3K27ME3       5.684943e-34 2.371190e-30
## BENPORATH_SUZ12_TARGETS          1.482686e-28 4.122856e-25
## HOX_GENES                        1.863102e-27 3.885499e-24
## BENPORATH_PRC2_TARGETS           6.092783e-27 1.016520e-23
## MEISSNER_BRAIN_HCP_WITH_H3K27ME3 7.157888e-23 9.951850e-20
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        Genes
## BENPORATH_EED_TARGETS            HOXC5 HOXD4 TP73 FLJ45983 DKK1 HOXD9 SIX6 EVX1 HIST1H3I NEUROD2 KCNQ1 TRPA1 SOX17 CHST8 HOXD3 CCDC140 DLX5 DCC ADCYAP1 GHSR PAX3 HOXA2 MKX BTG4 EMX2 SIM1 CRTAC1 ISL2 HIST1H4F FEZF2 DLK1 HOXC6 HSPA1L HOXB7 GBX2 HOXB8 GSC GATA3 MSX1 HOXC9 CYP26C1 DUOXA1 ALX4 HOXA3 TAP1 EYA4 DSC3 HOXA9 SFRP1 HOXD8 FOXF1 QRFPR VSX1 HCG9 HOXC11 TAC1 HTR1E OTP ALX1 HOXA10 CALCA HOXD12 RAB32 CIDEA HIST1H4L PDGFRA HOXA7 FOXG1 PTGDR HIST1H3C BARHL2 ZIC4 APBB1IP OLIG3 HSPA1A LBX1 PITX1 IGFBP3 EN1 CBLN1 ZIC1 PAX6 DUOX1 PRDM14 HOXB1 HOXD1 ADARB2 CMTM2 MAGI2 HOXA6 HIC1 HTR1A RASGRF1 GATA4 HOXD13 TBX5 LHX8 UCN HLX DLX1 NKX2-3 FLJ32063 PAX9 HOXC13 HOXA4 TLX1 SIM2 HPSE2 PDX1 LOC153684 TMEM26 HIST1H3E DRD5 HS3ST3B1 VAX1 HOXC4 MAB21L1 PRLHR
## BENPORATH_ES_WITH_H3K27ME3                  GBX2 OTP KLHL1 VAX1 PRLHR CALCA LRAT HSPA1A DUOXA1 HOXD8 OLIG3 LHX8 PCDHGC4 EVX1 DSC3 SOX17 HTR1A HLA-F SFRP1 PITX1 ZIC4 PAX3 HLX SIM2 ESR1 EXOC3L2 RFX4 ISL2 BARHL2 DUOX1 HOXA3 SIM1 CRTAC1 KLHL14 GATA4 HOXC11 HCG9 CBLN1 NKX2-3 FOXG1 UCN LOC153684 TAP1 FEZF2 DLK1 PAX6 RASGRF1 EYA4 HPSE2 ADARB2 HOXD1 HOXB8 ADCYAP1 GABRG3 TAC1 EGFLAM OSR2 CCNA1 EN1 HOXD13 HOXB7 ZIC1 KCNQ1 GHSR GATA3 PAX9 PTGDR ALX4 HOXA2 HOXD4 PDGFRA PTPRN2 TBX5 CMTM2 PDX1 HOXD12 HOXD9 HSPA1L MSX1 IGFBP3 NEUROD2 GPR150 DCC DLX1 CCDC140 MAB21L1 HOXB1 CHST8 CYP26C1 HOXD3 TLX1 FLJ45983 PRRT1 TP73 CIDEA FOXF1 DRD5 SIX6 HOXA6 HOXA9 HOXC4 THBS2 HOXA10 HOXC5 HIC1 DKK1 HOXA4 MKX DLX5 VSX1 HOXA7 HOXC6 FLJ32063 HS3ST3B1 ZFPM1 SLC6A2 LBX1 GSC
## BENPORATH_SUZ12_TARGETS                                                                                                                 HOXC9 CIDEA RASGRF1 CRTAC1 SFRP1 FLJ45983 PRLHR FOXF1 CBLN1 TMEM26 TLX1 DUOXA1 ADCYAP1 EMX2 GATA5 HOXC6 MAB21L1 GHSR GSC FOXG1 ALX4 MSX1 ISL2 HOXD12 ADARB2 DLX1 ZIC1 CCDC140 LOC153684 SPOCK2 QRFPR HOXD13 PAX6 HOXD1 HPSE2 HOXC11 DKK1 VSX1 LHX8 GPR158 PAX9 VAX1 HOXC4 NEUROD2 DCC PAX3 MKX CYP26C1 PITX1 UCN DRD5 HOXC5 GATA4 SLC6A2 PDX1 HOXC13 GPR83 CHST8 HOXB1 PTGDR TOX2 TP73 GBX2 TBX5 CRABP1 EN1 APBB1IP GATA3 FLJ32063 SIX6 ZAR1 TRPA1 HOXD4 ZIC4 CCNA1 HTR1A OSR2 HOXB8 PDGFRA HOXD9 CALCA ACCN1 SOX17 BARHL2 EGFLAM HOXB7 FEZF2 DSC3 ALDH1L1 HS3ST3B1 LBX1 GIPC2 HLX ALX1 NKX2-3 DUOX1 SIM2 UPB1 HOXD8 HOXD3 OTP CMTM2
## HOX_GENES                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 HOXA6 HOXC5 HOXD10 HOXB5 LBX1 HOXC6 HOXA5 HOXD9 HOXA2 HOXC11 HOXD13 HOXA4 HOXA11 HOXC13 PITX1 HOXD1 HOXD12 HOXA3 HOXA10 HOXD4 HOXB7 HOXC4 HOXB8 HOXA7 HOXD3 HOXB1 HOXB9 TLX1 HOXA9
## BENPORATH_PRC2_TARGETS                                                                                                                                                                                                                                                                 TP73 CCDC140 HOXD13 SOX17 PDX1 TBX5 OTP HOXC5 HOXD3 DKK1 PAX6 PDGFRA DLX1 HOXD12 ZIC1 CHST8 CRTAC1 CIDEA DSC3 HOXC11 ISL2 NEUROD2 VAX1 EN1 GBX2 ADARB2 MKX PAX9 PRLHR DCC TLX1 RASGRF1 GHSR UCN ADCYAP1 HOXD1 GATA4 SIM2 HTR1A HLX LOC153684 GSC LHX8 SIX6 FEZF2 HOXB7 HOXD8 PTGDR PAX3 DUOX1 DRD5 LBX1 NKX2-3 CALCA HPSE2 CMTM2 GATA3 ZIC4 FOXF1 HOXC4 HS3ST3B1 HOXB8 FLJ32063 MSX1 PITX1 FLJ45983 CBLN1 HOXC6 DUOXA1 BARHL2 MAB21L1 CYP26C1 HOXB1 SFRP1 HOXD4 VSX1 HOXD9 ALX4 FOXG1
## MEISSNER_BRAIN_HCP_WITH_H3K27ME3                                                                                                                                                                                                                                                                                                                                                                                                                                                       HOXB5 LOXL1 PRDM6 LHX8 HOXD1 GATA5 ALX1 EVX1 HOXC4 TNXB SIM1 OTP TLX1 PRDM14 OLIG3 HOXC13 DBX1 HOXD9 HOXB9 SLC6A2 SIM2 GSC HOXC11 TP73 HOXC5 KRT7 HOXA11 HOXD8 ALX4 PRSS50 HOXA4 HOXB7 GATA4 HOXD12 HOXD10 LBX1 GATA3 HOXA10 GBX2 BARHL2 PITX1 HOXA5 SIX6 LHX1 DLK1 HOXD13 HOXB4 PDX1
# Above is the GSEA result for differential methylation regions.
# Too many information may be printed, so we are not going to show the result here.

6.11差异甲基化互作热点

(Differential Methylated Interaction Hotspots)

ChAMP中的另一个功能是champ.EpiMod()。此函数使用了FEM包在用户特定的功能基因网络内推断差异甲基化的基因模块。该网络可以是例如蛋白质 - 蛋白质互作网络。因此,可以将champ.EpiMod()函数视为使用基因(通常为PPI网络)之间的关系网络的功能监督算法来识别其中大量基因与感兴趣的表型相关联的子网络( POI)。EpiMod算法可以以两种不同的模式运行:在探针级别,在这种情况下,最差异的甲基化探针被分配给每个基因,或者在基因水平上,在这种情况下,使用优化的DNA将DNAm值分配给每个基因。FEM包被开发用来推断差异甲基化基因模块(which are also deregulated at the gene expression level),然而这里我们仅提供EpiMod version,其仅用于推导差异甲基化模块。更多高级的用户可参考FEM软件包的帮助文档获取更多信息。

用户可以使用以下代码来做champ.EpiMod()

myEpiMod <- champ.EpiMod(beta=myNorm,pheno=myLoad$pd$Sample_Group)

所有重要的模块将被返回。而PDF verdion图将在resultDir(函数中的参数)返回。如下所示:

image.jpeg

上面我们展示了基于我们的450K测试数据的champ.EpiMod()检测到的四个模块。但在实践中,每个模块将被单独保存在resultsDir目录中一个pdf文件中。

6.12拷贝数变异

对于Copy Number Variation分析,我们提供champ.CNA()来实现这个工作。该函数使用HumanMethylation450或HumanMethylationEPIC数据来识别拷贝数变异。该函数利用每个探头的强度值来计数拷贝数,并确定是否存在拷贝数变异。拷贝数使用CopyNumber包确定。

基本上有两种方法可进行CNA分析:您可以比较病例样本和对照样本之间的拷贝数; 或者您可以将每个样品的拷贝数与所有样品的平均拷贝数状态进行比较。对于第一种方法,您可以在您的表型参数中指定哪些样本组,或者 ChAMP 已经包括了一些血液对照样品本身,你可以用它们作为对照,然后计算拷贝数变异,所有你需要做的就是分配参数control=TRUEcontrolGroup="champCtls"。对于其他方法,您可以分配control=FALSE

将产生两种图,每个样本的分析(sampleCNA=TRUE参数)和每个组(groupFreqPlots=TRUE参数)的分析。像之前的champ.QC()函数一样,有两个参数可以在绘图上分配。如果champ.CNA()将用R会话绘制,则需控制Rplot参数,如果champ.CNA()会将PDF文件保存到resultsDir,则需控制PDFplot参数。默认情况下,只有PDFplot为TRUE。

注意新版与旧版ChAMP不同。如果你想要对强度数据矫正批次效应,新版本则不会自动运行批次矫正,请使用champ.runCombat(),强度数据的用法与beta或M矩阵是完全相同的,唯一的区别是需要使用myLoad$intensity和设置来替换beta值,并蛇者参数logitTrans=FALSE

Feber (参考文献第10篇)比较了使用该方法从Illumina CytoSNP阵列和Affymetrix SNP 6.0阵列拷贝数据的结果,发现使用该方法对450k数据进行扩增和删失区域的识别是有效的。

用户可以使用以下代码来计算拷贝数差异:

myCNA <- champ.CNA(intensity=myLoad$intensity,pheno=myLoad$pd$Sample_Group)

champ.CNA()函数的主要用途是生成样本和组的表型图(The main usage of champ.CNA() function is to generate plot for sample and groups in pheno.)。如下所示:


image.png

6.13细胞类型异质性

DNA甲基化谱通常产生于有许多细胞类型混合的样品中。鉴于DNAm具有高度细胞型特异性,因此通常需要推断DMP / DMR哪些不是由细胞类型成分的潜在变化所导致的。目前已经提出了许多方法来处理细胞型异质性问题,如RefbaseEWAS。

RefbaseEWAS method uses a reference-database of DNAm profiles of the major cell-types thought to be present in the tissue of interest, and uses this reference to infer sample-specific cell-type proportions using a form of constrained multivariate regression, known as constrained projection or quadratic programming. These cell-type proportions can subsequently be used as covariates to infer phenotype-associated changes which are not driven by changes in cell-type composition. We included this refbase method into our ChAMP for users who want to detect cell proportion and remove cell type influence on their WHOLE BLOOD data.

在ChAMP中,我们包括一个全血参考数据库,一个用于27K,另一个用于450K。在使用champ.refbase()函数之后,将返回每个样本中的细胞类型异质性校正β基因和细胞类型特异性比例。请记住champ.refbase()只能用于血液样本数据集。

myRefBase <- champ.refbase(beta=myNorm,arraytype="450K")
# Our test data set is not whole blood. So it should not be run here.

请注意,我们的示例数据450K阵列不是全血数据,所以champ.refbase()无法在此数据集上使用。(以上代码仅在这里作为演示)。将来,我们将为champ.refbase函数添加更多参考数据。

7总结

我们简要介绍了新的ChAMP包中包含的各种流程和函数。在未来,我们希望进一步改进它,实现更多的功能。如果您有任何问题,错误报告或改进ChAMP的建议,请发送电子邮件至champ450k@gmail.com。Feeback是改进软件包的关键,该软件包试图从许多其他工具中无缝集成功能和并保证灵活性。

8引用ChAMP

ChAMP是一个包含许多其他来源的算法和函数的包装器。ChAMP中使用的大多数函数都是直接从其他软件包中引入,或者是从特定的已发布算法中获取的。因此,下面我们提供指导文件,在使用这些函数时,应该在ChAMP旁边引用哪些论文:

如果您使用ChAMP包,请引用:

Morris, T. J., Butcher, L. M., Feber, A., Teschendorff, A. E., Chakravarthy, A. R., Wojdacz, T. K., and Beck, S. (2014). Champ: 450k chip analysis methylation pipeline pg - 428-30. Bioinformatics, 30(3), 428-30.

The loading part of ChAMP (champ.load) uses minfi package, so please cite:

Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD and Irizarry RA (2014). Minfi: A flexible and comprehensive Bioconductor package for the analysis of Infinium DNA Methylation microarrays. Bioinformatics, 30(10), pp. 1363-1369. doi: 10.1093/bioinformatics/btu049.

Jean-Philippe Fortin, Timothy Triche, Kasper Hansen. Preprocessing, normalization and integration of the Illumina HumanMethylationEPIC array. bioRxiv 065490; doi: https://doi.org/10.1101/065490

Filtering out probes with SNPs was based following recommendations presented in this paper:

Zhou W, Laird PW and Shen H: Comprehensive characterization, annotation and innovative use of Infinium DNA Methylation BeadChip probes. Nucleic Acids Research 2016

FunctionalNormliazation and SWAN normalization methods are based on:

Jean-Philippe Fortin, Timothy Triche, Kasper Hansen. Preprocessing, normalization and integration of the Illumina HumanMethylationEPIC array. bioRxiv 065490; doi: https://doi.org/10.1101/065490

Maksimovic J et a. SWAN: Subset-quantile within array normalization for illumina infinium humanmethylation450 beadchips. Genome Biol. 2012;13(6):R44.

The default type-2 probe correction method in ChAMP is BMIQ, based on:

Teschendorff AE et al . A beta-mixture quantile normalization method for correcting probe design bias in illumina infinium 450 k dna methylation data. Bioinformatics. 2013;29(2):189-196.

Another type-2 probe correction method is PBC, which was presented in:

Dedeurwaerder S et a. Evaluation of the infinium methylation 450K technology. Epigenomics. 2011;3(6):771-784.

champ.SVD() function is based on:

Teschendorff AE et a. An epigenetic signature in peripheral blood predicts active ovarian cancer. PLoS One. 2009;4(12):e8274.

ChAMP uses the batch correction method Combat as implanted in the SVA package. Please cite:

Johnson WE et a. Adjusting batch effects in microarray expression data using empirical bayes methods. Biostatistics. 2007;8(1):118-127.

Leek JT, Johnson WE, Parker HS, Fertig EJ, Jaffe AE, Storey JD, Zhang Y and Torres LC (2017). sva: Surrogate Variable Analysis. R package version 3.24.4.

Above are all functions related to Data Normalization. Below are functions and methods for downstream analyses.

If using limma to find DMPs, please cite:

Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3. Epub 2004 Feb 12. PubMed PMID: 16646809.

Wettenhall JM, Smyth GK. limmaGUI: a graphical user interface for linear modeling of microarray data. Bioinformatics. 2004 Dec 12;20(18):3705-6. Epub 2004 Aug 5. PubMed PMID: 15297296.

For DMR detection, if you used Bumphunter method, please cite:

Jaffe AE et a. Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies. Int J Epidemiol. 2012;41(1):200-209.

If you used ProbeLasso function, please cite:

Butcher LM, Beck S. Probe lasso: A novel method to rope in differentially methylated regions with 450K dna methylation data. Methods. 2015;72:21-28.

Or if you used DMRcate function, please cite:

Peters TJ, Buckley MJ, Statham AL, et al. De novo identification of differentially methylated regions in the human genome. Epigenetics & Chromatin. 2015;8(1):1-16.

In this version of ChAMP, we provided a function to find Differential Methylation Blocks (DMB). The reference should be:

Hansen KD, Timp W, Bravo HC, et al. Increased methylation variation in epigenetic domains across cancer types. Nat Genet. 2011;43(8):768-775.

Timp W, Bravo HC, McDonald OG, et al. Large hypomethylated blocks as a universal defining epigenetic alteration in human solid tumors. Genome Med. 2014;6(8):61.

For GSEA method, if you used ‘gometh’ method, please cite below paper:

Paul Geeleher, Lori Hartnett, Laurance J. Egan, Aaron Golden, Raja Affendi Raja Ali, and Cathal Seoighe Gene-set analysis is severely biased when applied to genome-wide methylation data Bioinformatics 2013 : btt311v2-btt311.

Young MD, Wakefield MJ, Smyth GK and Oshlack A (2010). “Gene ontology analysis for RNA-seq: accounting for selection bias.” Genome Biology, 11, pp. R14.

Phipson, B, Maksimovic, J, Oshlack, A (2016). missMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics, 32, 2:286-8.

If you used the FEM package to find Differential Methylated Modules in gene networks, please cite:

Jiao Y, Widschwendter M, Teschendorff AE. A systems-level integrative framework for genome-wide dna methylation and gene expression data identifies differential gene expression modules under epigenetic control. Bioinformatics. 2014;30(16):2360-2366.

The CNA detection function was provided by Andy Feber of UCL, you may citate related paper below:

Feber A, Guilhamon P, Lechner M, et al. Using high-density dna methylation arrays to profile copy number alterations. Genome Biol. 2014;15(2):R30.

We included two methods for cell proportion correction. If you used the Refbase method in ChAMP, please cite:

Houseman EA, Accomando WP, Koestler DC, et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinformatics. 2012;13(1):1-16.

帮助文档原文链接:https://bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html#section-hydroxymethylation-analysis

扫描下方二维码关注生信客部落公众号:


生信客部落
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,884评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,347评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,435评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,509评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,611评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,837评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,987评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,730评论 0 267
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,194评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,525评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,664评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,334评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,944评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,764评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,997评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,389评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,554评论 2 349

推荐阅读更多精彩内容