美文网首页DNA甲基化R语言训练
甲基化样本和CpG位点QC的总流程（450k和850k）

甲基化样本和CpG位点QC的总流程（450k和850k）

作者: 橙子牛奶糖 | 来源:发表于2020-06-23 22:18 被阅读0次

    这篇应该是甲基化QC的最后一篇啦。

    感谢健明带入门。

    我前面已经写完两篇：

    QC1:甲基化数据QC：使用甲基化数据计算样本间的相关性

    QC2:甲基化数据QC: 使用甲基化数据推测SNP基因型（ewastools工具）

    下面补充一下对甲基化样本和CpG位点QC的总流程：

    1、导入、加载安装包

    if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GenomeInfoDbData")
BiocManager::install("IlluminaHumanMethylation450kmanifest")
BiocManager::install("IlluminaHumanMethylation450kanno.ilmn12.hg19")
BiocManager::install("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
BiocManager::install("IlluminaHumanMethylationEPICmanifest")
BiocManager::install("methylationArrayAnalysis")
BiocManager::install("limma")
BiocManager::install("minfi")
BiocManager::install("missMethyl")
BiocManager::install("minfiData")
BiocManager::install("Gviz")
BiocManager::install("DMRcate")
install.packages("knitr")
install.packages("RColorBrewer") 
library(knitr)
library(limma)
library(minfi)
library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(IlluminaHumanMethylation450kmanifest)
library(RColorBrewer)
library(missMethyl)
library(minfiData)
library(Gviz)
library(DMRcate)
library(stringr)
library("methylationArrayAnalysis")

    2、加载数据

    dataDirectory <- system.file("extdata", package = "methylationArrayAnalysis")
list.files(dataDirectory, recursive = TRUE)

    3、加载甲基化注释包

    ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
head(ann450k)
ann850k <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
head(ann850k)

    4、加载样本信息数据

    targets <- read.metharray.sheet(dataDirectory, pattern="SampleSheet.csv")
targets
    image

    5、读取甲基化的原始数据idat

    rgSet <- read.metharray.exp(targets=targets)

    6、将样本名添加到甲基化数据中

    targets$ID <- paste(targets$Sample_Group,targets$Sample_Name,sep=".")
sampleNames(rgSet) <- targets$ID
rgSet

    开始QC~

    7、甲基化cgp位点P值过滤

    原理：对每一个样本的每一个Cpg位点的总信号(M+U)和背景信号进行比较，可以得到P值。一般认为，越低的P值表示该位点越可靠，P值大于0.01的cpg位点，是质量比较差的位点；

    检测P值：

    detP <- detectionP(rgSet)
head(detP)

    结果如下图所示：
    红色框框为样本名，蓝色框框为为一个cpg位点的P值。


    image

    画每个样本cpg位点的平均P值

    pal <- brewer.pal(8,"Dark2")
par(mfrow=c(1,2))
barplot(colMeans(detP), col=pal[factor(targets$Sample_Group)], las=2, 
        cex.names=0.8, ylab="Mean detection p-values")
abline(h=0.05,col="red")
legend("topleft", legend=levels(factor(targets$Sample_Group)), fill=pal,
       bg="white")

    如下图所示：

    image

    可以看到，只有最后一个样本的cpg平均P值是超过0.05，也就是说，这个样本的质量是比较差的，后续应该被剔除掉。

    导出质量报道：
    qcReport(rgSet, sampNames=targetsID, sampGroups=targetsSample_Group,
    pdf="qcReport.pdf")

    7.1、剔除甲基化中高P值样本

    keep <- colMeans(detP) < 0.05
rgSet <- rgSet[,keep]
rgSet

    这里对P值设定的阈值是大于0.05。
    我们只保留cpg平均P值小于0.05的样本，对于P值大于0.05的样本（比如本例的birth.11）应被剔除。
    剔除以后，11个样本就只剩下10个样本：

    image

    7.2、剔除样本信息中高P值的样本

    targets <- targets[keep,]
targets[,1:5]

    7.3、剔除P值中高P值的样本

    detP <- detP[,keep]
dim(detP)

    8、甲基化标准化

    甲基化标准化是为了较少样本间的差异。

    有两种包可以进行甲基化的标准化工作，分别为preprocessFunnorm和preprocessQuantile。

    但这两个包的用途是不一样的。

    preprocessFunnorm包是针对甲基化数据来源于有明显分层的样本，比如癌症样本和正常样本，皮肤组织样本和大脑组织的样本。像这种明显有不同来源的样本建议用preprocessFunnorm包进行标准化。

    preprocessQuantile包则针对没有明显分层的样本，比如都是健康人群，都是来自血液样本这种情况。像这种单一来源的样本建议用preprocessQuantile包进行标准化。

    使用preprocessQuantile包进行标准化：

    mSetSq <- preprocessQuantile(rgSet)

    比较标准化前后的样本的beta值分布

    par(mfrow=c(1,2))
densityPlot(rgSet, sampGroups=targets$Sample_Group,main="Raw", legend=FALSE)
legend("top", legend = levels(factor(targets$Sample_Group)), 
       text.col=brewer.pal(8,"Dark2"))
densityPlot(getBeta(mSetSq), sampGroups=targets$Sample_Group,
            main="Normalized", legend=FALSE)
legend("top", legend = levels(factor(targets$Sample_Group)), 
       text.col=brewer.pal(8,"Dark2"))

    画出来的图如下所示，左边是未进行标准化的，右边是标准化以后的。

    可见，进行标准化后，样本间的差异会缩小。

    image

    9、查找标准化后数据可能存在的差异来源

    这一步是为EWAS做准备的，我们前面进行了标准化，但标准化的数据不代表就可以完全去除样本批次效应、细胞类型等差异。

    如果样本间存在批次效应等可能的混淆因素，在后续进行EWAS分析时极大可能会产生假阳性。

    因此我们需要通过主成分分析对标准化的数据进行可视化。确定可能的混淆因素，并在EWAS分析时进行校正。

    9.1、通过主成分1、2确认样本间的差异来源

    par(mfrow=c(1,2))
plotMDS(getM(mSetSq), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Group)])
legend("top", legend=levels(factor(targets$Sample_Group)), text.col=pal,
       bg="white", cex=0.7)

