ChAMP分析甲基化数据：标准流程

医学和生信笔记

发布于 2022-11-15 12:44:34

8070

发布于 2022-11-15 12:44:34

文章被收录于专栏：医学和生信笔记

之前学习了一下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

奇异值分解。

可用于找出比较重要的变化。

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

还有图形化界面：

DMP.GUI(DMP = myDMP[[1]],
        beta = myNorm,
        pheno = myLoad$pd$Sample_Type
        )

差异甲基化区域分析

Differentially Methylated Regions (DMRs)

myDMR <- champ.DMR(beta = myNorm,
                   pheno = myLoad$pd$Sample_Type,
                   arraytype="EPIC", # 注意选择！！
                   method = "Bumphunter"
                   )

提供图形化界面探索结果：

DMR.GUI(DMR=myDMR)

Differential Methylation Blocks

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

同样是提供图形化界面查看结果：

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，只要提供β矩阵和分组信息即可：

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的数据也是一模一样的流程！

上面所有的步骤，这个包还提供了一个一步法完成：

champ.process(directory = myDir)

不过还是分步运行更能获得自己想要的结果哦。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-09-08，如有侵权请联系 cloudcommunity@tencent.com 删除

文件存储

本文分享自医学和生信笔记微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

文件存储

登录后参与评论

0 条评论

热度