甲基化晶片數據的一些質控指標

前面我們介紹了一些背景知識，主要是理解什麼是DNA甲基化，為什麼要檢測它，以及晶片和測序兩個方向的DNA甲基化檢測技術。具體介紹在：甲基化的一些基礎知識，也了解了甲基化晶片的一般分析流程 。

然後下載了自己感興趣的項目的每個樣本的idat原始文件，也可以簡單通過minfi包或者champ處理它們拿到一個對象。

成功下載了數據而且導入了R裡面，按照道理應該是要直奔主題搞差異分析啦，但是呢，我強調過很多次，甲基化信號值矩陣是有它的特殊性，雖然分析流程與mRNA那樣的表達晶片總體上是一致，幾個細節還是要注意的。

再次強調，需要區分：beta matrix, M matrix, intensity matrix from same .idat file ，自行查看我們生信技能樹往期教程哦！還有 beadcount, detection P matrix, or Meth UnMeth matrix 這些名詞也需要了解。

    beta        One single beta matrix to do filtering. (default = myImport$beta). 
  M        One single M matrix to do filtering. (default = NULL).  
  pd        pd file related to this beta matrix, suggest provided, because maybe filtering would be on pd file. (default = myImport$pd)    
  intensity        intensity matrix. (default = NULL). 
  Meth        Methylated matrix. (default = NULL).    
  UnMeth        UnMethylated matrix. (default = NULL).  
  detP        Detected P value matrix for corresponding beta matrix, it MUST be 100% corresponding, which can be ignored if you don't have.(default = NULL)   
  beadcount        Beadcount information for Green and Red Channal, need for filterBeads.(default = NULL)

就是因為概念名詞太多，所以很多人就不願意去仔細理解甲基化晶片數據。

從ExpressionSet對象拿到甲基化信號值矩陣

通常我們是從GEO資料庫下載甲基化信號值矩陣文件，使用getGEO函數導入成為ExpressionSet對象，如下：

require(GEOquery)
require(Biobase)
eset <- getGEO("GSE68777",destdir = './',AnnotGPL = T,getGPL = F)
beta.m <- exprs(eset[[1]])
beta.m[1:4,1:4]
save(eset,file = 'GSE68777_geoquery_eset.Rdata')

## 順便把臨床信息製作一下，下面的代碼，具體每一個項目都是需要修改的哦
load(file = 'GSE68777_geoquery_eset.Rdata')
pD.all <- pData(eset[[1]])
pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1.1", "characteristics_ch1.2")]
head(pD)
names(pD)[c(3,4)] <- c("group", "sex")
pD$group <- sub("^diagnosis: ", "", pD$group)
pD$sex <- sub("^Sex: ", "", pD$sex)
head(pD)
save(pD,file = 'GSE68777_geoquery_pd.Rdata')

其中getGEO函數拿到的是ExpressionSet對象，所以可以使用exprs函數來獲取甲基化信號值矩陣，那個beta.m就是後續需要質控的。

從minfi的對象拿到甲基化信號值矩陣

使用minfi包的read.metharray.exp函數讀取，前面下載的該數據集的RAW.tar 裡面的各個樣本的idat文件，就被批量加載到R裡面，代碼如下：

# BiocManager::install("minfi",ask = F,update = F)
library("minfi")
# minfi 無法讀取壓縮的idat文件，所以需要解壓
list.files("idat", pattern = "idat")
rgSet <- read.metharray.exp("idat")
rgSet
save(rgSet,file = 'GSE68777_minfi_rgSet.Rdata')

其中 idat 是一個文件夾，裡面有很多idat的晶片原始數據哦。

你需要學習 rgSet 所表示的對象 class: RGChannelSet ，才能進行後續處理，其實也就是學習如何從裡面拿到甲基化信號矩陣和表型信息。

從ChAMP的對象拿到甲基化信號值矩陣

同樣的是可以讀取數據集的RAW.tar 裡面的各個樣本的idat文件，唯一的區別是需要對你的項目製作一個csv表型文件，示例如下：