plotMDS(getM(mSetSq), top=1000, gene.selection="common",  
        col=pal[factor(targets$Sample_Source)])
legend("top", legend=levels(factor(targets$Sample_Source)), text.col=pal,
       bg="white", cex=0.7)

    如下图所示，可以很明显的看到这10个样本被分为三个聚类。

    image

    说明即便是前期进行了标准化处理后，样本间还是存在差异，比如样本act_naive.5和naive.1就很明显的在不同的聚类中。

    这种差异在进行EWAS分析时是我们不愿意看到的，因此后期进行EWAS分析时，应考虑将他们纳入协变量中。

    9.2、通过主成分1、2、3、4确认样本间的其他差异来源

    par(mfrow=c(1,3))
plotMDS(getM(mSetSq), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Group)], dim=c(1,3))
legend("top", legend=levels(factor(targets$Sample_Group)), text.col=pal, 
       cex=0.7, bg="white")

plotMDS(getM(mSetSq), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Group)], dim=c(2,3))
legend("topleft", legend=levels(factor(targets$Sample_Group)), text.col=pal,
       cex=0.7, bg="white")

plotMDS(getM(mSetSq), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Group)], dim=c(3,4))
legend("topright", legend=levels(factor(targets$Sample_Group)), text.col=pal,
       cex=0.7, bg="white")

    对主成分1,2,3,4画图后如下所示：

    image

    可以看到，不同样本按照不同颜色很好地被分层了，说明主成分3,4反应的是细胞类型的差异。

    同样的，细胞类型差异在进行EWAS分析时也是我们并不愿意看到的，因此他们在进行EWAS分析时应被一起纳入协变量中校正掉。

    10、探针过滤

    前面我们根据cpg的P值结果对样本进行了过滤。

    现在我们需要对探针进行过滤。

     detP <- detP[match(featureNames(mSetSq),rownames(detP)),] #匹配ID
 keep <- rowSums(detP < 0.01) == ncol(mSetSq) #对P值小于0.01的探针进行计数
 table(keep) #统计有多少个探针P值小于0.01
