经过预处理之后的数据,就可以进行差异分析了。对于甲基化芯片而言,有两个方面的差异分析
DMP 差异甲基化探针
DMR 差异甲基化区域
在ChAMP
包中,champ.DMP
函数用于分析差异甲基化探针,champ.DMR
函数用于分析差异甲基化区域。本章我们先看下差异探针的分析
champ.DMP
函数的用法示例
myDMP <- champ.DMP()
在差异分析时,我们需要两个输入数据,一个就是探针的表达谱数据,beta matrix, 另外一个就是样本的分组信息。
在champ.DMP
函数中,默认myNorm
作为归一化之后的beta matrix,对于样本的分组信息,ChAMP
默认从Samplesheet.csv文件中读取,在数据导入成功后,myLoad$pd
代表的就是SampleSheet.csv文件的信息,所以myLoad$pd$Sample_Group
代表样本的分组信息。
在差异分析时,最关键的就是差异分组问题。在实验设计阶段,有很多类型的分组设计,比如最常见的case_vs_control, 两个group的分组;多个组织,比如3个组织,共3个group; 时间序列,比如药物处理后的几个时间点。不同的实验设计,在差异分析时,想要关注的差异点自然不同,在分析时也要采取不同的分析策略。
对于ChAMP
来说,上述的几种分组设计都是支持。
champ.DMP
计算过程分为以下3步:
1. 检测样本分组,确定分组比较
测试数据分成T和C两组,每组各4个样本、
在这一步,需要确定两个因素:
分组的类型
主要识别分组变量是离散型还是连续型,如果是字符串型character
, 就是离散型,如果是数值型numeric
, 就是连续型。不同的类型,采取的差异分析策略不同。
测试数据是字符型的两个group, 具体的输入信息如下分组的个数
确定group的个数,2个group 肯定是两者之间进行差异分析,但是当group 个数3个或以上时,就需要确定如何分组比较。 默认情况下两两之间都进行差异分析,如果你不需要这么多的差异结果,可以通过compare.group
参数指定,compare.group
参数的值是一个list, list 中的每个元素是一个长度为2的向量,指定了用于差异分析的两个group
[ Section 1: Check Input Pheno Start ]
You pheno is character type.
Your pheno information contains following groups. >>
C:4 samples.
T:4 samples.
[The power of statistics analysis on groups contain very few samples may not strong.]
pheno contains only 2 phenotypes
compare.group parameter is NULL, two pheno types will be added into Compare List.
C_to_T compare group : C, T
[ Section 1: Check Input Pheno Done ]
进行差异探针分析
通过调用limma
函数进行差异分析,默认通过BH
方法进行多重建设检验的校正,p.adjust < 0.05 的认为是差异探针
可以通过adjPVal
参数修改p.adjust的阈值,当然也可以修改adjust.method
参数的值,调整多重假设检验校正的算法,默认值为BH
, 可选值包括 “none”, “BH”, “BY”, “holm”。
[ Section 2: Find Differential Methylated CpGs Start ]
Start to Compare : C, T
Contrast Matrix
Contrasts
Levels pT-pC
pC -1
pT 1
You have found 4283 significant MVPs with a BH adjusted P-value below 0.05.
Calculate DMP for C and T done.
[ Section 2: Find Numeric Vector Related CpGs Done ]
添加差异探针的注释信息
之前的分析都是针对探针的beta matrix 进行的分析,找的差异探针之后,我们肯定希望知道这个探针对应的基因,染色体位置等注释信息。这一步实际就是在已有的差异结果的基础上,追加探针的注释信息。
[ Section 3: Match Annotation Start ]
[ Section 3: Match Annotation Done ]
结果展示
str(myDMP)
List of 1
$ C_to_T:’data.frame’: 4283 obs. of 23 variables:
..$ logFC : num [1:4283] 0.724 …
..$ AveExpr : num [1:4283] 0.398 …
..$ t : num [1:4283] 28.2 …
..$ P.Value : num [1:4283] 1.74e-08…
..$ adj.P.Val : num [1:4283] 0.00703…
..$ B : num [1:4283] 7.6 7.05 …
..$ C_AVG : num [1:4283] 0.0358 …
..$ T_AVG : num [1:4283] 0.759 …
..$ deltaBeta : num [1:4283] 0.724 …
..$ CHR : Factor w/ 25 levels “”,”1”,”10”,”11”,..:
..$ MAPINFO : int [1:4283] 141516291 …
..$ Strand : Factor w/ 3 levels “”,”F”,”R”: 3 2 …
..$ Type : Factor w/ 2 levels “I”,”II”: 2 2 …
..$ gene : Factor w/ 20622 levels “”,”A1BG” …
..$ feature : chr [1:4283] “Body” “Body” …
..$ cgi : Factor w/ 4 levels “island”,”opensea”,..: …
..$ feat.cgi : Factor w/ 28 levels “1stExon-island”,..: 13 …
..$ UCSC_CpG_Islands_Name: Factor w/ 27177 levels “”,”chr1:10003165-10003585”,..: 18278 …
..$ DHS : logi [1:4283] NA NA NA NA NA NA …
..$ Enhancer : logi [1:4283] TRUE TRUE NA NA NA NA …
..$ Phantom : Factor w/ 11912 levels “”,”high-CpG:100009860- …
..$ Probe_SNPs : Factor w/ 56004 levels “”,”rs10000615”,..:
..$ Probe_SNPs_10 : Factor w/ 35790 levels “”,”rs10000804” …
myDMP
就是最终的差异分析结果,是一个list
对象,list中的每个元素是两个group之间差异分析的结果。
测试数据只有两个分组,所以list 中只有一个元素。差异分析的结果是一个data.frame
对象,可以分成3个部分。
从logFC
到B
的部分是limma 差异输出结果, C_AVG
到deltaBeta
是每组表达量的均值,deltaBate
是两组均值的差,CHR
到Probe_SNPs_10
是探针的注释信息。