查看原文
其他

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

阿越就是我 医学和生信笔记 2023-06-15
关注公众号,发送R语言python,可获取资料

💡专注R语言在🩺生物医学中的使用


之前学习了一下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需要rgSetmset,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)

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




获取更多信息,欢迎加入🐧QQ交流群:613637742


医学和生信笔记,专注R语言在临床医学中的使用、R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。


往期推荐



dplyr中的行操作

dplyr中的across操作

dplyr强大的分组汇总

使用dplyr进行数据分析:入门篇

生信初学者基础知识资源推荐

初学者使用R语言读取、写出文件(csv/txt/excel/rdata等)的注意事项

R语言可视化聚类

R语言画好看的聚类

又是聚类分析可视化!



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存