mSetSqFlt <- mSetSq[keep,] #保留P值在所有样本中均小于0.01的探针。
mSetSqFlt

    11、移除包含性染色体的探针

    keep <- !(featureNames(mSetSqFlt) %in% ann450k$Name[ann450k$chr %in% c("chrX","chrY")])
 table(keep)
 mSetSqFlt <- mSetSqFlt[keep,]

    12、移除SNP探针

    mSetSqFlt <- dropLociWithSnps(mSetSqFlt)
mSetSqFlt

    13.1、移除匹配在多个基因组上的探针（如果是450k）

    这一步是针对450k的数据，如果你的数据是850k，略过这一步，请看下面850k的工作。

    xReactiveProbes <- read.csv(file=paste(dataDirectory,
                                       "48639-non-specific-probes-Illumina450k.csv",
                                       sep="/"), stringsAsFactors=FALSE)
keep <- !(featureNames(mSetSqFlt) %in% xReactiveProbes$TargetID)
table(keep)
mSetSqFlt <- mSetSqFlt[keep,] 
mSetSqFlt

    13.2移除匹配在多个基因组上的探针（如果是850k）

    这一步是针对850k的数据，如果你的数据是450k，略过这一步，请看上面450k的工作。

    if (! ("devtools" %in% installed.packages()) install.packages("devtools")
devtools::install_github("markgene/maxprobes")
library(maxprobes) 
xloci <- maxprobes::xreactive_probes(array_type = "EPIC")
length(xloci)
mSetSqFlt <- maxprobes::dropXreactiveLoci(mSetSqFlt)

    14.1、重新评估是否已经消除样本间的差异（方法一：minfi包）

    par(mfrow=c(1,2))
plotMDS(getM(mSetSqFlt), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Group)], cex=0.8)
legend("right", legend=levels(factor(targets$Sample_Group)), text.col=pal,
       cex=0.65, bg="white")

plotMDS(getM(mSetSqFlt), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Source)])
legend("right", legend=levels(factor(targets$Sample_Source)), text.col=pal,
       cex=0.7, bg="white")
    image 
    par(mfrow=c(1,3))
plotMDS(getM(mSetSqFlt), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Source)], dim=c(1,3))
legend("right", legend=levels(factor(targets$Sample_Source)), text.col=pal,
       cex=0.7, bg="white")

plotMDS(getM(mSetSqFlt), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Source)], dim=c(2,3))
legend("topright", legend=levels(factor(targets$Sample_Source)), text.col=pal,
       cex=0.7, bg="white")

plotMDS(getM(mSetSqFlt), top=1000, gene.selection="common", 
        col=pal[factor(targets$Sample_Source)], dim=c(3,4))
legend("right", legend=levels(factor(targets$Sample_Source)), text.col=pal,
       cex=0.7, bg="white")
    image

    在这里，我们可以看到样本的组别、来源差异相比未标准化和未过滤前，已经减少了很多。

    14.2、重新评估是否已经消除样本间的差异（方法二：ChAMP包）

    bVals <- getBeta(mSetSqFlt)
