之前学习了一下ChAMP包读取IDAT文件，真的是太贴心！
上次主要演示了ChAMP包需要的样本信息csv文件的制作以及IDAT数据读取过程。
ChAMP分析甲基化数据：样本信息csv的制作和IDAT读取
今天继续走完后面的流程，很多日志文件我没放上来。
数据质控
读取数据之后需要进行一些质控。
直接一个函数搞定：champ.QC()。
champ.QC(beta = myLoad$beta,
pheno = myLoad$pd$Sample_Type # 注意列名要选对
)
会生成3张图，放在CHAMP_QCimages这个文件夹下。
- MDS plot：根据前1000个变化最大的位点看样品相似性。
- densityPlot：每个样品的beta分布曲线，比较离群的可能是质量比较差的样本。
- 聚类图
标准化
使用champ.norm()函数实现，提供4种方法：
- BMIQ,
- SWAN,
- PBC,
- FunctionalNormliazation
FunctionalNormliazation需要rgSet对象，SWAN需要rgSet和mset，PBC和BMIQ只需要beta 矩阵，FunctionalNormliazation和SWAN需要在读取数据时使用method = "minfi"。
myNorm <- champ.norm(beta = myLoad$beta,
arraytype = "EPIC",
cores = 8
)
SVD
奇异值分解。
可用于找出比较重要的变化。
代码语言：javascript复制champ.SVD(beta = myNorm |> as.data.frame(), # 这里需要注意
pd=myLoad$pd)
这个图中只有Sample_Type这个变量比较重要，如果你在csv文件中提供更多样本信息，这个图会看起来像热图。
save(myNorm,myLoad,file = "EPIC.rdata")
矫正批次效应
借助了sva包的Combat函数实现。
不是所有的都需要，根据自己的实际情况来。
注意使用时需要指定列名。
如果使用M值，需要指定logitTrans = T。
# 我们这个数据没有批次效应，就不运行这一步了
myCombat <- champ.runCombat(beta = myNorm,
pd = myLoad$pd,
batchname = c("Slide"), # pd文件中哪一列是批次信息
variablename = c("Sample_Type") # 指定分组列名，默认是Sample_Group
)
甲基化差异分析
甲基化有3个层次的差异分析：DMP，DMR，DMB：
- DMP代表找出Differential Methylation Probe（差异化CpG位点）
- DMR代表找出Differential Methylation Region（差异化CpG区域）
- Block代表Differential Methylation Block（更大范围的差异化region区域）
甲基化位点差异分析
DMP, Differentially Methylated Probes.
借助了limma包进行差异分析。支持多个分组的比较，如果分组大于2组，会自动进行两两比较。也支持数值型变量，比如年龄这种。
myDMP <- champ.DMP(beta = myNorm,
pheno = myLoad$pd$Sample_Type,
arraytype = "EPIC" # 别忘了改这里！！
)
结果就是myDMP[[1]],我们查看前6个，信息很多，一共20列：
head(myDMP[[1]])
logFC AveExpr t P.Value adj.P.Val B normal_AVG cancer_AVG deltaBeta CHR MAPINFO
cg16601494 -0.6402505 0.3910975 -28.66283 1.344757e-19 9.927870e-14 33.28137 0.07097231 0.7112228 0.6402505 1 1475737
cg09296001 -0.5763521 0.3706256 -27.71457 2.870274e-19 1.059511e-13 32.67560 0.08244952 0.6588016 0.5763521 7 127672564
cg22697045 0.3923163 0.6356132 25.87853 1.339179e-18 2.882866e-13 31.41759 0.83177142 0.4394551 -0.3923163 11 47359223
cg17301223 -0.6012424 0.4088559 -25.70149 1.561968e-18 2.882866e-13 31.28995 0.10823474 0.7094771 0.6012424 8 145106438
cg03241244 -0.5798955 0.3608224 -25.15421 2.529505e-18 3.734890e-13 30.88793 0.07087468 0.6507702 0.5798955 12 122017052
cg26256223 -0.5925507 0.3537249 -23.99254 7.275558e-18 7.140446e-13 29.99558 0.05744958 0.6500003 0.5925507 8 145106582
Strand Type gene feature cgi feat.cgi UCSC_Islands_Name SNP_ID SNP_DISTANCE
cg16601494 R I C1orf70 5'UTR shore 5'UTR-shore chr1:1476093-1476669 rs561063770 23
cg09296001 R I SND1 Body island Body-island chr7:127671158-127672853
cg22697045 F I MYBPC3 Body shore Body-shore chr11:47359953-47360216 rs3729950 21
cg17301223 R I OPLAH Body island Body-island chr8:145103285-145108027 rs566592895;rs148937107 51;47
cg03241244 F I KDM2B Body island Body-island chr12:122016170-122017693 rs112471261;rs531357599 21;33
cg26256223 R II OPLAH Body island Body-island chr8:145103285-145108027 rs553503071 20
还有图形化界面：
代码语言：javascript复制DMP.GUI(DMP = myDMP[[1]],
beta = myNorm,
pheno = myLoad$pd$Sample_Type
)
差异甲基化区域分析
Differentially Methylated Regions (DMRs)
代码语言：javascript复制myDMR <- champ.DMR(beta = myNorm,
pheno = myLoad$pd$Sample_Type,
arraytype="EPIC", # 注意选择！！
method = "Bumphunter"
)
提供图形化界面探索结果：
代码语言：javascript复制DMR.GUI(DMR=myDMR)
Differential Methylation Blocks
代码语言：javascript复制myBlock <- champ.Block(beta=myNorm,
pheno=myLoad$pd$Sample_Type,
arraytype="EPIC"
)
head(myBlock$Block)
chr start end value area cluster indexStart indexEnd L clusterL p.value fwer p.valueArea fwerArea
Block_1 5 142602343 171118422 0.1663742 512.5989 74 225817 228897 3081 8726 0 0 0.000000e 00 0.000
Block_2 12 73523354 112200054 0.1249276 496.5871 168 67195 71169 3975 7196 0 0 0.000000e 00 0.000
Block_3 2 18435 25360029 0.1308144 423.4463 20 137105 140341 3237 10107 0 0 0.000000e 00 0.000
Block_4 7 22371030 43698429 0.1640385 418.1342 92 250282 252830 2549 4868 0 0 0.000000e 00 0.000
Block_5 11 118313315 134945796 0.1556909 405.2634 161 56593 59195 2603 4782 0 0 2.413494e-06 0.002
Block_6 2 48641899 85085187 0.1250336 401.1079 20 143636 146843 3208 10107 0 0 2.413494e-06 0.002
同样是提供图形化界面查看结果：
代码语言：javascript复制Block.GUI(Block=myBlock,
beta=myNorm,
pheno=myLoad$pd$Sample_Type,
runDMP=TRUE,
compare.group=NULL,
arraytype="EPIC")
save(myDMP,myDMR,myBlock,file = "./gse149282/gse149282_dmprb.rdata")
富集分析
提供GSEA分析的函数，这一步完全可以使用更加专业的clusterprofiler。
myGSEA <- champ.GSEA(beta = myNorm,
DMP = myDMP[[1]],
DMR = myDMR,
arraytype = "EPIC",
adjPval = 0.05,
method = "gometh"
)
head(myGSEA$DMP)
ONTOLOGY TERM N DE P.DE FDR
GO:0071944 CC cell periphery 5604 5332 1.072641e-74 2.430284e-70
GO:0005886 CC plasma membrane 5139 4885 1.605572e-65 1.818872e-61
GO:0032501 BP multicellular organismal process 6970 6478 2.138994e-40 1.615440e-36
GO:0007154 BP cell communication 6026 5607 3.068516e-35 1.691964e-31
GO:0031224 CC intrinsic component of membrane 5082 4739 3.733865e-35 1.691964e-31
GO:0016020 CC membrane 8634 7958 6.456365e-35 2.438031e-31
head(myGSEA$DMR)
ONTOLOGY TERM N DE P.DE FDR
GO:0003700 MF DNA-binding transcription factor activity 1336 195 1.108904e-25 2.512443e-21
GO:0000981 MF DNA-binding transcription factor activity, RNA polymerase II-specific 1290 189 6.300238e-25 7.137224e-21
GO:0000977 MF RNA polymerase II transcription regulatory region sequence-specific DNA binding 1334 190 3.119916e-23 2.356265e-19
GO:1990837 MF sequence-specific double-stranded DNA binding 1479 200 2.737843e-21 1.550783e-17
GO:0000976 MF transcription cis-regulatory region binding 1426 193 1.390539e-20 6.301089e-17
GO:0001067 MF transcription regulatory region nucleic acid binding 1428 193 1.801922e-20 6.804360e-17
还可以用贝叶斯方法，这种方法不需要DMP或者DMR，只要提供β矩阵和分组信息即可：
代码语言：javascript复制myebayGSEA <- champ.ebGSEA(beta = myNorm,
pheno = myLoad$pd$Sample_Type,
arraytype = "EPIC"
)
拷贝数变异分析
借助HumanMethylation450/HumanMethylationEPIC data实现这个功能，所以可能不能分析27K数据。
提供2种方法分析拷贝数变异。一种是分析两个状态间的拷贝数差异，比如说normal和cancer，另一种是分析单个样本的拷贝数和总体平均拷贝数。
使用第一种方法需要指定和谁比，如果是血液样本，这个包自带了一些血液样本可以比较，不需要指定controlGroup，但是如果不是血液样本，需要指定和谁比。
使用第二种方法只要选择control=FALSE即可。
# 非常耗时！
myCNA <- champ.CNA(intensity = myLoad$intensity,
pheno = myLoad$pd$Sample_Type,
arraytype = "EPIC",
controlGroup = "normal" # 指定和normal比
)
这个过程很慢，会在CHAMP_CNA文件夹中生成很多图。
标准化流程就是这么多，在ChAMP中都是一个函数搞定，基因注释等都是自动完成的，太方便了！
EPIC数据的甲基化分析在ChAMP中非常简单，就是这几步：
# 数据读取
myDir="./gse149282/GSE149282_RAW/"
myLoad <- champ.load(myDir, arraytype="EPIC")
# 数据预处理
champ.QC()
myNorm <- champ.norm(arraytype="EPIC")
champ.SVD()
myCombat <- champ.runCombat() # 可选
# 差异分析
myDMP <- champ.DMP(arraytype="EPIC")
myDMR <- champ.DMR(arraytype="EPIC")
myBlock <- champ.Block(arraytype="EPIC")
# 富集分析
myGSEA <- champ.GSEA(arraytype="EPIC")
# 拷贝数分析
myCNA <- champ.CNA(arraytype = "EPIC")
450K的数据也是一模一样的流程！
上面所有的步骤，这个包还提供了一个一步法完成：
代码语言：javascript复制champ.process(directory = myDir)
不过还是分步运行更能获得自己想要的结果哦。