[Header],,,,,,,
Investigator Name,,,,,,,
Project Name,,,,,,,
Experiment Name,,,,,,,
Date,3/18/2012,,,,,,
,,,,,,,
[Data],,,,,,,
Sample_Name,Sample_Plate,Sample_Group,Pool_ID,Project,Sample_Well,Sentrix_ID,Sentrix_Position
C1,,C,,,E09,7990895118,R03C02
C2,,C,,,G09,7990895118,R05C02
C3,,C,,,E02,9247377086,R01C01
C4,,C,,,F02,9247377086,R02C01
T1,,T,,,B09,7766130112,R06C01
T2,,T,,,C09,7766130112,R01C02
T3,,T,,,E08,7990895118,R01C01
T4,,T,,,C09,7990895118,R01C02

大多數人的R語言學的不咋地，所以可以考慮在Excel裡面整理這個csv表格咯。但是作者一直強調：the user MUST understand their pd file well to achieve a successful analysis.

最重要的是 Sample_Group 列，表明你需要把你的甲基化信號矩陣如何分組後續進行差異分析。

其次是 Sentrix_ID,Sentrix_Position兩列，決定你的idat文件名前綴。我代碼如下：

rm(list = ls())   
options(stringsAsFactors = F)

# BiocManager::install("ChAMP",ask = F,update = F)
library("ChAMP")
require(GEOquery)
require(Biobase)
# 臨床信息來自於前面的getGEO函數導 
load(file = 'GSE68777_geoquery_pd.Rdata')
head(pD)

fs=list.files("idat", pattern = "idat") 
library(stringr)
fsd=unique(str_split(fs,'_',simplify = T)[,1:3])
pd=cbind(fsd,pD[fsd[,1],])

cln=strsplit('Sample_Name,Sample_Plate,Sample_Group,Pool_ID,Project,Sample_Well,Sentrix_ID,Sentrix_Position',
         ',')[[1]]
pd=pd[,c(4,2,6,7,3,1,1,4)]
colnames(pd)=cln
pd[,2]=NA;pd[,4]=NA;pd[,5]=NA
write.table(pd,file = 'idat/sample_sheet.csv',row.names = F,
            col.names = T,quote = F,sep = ',')

myLoad <- champ.load('idat/',arraytype="450K")
save(myLoad,file = 'GSE68777_champ_load.Rdata')

看起來比前面兩個方法都複雜很多哦，主要是因為我使用R來製作樣本信息表格啦！

我們現在存儲了3個數據對象，接下來的質控就針對這3個分別操作哦！

156M Feb  8 15:34 GSE68777_champ_load.Rdata
141M Feb  8 13:45 GSE68777_geoquery_eset.Rdata
607B Feb  8 15:15 GSE68777_geoquery_pd.Rdata
114M Feb  7 23:30 GSE68777_minfi_rgSet.Rdata

因為數據都很大，所以不方便GitHub共享，大家可以自行下載idat的晶片原始數據文件，然後走我裡面的代碼，肯定是成功的。

除非是你三五年後看到這個教程，有可能R包更新導致某些函數會失效。當然了，這也就是給你提個醒咯，函數和代碼是有可能失效的哈。

質控的指標如果是拿到甲基化信號值矩陣表達矩陣

如果是mRNA表達矩陣，我們通常是看3張圖，代碼裡面我挑選了top1000的sd基因繪製熱圖，然後就可以分辨出來自己處理的數據集裡面的樣本分組是否合理啦。其實這個熱圖差不多等價於PCA分析的圖，被我稱為表達矩陣下遊分析標準3圖！詳見：你確定你的差異基因找對了嗎？ ，就是下面的3張圖：

左邊的熱圖，說明我們實驗的兩個分組，normal和npc的很多基因表達量是有明顯差異的

中間的PCA圖，說明我們的normal和npc兩個分組非常明顯的差異

右邊的層次聚類也是如此，說明我們的normal和npc兩個分組非常明顯的差異

PS：如果你的轉錄組實驗分析報告沒有這三張圖，就把我們生信技能樹的這篇教程甩在他臉上，讓他瞧瞧，學習下轉錄組數據分析。

同樣的道理，針對甲基化信號值矩陣表達矩陣也可以繪製標準3張圖咯。首先看看簡單的boxplot和密度圖，mds圖。

# 使用 wateRmelon 進行 歸一化
library("wateRmelon")
beta.m=beta.m[rowMeans(beta.m)>0.005,]
pdf(file="rawBox.pdf")
boxplot(beta.m,col = "blue",xaxt = "n",outline = F)
dev.off()
beta.m = betaqn(beta.m)
pdf(file="normalBox.pdf")
boxplot(beta.m,col = "red",xaxt = "n",outline = F)
dev.off()