champ.SVD(beta = bVals ,
              rgSet=NULL,
              pd=targets,
              RGEffect=FALSE,
              PDFplot=TRUE,
              Rplot=TRUE,
              resultsDir="./CHAMP_SVDimages/")
    image

    计算原理是：先对甲基化beta值做主成分分析，对每个主成分和变量（比如本例中的sample_label,sample_group,sample_source，ID，Array等）进行kruskal.test检验，确定两组或多组的中位数是否存在差异。

    如果存在差异，说明变量和甲基化beta值存在相关性，也就是说，变量不能被很好的校正掉，那么，将这些没有被很好校正掉的甲基化数值进行后续分析的话，就很容易产生假阳性。

    从上面截图可以看到，sample_label,sample_group,sample_source这几个变量与甲基化主成分显著相关，后续做EWAS分析时应将他们作为协变量纳入分析中，或者用ChAMP包的champ.runCombat函数或者minfi包的sva函数将这些变量进行校正。

    15、提取M值和beta值

    提取M值（mVals）和beta值（bVals）:

    mVals <- getM(mSetSqFlt)
head(mVals[,1:5])
bVals <- getBeta(mSetSqFlt)
head(bVals[,1:5])

    对M值（mVals）和beta值（bVals）进行画图：

    par(mfrow=c(1,2))
densityPlot(bVals, sampGroups=targets$Sample_Group, main="Beta values", 
            legend=FALSE, xlab="Beta values")
legend("top", legend = levels(factor(targets$Sample_Group)), 
       text.col=brewer.pal(8,"Dark2"))
densityPlot(mVals, sampGroups=targets$Sample_Group, main="M-values", 
            legend=FALSE, xlab="M values")
legend("topleft", legend = levels(factor(targets$Sample_Group)), 
       text.col=brewer.pal(8,"Dark2"))
    image

    收获美美的双峰！

    上面的教程大部分是基于minfi包展开的。

    实际上，除了minfi包，ChAMP包也可以完成这个工作，ChAMP包更简单。直接四个函数搞定。

    如下截图所示。

    这里我就不展开讲了，原理跟minfi包一样的，只不过ChAMP包把它封装好了。

    champ.filter(beta=myImport$beta,
             M=NULL,
             pd=myImport$pd,
             intensity=NULL,
             Meth=NULL,
             UnMeth=NULL,
             detP=NULL,
             beadcount=NULL,
             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,
             fixOutlier = TRUE,
             arraytype = "EPIC")

champ.QC(beta = myLoad$beta,
             pheno=myLoad$pd$Sample_Group,
             mdsPlot=TRUE,
             densityPlot=TRUE,
             dendrogram=TRUE,
             PDFplot=TRUE,
             Rplot=TRUE,
             Feature.sel="None",
             resultsDir="./CHAMP_QCimages/")

champ.norm(beta=myLoad$beta,
               rgSet=myLoad$rgSet,
               mset=myLoad$mset,
               resultsDir="./CHAMP_Normalization/",
               method="BMIQ",
               plotBMIQ=FALSE,
               arraytype="EPIC",
               cores=3)

champ.SVD(beta = myNorm,
              rgSet=NULL,
              pd=myLoad$pd,
              RGEffect=FALSE,
              PDFplot=TRUE,
              Rplot=TRUE,
              resultsDir="./CHAMP_SVDimages/")

    甲基化QC工作到此结束啦~

    16、总结

    minfi流程多、繁琐，胜在轻巧，按着流程走，一般不会出现什么报错。

    ChAMP包方便，但如果数据多的话，对电脑的配置要求也很高，我跑3000个样本时，256G，32cpu核是带不动的，经常跑着跑着就被kill了。用几百个样本跑时，就很顺利。

    17、致谢

    感谢健明分享的甲基化分析入门练习：甲基化芯片的一般分析流程

    建议各位刚入门甲基化的同学们可以看看健明在B站的视频，讲的很详细。

    image

    相关文章

      网友评论

        本文标题：甲基化样本和CpG位点QC的总流程（450k和850k）

        本文链接：https://www.haomeiwen.com/subject/qqvkfktx.html

        延伸阅读

        深度阅读

        • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

        栏目导航

        热点阅读

        DNA甲基化

        R语言训练

        关于我们|服务条款|联系我们|甲基化样本和CpG位点QC的总流程（450k和850k）|投稿指南|网站地图|RSS订阅|排版工具|手机版

        提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

        Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

        备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

        本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

        所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！