這篇應該是甲基化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
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值。
畫每個樣本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")
如下圖所示：
可以看到，只有最后一個樣本的cpg平均P值是超過0.05，也就是說，這個樣本的質量是比較差的，后續應該被剔除掉。
導出質量報道：
qcReport(rgSet, sampNames=targets\(ID, sampGroups=targets\)Sample_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個樣本：
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"))
畫出來的圖如下所示，左邊是未進行標准化的，右邊是標准化以后的。
可見，進行標准化后，樣本間的差異會縮小。
9、查找標准化后數據可能存在的差異來源
這一步是為EWAS做准備的，我們前面進行了標准化，但標准化的數據不代表就可以完全去除樣本批次效應、細胞類型等差異。
如果樣本間存在批次效應等可能的混淆因素，在后續進行EWAS分析時極大可能會產生假陽性。
因此我們需要通過主成分分析對標准化的數據進行可視化。確定可能的混淆因素，並在EWAS分析時進行校正。
10、通過主成分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個樣本被分為三個聚類。
說明即便是前期進行了標准化處理后，樣本間還是存在差異，比如樣本act_naive.5和naive.1就很明顯的在不同的聚類中。
這種差異在進行EWAS分析時是我們不願意看到的，因此后期進行EWAS分析時，應考慮將他們納入協變量中。
11、通過主成分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畫圖后如下所示：
可以看到，不同樣本按照不同顏色很好地被分層了，說明主成分3,4反應的是細胞類型的差異。
同樣的，細胞類型差異在進行EWAS分析時也是我們並不願意看到的，因此他們在進行EWAS分析時應被一起納入協變量中校正掉。
12、探針過濾
前面我們根據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
13、移除包含性染色體的探針
keep <- !(featureNames(mSetSqFlt) %in% ann450k$Name[ann450k$chr %in% c("chrX","chrY")])
table(keep)
mSetSqFlt <- mSetSqFlt[keep,]
14、移除SNP探針
mSetSqFlt <- dropLociWithSnps(mSetSqFlt)
mSetSqFlt
15.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
15.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)
16.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")
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")
在這里，我們可以看到樣本的組別、來源差異相比未標准化和未過濾前，已經減少了很多。
16.2、重新評估是否已經消除樣本間的差異（方法二：ChAMP包）
bVals <- getBeta(mSetSqFlt)
champ.SVD(beta = bVals ,
rgSet=NULL,
pd=targets,
RGEffect=FALSE,
PDFplot=TRUE,
Rplot=TRUE,
resultsDir="./CHAMP_SVDimages/")
計算原理是：先對甲基化beta值做主成分分析，對每個主成分和變量（比如本例中的sample_label,sample_group,sample_source，ID，Array等）進行kruskal.test檢驗，確定兩組或多組的中位數是否存在差異。
如果存在差異，說明變量和甲基化beta值存在相關性，也就是說，變量不能被很好的校正掉，那么，將這些沒有被很好校正掉的甲基化數值進行后續分析的話，就很容易產生假陽性。
從上面截圖可以看到，sample_label,sample_group,sample_source這幾個變量與甲基化主成分顯著相關，后續做EWAS分析時應將他們作為協變量納入分析中，或者用ChAMP包的champ.runCombat函數或者minfi包的sva函數將這些變量進行校正。
17、提取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"))
收獲美美的雙峰！
上面的教程大部分是基於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工作到此結束啦~
18、總結
minfi流程多、繁瑣，勝在輕巧，按着流程走，一般不會出現什么報錯。
ChAMP包方便，但如果數據多的話，對電腦的配置要求也很高，我跑3000個樣本時，256G，32cpu核是帶不動的，經常跑着跑着就被kill了。用幾百個樣本跑時，就很順利。
19、致謝
感謝健明分享的甲基化分析入門練習：甲基化芯片的一般分析流程
建議各位剛入門甲基化的同學們可以看看健明在B站的視頻，講的很詳細。