# 然後進行簡單的QC
head(pD)
pdf(file="densityBeanPlot.pdf")
par(oma=c(2,10,2,2))
densityBeanPlot(beta.m, sampGroups = pD$group)
dev.off()
pdf(file="mdsPlot.pdf")
mdsPlot(beta.m, numPositions = 1000, sampGroups = pD$group)
dev.off()
# 後續針對 beta.m 進行差異分析, 比如 minfi 包
grset=makeGenomicRatioSetFromMatrix(beta.m,what="Beta")
M = getM(grset)
# 因為甲基化晶片是450K或者850K，幾十萬行的甲基化位點，統計檢驗通常很慢。
dmp <- dmpFinder(M, pheno=pD$group, type="categorical")
dmpDiff=dmp[(dmp$qval<0.05) & (is.na(dmp$qval)==F),]

上面代碼的4張圖，僅僅是展現了wateRmelon進行歸一化前後表達量的均一性，那個mds圖，就是pca圖的類似，也能說明這個項目的分組並不是數據最大的差異，好奇怪哦。

然後看看我們生信技能樹獨創的表達矩陣標準3張圖，其中pca圖類似於上面的mds圖哈，具體統計學原理我這裡略過。

group_list=pD$group
if(T){


  dat=t(beta.m)
  dat[1:4,1:4] 
  library("FactoMineR")#畫主成分分析圖需要加載這兩個包
  library("factoextra")  
  # 因為甲基化晶片是450K或者850K，幾十萬行的甲基化位點，所以PCA不會太快
  dat.pca <- PCA(dat , graph = FALSE) 
  fviz_pca_ind(dat.pca,
               geom.ind = "point", # show points only (nbut not "text")
               col.ind = group_list, # color by groups
               # palette = c("#00AFBB", "#E7B800"),
               addEllipses = TRUE, # Concentration ellipses
               legend.title = "Groups"
  )
  ggsave('all_samples_PCA.png')

  dat=beta.m
  dat[1:4,1:4] 
  cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，從小到大排序，取最大的1000個
  library(pheatmap)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #對那些提取出來的1000個基因所在的每一行取出，組合起來為一個新的表達矩陣
  n=t(scale(t(dat[cg,]))) # 'scale'可以對log-ratio數值進行歸一化
  n[n>2]=2 
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n)  
  pheatmap(n,show_colnames =F,show_rownames = F,
           annotation_col=ac,filename = 'heatmap_top1000_sd.png')
  dev.off()

  exprSet=beta.m
  pheatmap::pheatmap(cor(exprSet)) 
  # 組內的樣本的相似性應該是要高於組間的！
  colD=data.frame(group_list=group_list)
  rownames(colD)=colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet),
                     annotation_col = colD,
                     show_rownames = F,
                     filename = 'cor_all.png')
  dev.off() 
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  # M=cor(log2(exprSet+1)) 
  M=cor(exprSet)
  pheatmap::pheatmap(M,annotation_col = colD)
  pheatmap::pheatmap(M,
                     show_rownames = F,
                     annotation_col = colD,
                     filename = 'cor_top500.png')
  dev.off() 

}

因為這個甲基化晶片的信號值矩陣問題，我們的圖如下：

如果你多做一點項目，就會有自己的認知啦。之所以你拿450K全部晶片來做PCA，熱圖，都是生物學分組不明顯，因為性別差異，在甲基化晶片效應非常明顯，如下圖：

可以很明顯的看到，top1000的sd的甲基化探針，主要是在性染色體上面差異巨大！

但是呢，你有沒有覺得很煩躁，每次都有做這麼多圖，最後其實也用不上，僅僅是自己看看罷了。

如果是minfi的對象

主要是參考：http://master.bioconductor.org/packages/release/workflows/vignettes/methylationArrayAnalysis/inst/doc/methylationArrayAnalysis.html

同樣是進行一系列質控，可以把450K晶片過濾到420K左右：

detP <- detectionP(rgSet)
head(detP)
mSetSq <- preprocessQuantile(rgSet)  
detP <- detP[match(featureNames(mSetSq),rownames(detP)),]  
keep <- rowSums(detP < 0.01) == ncol(mSetSq)  
table(keep)
# 可以看到，這個步驟過濾的探針很有限
mSetSqFlt <- mSetSq[keep,]
mSetSqFlt
# 繼續過濾性別相關探針，非常重要
ann450k <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
head(ann450k)
keep <- !(featureNames(mSetSqFlt) %in% ann450k$Name[ann450k$chr %in% 
                                                      c("chrX","chrY")])
