minfi 分析甲基化芯片数据 - 质量过滤篇

生信修炼手册

发布于 2020-05-09 17:22:39

1.9K00

代码可运行

文章被收录于专栏：生信修炼手册生信修炼手册

运行总次数：0

代码可运行

对于原始的芯片数据，在分析之前，我们首先要做的就是质量过滤，主要是探针水平的过滤，包含以下三个方面；

过滤掉可信度较低的探针；
过滤掉覆盖了snp位点的探针；
过滤掉位于性染色体上的探针；

过滤掉可信度较低的探针

对于探针信号的可信度，通过对应的pvalue 来标识。在minfi 中，通过detectP 这个函数计算探针对应的p值。用法示例

# import data 第一个参数为Sampel_Sheet.csv 文件所在的目录
targets <- read.metharray.sheet("./", pattern="HumanMethylation450_Demo_Sample_Sheet.csv")
rgSet <- read.metharray.exp(targets=targets)# probe pvalue
probeP <- detectionP(rgSet)

通常情况下，我们认为p > 0.01 的探针信号是不可信的，在实际操作中，只要在1个样本中其p值大于0.01. 我们就过滤掉这个探针。

keep <-  apply(probeP, 1 , function(t){all(t < 0.01)})
rgSet <- rgSet[keep,]

在实验操作不当时，会出现某个样本质量特别差的情况。通常认为所有探针pvalue 均值大于0.05的样本为质量差的样本，这样的样本也需要过滤掉。

keep <- apply(probeP, 2, mean) < 0.05
rgSet <- rgSet[,keep]

对于探针信号的可信度，直接读取原始数据就可以计算出来。但是对于探针的snp信息和染色体信息，就需要对应的探针注释文件。在bioconductor 中，IlluminaHumanMethylation450kanno.ilmn12.hg19 包提供了450K 芯片的探针注释信息；IlluminaHumanMethylationEPICanno.ilm10b4.hg19 包提供了850K 芯片的探针数据信息。minfi 在读取原始数据时，会根据探针数量自动判断芯片的类型，从而去对应的包中提取芯片注释信息。在进行snp和性染色探针过滤之前，我们首先要得到芯片的注释信息。

mSet <- preprocessRaw(rgSet)
Gset <- mapToGenome(mSet)
annotation <- getAnnotation(Gset)

过滤掉覆盖了snp位点的探针

如果探针覆盖到了snp位点，有可能影响其杂交情况。对于这部分探针，通常也会选择过滤掉。这时候就需要探针对应的snp注释信息。对于甲基化芯片而言，有3种SNP 注释信息

> snps <- getSnpInfo(Gset)
> head(snps)
DataFrame with 6 rows and 6 columns
             Probe_rs Probe_maf      CpG_rs   CpG_maf      SBE_rs   SBE_maf
          <character> <numeric> <character> <numeric> <character> <numeric>
cg13869341          NA        NA          NA        NA          NA        NA
cg14008030          NA        NA          NA        NA          NA        NA
cg12045430          NA        NA          NA        NA          NA        NA
cg20826792          NA        NA          NA        NA          NA        NA
cg00381604          NA        NA          NA        NA          NA        NA
cg20253340          NA        NA          NA        NA          NA        NA

CpG_rs

CpG_rs 指的是这个CpG位点同时也是一个snp位点。对于探针而言，CpG位点在基因组上的位置对应pos字段的信息，以第一行位置， rs71507462 在hg19的染色体位置就是762592。

> head(annotation[!is.na(annotation$CpG_rs), c(1,2,3,14,15)])
DataFrame with 6 rows and 5 columns
                  chr       pos      strand      CpG_rs   CpG_maf
          <character> <integer> <character> <character> <numeric>
cg06402284        chr1    762592           -  rs71507462  0.274931
cg03213877        chr1    763116           - rs191770232  0.020600
cg20788133        chr1    765028           -   rs2519025  0.092322
cg09139287        chr1    787398           -   rs2905055  0.336368
cg00168193        chr1    790667           + rs149504348  0.038900
cg01062849        chr1    838731           +  rs28477624  0.040256

SBE_rs

SBE_r指的是CpG位点下游的第一个碱基是snp位点，rs2905055 在hg19的位置是787399

> head(annotation[!is.na(annotation$SBE_rs), c(1,2,3,16,17)])
DataFrame with 6 rows and 5 columns
                  chr       pos      strand      SBE_rs   SBE_maf
          <character> <integer> <character> <character> <numeric>
cg09139287        chr1    787398           -   rs2905055  0.336368
cg23112672        chr1    876199           -  rs71628924  0.012351
cg10644916        chr1    899000           +   rs3813184  0.087352
cg10625579        chr1    942635           -  rs76233940  0.066406
cg13917577        chr1    957119           -   rs3121553  0.282192
cg11853970        chr1    986037           -  rs28758798  0.000456

Probe_rs

Probe_rs 指的是探针区域覆盖到的snp位点。探针长度为50bp左右，rs77418980在hg19的位置为 91536，在91550上游50bp范围内

> head(annotation[!is.na(annotation$Probe_rs), c(1,2,3,12,13)])
DataFrame with 6 rows and 5 columns
                  chr       pos      strand    Probe_rs Probe_maf
          <character> <integer> <character> <character> <numeric>
cg03130891        chr1     91550           -  rs77418980  0.305556
cg24335620        chr1    135252           - rs147502335  0.012800
cg24669183        chr1    534242           -   rs6680725  0.108100
cg14057946        chr1    713985           -  rs74512038  0.076823
cg11422233        chr1    714002           -  rs74512038  0.076823
cg16047670        chr1    714012           + rs114983708  0.115264

使用dropLociWithSnps 函数可以根据snp类型和频率进行顾虑，从探针的检测原理来说，CpG_rs和SBE_rs是对探针杂交影响最大的snp类型，所以常用这两种snp进行过滤。用法示例

GRset <- dropLociWithSnps(Gset, snps=c("Probe","SBE","CpG"), maf=0)

过滤掉位于性染色体上的探针

当样本中混合了两种性别时，无法正确的衡量性染色体上的甲基化水平的差异，所以需要去掉性染色体上的探针

sex_probe <- rownames(annotation)[annotation$chr %in% c("chrX", "chrY")]
keep <-  !(featureNames(GRset) %in% sex_probe)
GRset <- GRset[keep, ]

总结：

探针信号的可信度通过p值进行衡量，p值越小可信度越高，通常过滤掉p>0.01的探针。
由于snp影响探针的杂交情况，所以需要去除覆盖了snp位点的探针。探针上的SNP分为Probe_rs, CpG_rs, SBE_rs 3种类型，其中CpG_rs, SBE_rs 对甲基化水平的检测影响最大，在过滤snp时还可以结合snp频率进行过滤。
当分组比较的样本混合了不同性别，需要去除位于性染色体上的探针。

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

原始发表：2018-03-22，如有侵权请联系 cloudcommunity@tencent.com 删除

硬件开发

本文分享自生信修炼手册微信公众号，前往查看

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

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

硬件开发

登录后参与评论

0 条评论

热度