イルミナのDNAメチルのmicroarray：Illumina HumanMethylation 450 / EPIC Beadchip（pdf；サイト）の解析。


１。pdファイルの準備
元となるidatファイルをR packageのminfiで読み込んで、ChAMP で処理して、という形で順次解析していくが、それに当たっては各症例のデータ表（Disease or Control、年齢、性別、データ所在フォルダ、などメタデータに相当）であるpdファイルをきちんと形式を守って準備しておく必要がある。

（pdファイル例）


上記は一例で、data.frame形式、一列が一つのmicroarrayサンプルに相当しており、Sample_Groupが状態ラベルを表し（ここではAD（＝Alzheimer）or NC（= Normal Control））、SlideはBatch番号に相当する。元々のidatファイルは　上記では"~/Desktop/Infinium450k_ChAMP"というフォルダに収められていて、このフォルダ内のSlide名フォルダ内に、Slide名_Array名_Grn.idat あるいは　Slide名_Array名_Red.idat という名前のidatファイルが入っているものとする。（＊例えば、"~/Desktop/Infinium450k_ChAMP/8667045130/8667045130_R02C01_Grn.idat"のように）
（なお、idatファイルの格納構造はpdファイルのBasename列できちんと指定できていればなんでも良い。"~/Desktop/idat/~"とかでも構わない。）


２。idatファイルの読み込み

pdファイルを"pd_file"という名前のdata.frameにしてRstudio内のオブジェクトとして持っておくとする。（＊pdファイル例のように、"Slide名_Array名" というデータIDを作って、それを行名にもしておくこと）

まず{minfi}など、必要なパッケージを呼び出す。（まだいれていない場合はinstallしておく）

> library(minfi); library(minfiData); library(sva)

まずidatファイルからRGSetに変換する（＊ pd_file$Basenameでidatファイルの場所が指定されていることが前提）

> RGSet <- read.metharray.exp(targets = pd_file )

次にMSetに変換する（preprocessRawなど色々設定はあるので解説参照）

> MSet <- preprocessNoob(RGSet)

クオリティcheckを行う

> plotQC( getQC(MSet) )

そうすると以下のようなplotが出てくる。bad qualityのものは点線より左下に赤くplotされる。
概ね、横軸＋縦軸 が21以下くらいだとbad判定される。
ここでは全部quality OKなので次に進む。

もし引っかかったら、

qc_df <- (as.data.frame(getQC(MSet)))
qc_df$sum <- qc_df[,1]+qc_df[,2]
qc_df$good <- ifelse(qc_df$sum >= 21 ,1L,0L)
pd_file <- pd_file[rownames(subset(qc_df, good==1)), ]

としてgood qualityのものだけpd_fileに残して、読み込むところからやり直せば良い。

続いて、RSetを経由して、β値行列を得る。

> RSet <- ratioConvert(MSet, what = "both", keepCN = TRUE)
> beta_matrix <- getBeta(RSet)

これで解析元となるβ値行列が得られた。



３。前処理

ここからは{ChAMP}を利用する。
まず、SNPや性染色体に当たるCpGサイト、detection Pの悪いCpGサイトなどをフィルタリングで除外する。（＊arraytypeは"450K" or "EPIC" のみ対応）

> library(ChAMP)
> myFilter <- champ.filter(beta=beta_matrix, pd=pd_file , arraytype = "450K")
> myFilter_df <- as.data.frame(myFilter$beta)
> myFilter_pd <- myFilter$pd

続いて、normalizeを行う。
> myNorm <- champ.norm( myFilter_df , method="BMIQ", plotBMIQ=F, arraytype="450K" ) 
methodは"PBC", "SWAN"など複数あり、Pipeline参照のこと。基本的には結果がすごく大きく変わることはないが、BMIQでやるのが一般的ではある様子。ただ時間がかかるので、結果を見るのを優先したい時はPBCなどが速くて良い。

さらに、batch-effectについて検証する。

> champ.SVD(beta = myNorm , pd= myFilter_pd[,c("Slide","sex")]) 

"Slide"というバッチによって結果が大きく異なるかだけみる。有意な結果なら補正が必要で、

library(sva)
myCombat <- ComBat(dat=myNorm , 
                   batch= myFilter_pd$Slide,
                   mod =  myFilter_pd[,c("status","genderF","age")]  )

上記でComBatをかけると、batch-effectを補正しつつ、modで指定した項目についてunexpectedlyに過剰補正してしまわないようにできる。年齢・性別を入れるかはともかく、disease statusについては指定しないと群間差が小さめの場合、差が消えてしまったりするので指定は必須。（逆にメチルデータで予測したい場合は指定してはいけないということになる。この場合はfrozen sva (fsva)などで対応できるかもしれない。。。）

続いて、群間のデータのばらつきを可視化で把握する。

> mdsPlot(  myCombat , sampGroups= myFilter_pd$status
          ,numPositions = 500, pch=20, legendPos = "bottomleft"  , main = "MDS plot"  )

４。群間比較







 # <pre  class="prettyprint">    </pre>