table(keep)
mSetSqFlt <- mSetSqFlt[keep,]

# remove probes with SNPs at CpG site
mSetSqFlt <- dropLociWithSnps(mSetSqFlt)
mSetSqFlt
# BiocManager::install("methylationArrayAnalysis",ask = F,update = F)
dataDirectory <- system.file("extdata", package = "methylationArrayAnalysis")
# list the files
# exclude cross reactive probes 
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
# 全部過濾完成後，還剩下420K的甲基化位點

load(file = 'GSE68777_geoquery_pd.Rdata')
head(pD)
group_list=pD$group 
qcReport(rgSet,   sampGroups=group_list, 
         pdf="minifi_raw_qcReport.pdf")

save(mSetSqFlt,file =  'GSE68777_minfi_mSetSqFlt.Rdata')
# calculate M-values for statistical analysis
mVals <- getM(mSetSqFlt)
head(mVals[,1:5])
bVals <- getBeta(mSetSqFlt)
head(bVals[,1:5])

同樣的，有了甲基化信號值矩陣後，就可以走我們的標準3張圖。

如果是champ的對象

主要是參考：http://www.bioconductor.org/packages/release/bioc/vignettes/ChAMP/inst/doc/ChAMP.html

進行一系列質控即可，值得注意的是，champ在讀取idat晶片原始文件的時候，默認就做了6個步驟的過濾，如下：

First filter is for probes with detection p-value (default > 0.01).

Second, ChAMP will filter out probes with <3 beads in at least 5% of samples per probe.

Third, ChAMP will by default filter out all non-CpG probes contained in your dataset.

Fourth, by default ChAMP will filter all SNP-related probes.

Fifth, by default setting, ChAMP will filter all multi-hit probes.

Sixth, ChAMP will filter out all probes located in chromosome X and Y.

所以，通常一個450K的晶片，加載到R裡面就只有400K位點啦。可以拿到過濾後的信號值矩陣自己走我們標準的3張圖質控策略，也可以使用champ自帶的質控函數。

champ.norm 提供了四種方法：BMIQ, SWAN1, PBC2 and FunctionalNormliazation。默認的方法是BMIQ, 且BMIQ對850K的標準化方法更好一點!

代碼異常簡單：

load('GSE68777_champ_load.Rdata')
myLoad$pd
champ.QC() ## 輸出qc的圖表
myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
dim(myNorm) 

# 原來的450K經過質控過濾後是400K啦
beta.m=myNorm
dim(beta.m)
group_list=pD$group 

同樣的，有了甲基化信號值矩陣後，就可以走我們的標準3張圖。

如果是TCGA資料庫下載的甲基化信號值矩陣

其實跟從GEO資料庫下載甲基化信號值矩陣文件沒什麼區別哈，通常也推薦使用 ChAMP 流程咯。

A significant improvement for ChAMP is that users can also conduct full analysis even if they are not starting from raw .idat files. As long as you have a methylation beta matrix and the corresponding phenotypes (pd) file, you can conduct nearly all of the ChAMP analysis. This makes analysis easier for users who have beta-values only but not original idat files, for example if users obtain data from TCGA or GEO.

強烈建議你使用ChAMP 流程的測試例子，幾行代碼就搞定甲基化晶片數據分析全部環節。

data(EPICSimData)
CpG.GUI(arraytype="EPIC")
champ.QC() # Alternatively QC.GUI(arraytype="EPIC")
myNorm <- champ.norm(arraytype="EPIC")
champ.SVD()
# If Batch detected, run champ.runCombat() here.This data is not suitable.
myDMP <- champ.DMP(arraytype="EPIC")
DMP.GUI()
myDMR <- champ.DMR()
DMR.GUI()
myDMR <- champ.DMR(arraytype="EPIC")
DMR.GUI(arraytype="EPIC")
myBlock <- champ.Block(arraytype="EPIC")
Block.GUI(arraytype="EPIC") # For this simulation data, not Differential Methylation Block is detected.
myGSEA <- champ.GSEA(arraytype="EPIC")
myEpiMod <- champ.EpiMod(arraytype="EPIC")