我是如何发现850K甲基化芯片和EPIC的区别

生信技能树

发布于 2020-07-31 14:17:27

2.2K1

有粉丝求助，她做850K甲基化芯片数据处理的时候，使用champ流程，然后报错非常诡异，是Failed CpG Fraction，各种查资料都无法解决，我让她复制粘贴报错的关键信息，如下：

[ Section 3: Use Annotation Start ]

  Reading 850K Annotation >>

  Fetching NEGATIVE ControlProbe.
    Totally, there are 613 control probes in Annotation.
    Your data set contains 556 control probes.

  Generating Meth and UnMeth Matrix
    Extracting Meth Matrix...
      Totally there are 485512 Meth probes in 850K Annotation.
      Your data set contains 485512 Meth probes.
    Extracting UnMeth Matrix...
      Totally there are 485512 UnMeth probes in 850K Annotation.
      Your data set contains 485512 UnMeth probes.

  Generating beta Matrix
  Generating M Matrix
  Generating intensity Matrix
  Calculating Detect P value
  Counting Beads
[ Section 3: Use Annotation Done ]
---
中间省略
---
[ Section 2: Filtering Start >>

  Filtering Detect P value Start
    The fraction of failed positions per sample
    You may need to delete samples with high proportion of failed probes:

         Failed CpG Fraction.
sample1                 NaN
sample2                NaN 
---后面省略一些样本
Error in if (any(numfail >= SampleCutoff)) { : 
  missing value where TRUE/FALSE needed

也帮忙去各种检索，但确实没有好的解决方案，就让她发过来2个G的原始数据和代码，认真检查了好久，看起来就是我的教程的代码，一模一样啊！

myLoad <- champ.load("raw/",arraytype="850K")

而且我看了她关于"raw/"文件夹下的idat文件，以及制作好的'raw/sample_sheet.csv'文件，都是合格的。没办法，我只好看champ.load函数的帮助文档了：

champ.load(directory = getwd(),
           method="ChAMP",
           methValue="B",
           autoimpute=TRUE,
           filterDetP=TRUE,
           ProbeCutoff=0,
           SampleCutoff=0.1,
           detPcut=0.01,
           filterBeads=TRUE,
           beadCutoff=0.05,
           filterNoCG=TRUE,
           filterSNPs=TRUE,
           population=NULL,
           filterMultiHit=TRUE,
           filterXY=TRUE,
           force=FALSE,
           arraytype="450K")

刚开始一直看不出问题所在，但是最后注意到了：

arraytype 这个参数的选择是：
Choose microarray type is "450K" or "EPIC".(default = "450K")

也就是说，没有850K这个选项，有意思，于是我修改了代码，如下：

#myLoad <- champ.load("raw/",arraytype="850K")
myLoad <- champ.load("raw/",arraytype="EPIC")

确实解决了这个报错，成功运行champ流程，载入idat文件后的日志如下：

 Filtering probes with a detection p-value above 0.01.
    Removing 3813 probes.
    If a large number of probes have been removed, ChAMP suggests you to identify potentially bad samples

  Filtering BeadCount Start
    Filtering probes with a beadcount <3 in at least 5% of samples.
    Removing 22027 probes

  Filtering NoCG Start
    Only Keep CpGs, removing 2889 probes from the analysis.

  Filtering SNPs Start
    Using general EPIC SNP list for filtering.
    Filtering probes with SNPs as identified in Zhou's Nucleic Acids Research Paper 2016.
    Removing 95451 probes from the analysis.

  Filtering MultiHit Start
    Filtering probes that align to multiple locations as identified in Nordlund et al
    Removing 11 probes from the analysis.

  Filtering XY Start
    Filtering probes located on X,Y chromosome, removing 16655 probes from the analysis.

  Updating PD file

  Fixing Outliers Start
    Replacing all value smaller/equal to 0 with smallest positive value.
    Replacing all value greater/equal to 1 with largest value below 1..
[ Section 2: Filtering Done ]

 All filterings are Done, now you have 725072 probes and 24 samples.

很有意思哦，850K甲基化芯片和EPIC的差异是？我明明是在各种教程及文档，看到850K甲基化芯片和EPIC是同一个芯片的不同表述而已：

Illumina公司提供了一个更强大的甲基化分析平台：Illumina InfiniumMethylationEPIC BeadChip （DNA甲基化850K芯片），不但包含了原450K芯片90%以上的位点，并额外增加了增强子区的350,000个位点，可以对正常样本和FFPE样本单个CpG位点进行定量甲基化检测，该芯片是目前最适合甲基化图谱分析研究的全基因组DNA甲基化芯片。
850K芯片覆盖了全基因组853,307个CpG位点，全面覆盖CpG岛、启动子、编码区及增强子。覆盖CpG岛、RefSeq基因、ENCODE开放染色质、ENCODE转录因子结合位点、FANTOM5增强子区域。

这就是很神奇了，但我又不是公司客服，懒得去探索了。