Multifactor Dimensionality Reduction(MDR)という手法がある。意味のある組み合わせを考えたいが、組み合わせ爆発が起こるのでうまくbin化しながら考えるようなイメージだと思った。これを遺伝医学に応用した話。
Am J Reprod Immunol. 2013 Jul 31.
RではMDR, mbmdrなどでできる。元々はBioinformatics. 2003 Feb 12;19(3):376-82.
STEP1
CVするために n-fold CVとして今回は1/10にした。
STEP2
N個のlocusがfactorとしてある。
STEP3
遺伝子座について、
と
を取ってきて、3*3のカテゴリーでcase と controlの人数を考える。
STEP4
リスクが1以上なら high, 1以下なら low, 割り付けが存在しなかったら emptyのラベルに貼り替える。
STEP5
と
の組み合わせの時
のモデルのときの誤差をひたすら計算しておく。
STEP6
検証データで誤差を確認する。
library(mbmdr) data(simSNP) fit <- mbmdr(y=simSNP$Y, data=simSNP[, 3:12], order=2, family=binomial(link=logit), printStep1=TRUE) # print the results of step one
SNP1とSNP2を取ってきて計算した時の結果を示す。SNP1の遺伝子型がAA/Aa/aaのときとSNP2の遺伝子型がBB/Bb/bbのときの9パターンにおいて上のステップの通り計算し、Wald testのp値を返している。
Model: 2 1 SNP2 SNP1 cases controls beta p.value category 0 0 0 9 -2.9904 0.0510640 L 1 0 49 19 1.1286 0.0001103 H 2 0 0 15 -3.5117 0.0182311 L 0 1 50 19 1.1554 0.0000727 H 1 1 0 57 -5.0794 0.0004079 L 2 1 43 30 0.4396 0.0938642 H 0 2 0 14 -3.4397 0.0211741 L 1 2 58 26 1.0056 0.0001219 H 2 2 0 11 -3.1919 0.0348002 L
permutation testでp値が小さかった順にSNPの組み合わせを返す。
print(fit)
SNP1 SNP2 NH betaH WH PH NL betaL WL PL MIN.P SNP2 SNP1 4 6.1135 18.285 1.902e-05 5 -6.1135 18.285 1.902e-05 1.902e-05 SNP8 SNP2 1 0.5195 2.771 9.598e-02 2 -0.7656 7.972 4.750e-03 4.750e-03 SNP8 SNP7 1 0.6271 3.592 5.806e-02 1 -0.8908 6.901 8.614e-03 8.614e-03 SNP10 SNP4 2 0.6802 6.182 1.291e-02 1 -0.7414 3.099 7.836e-02 1.291e-02 SNP10 SNP9 0 NA NA NA 1 -1.1196 6.156 1.309e-02 1.309e-02 SNP10 SNP1 0 NA NA NA 1 -0.7777 5.435 1.974e-02 1.974e-02 SNP6 SNP1 0 NA NA NA 1 -0.5302 4.930 2.640e-02 2.640e-02 SNP8 SNP6 1 0.7379 4.845 2.773e-02 0 NA NA NA 2.773e-02 SNP8 SNP4 1 0.8045 2.940 8.640e-02 1 -0.6804 4.518 3.355e-02 3.355e-02 SNP6 SNP2 0 NA NA NA 1 -0.6804 4.518 3.355e-02 3.355e-02 SNP7 SNP3 1 1.0338 4.457 3.475e-02 0 NA NA NA 3.475e-02 SNP7 SNP2 0 NA NA NA 1 -0.5757 3.609 5.747e-02 5.747e-02 SNP6 SNP4 1 0.8744 3.540 5.989e-02 0 NA NA NA 5.989e-02 SNP8 SNP1 1 0.5572 3.503 6.126e-02 0 NA NA NA 6.126e-02 SNP10 SNP2 1 0.5390 3.132 7.677e-02 0 NA NA NA 7.677e-02 SNP9 SNP6 1 0.8045 2.940 8.640e-02 1 -0.5849 3.089 7.883e-02 7.883e-02 SNP4 SNP2 1 0.5214 3.037 8.137e-02 0 NA NA NA 8.137e-02 SNP5 SNP1 1 0.4901 2.867 9.040e-02 0 NA NA NA 9.040e-02
論文では、決定木を使って疾患リスクを推定している。
library("rpart") classification_tree <- rpart(Y ~ . -X, data = simSNP) #finding the best optimal size of tree cp <- classification_tree$cptable[which.min(classification_tree$cptable[, "xerror"]), "CP"]# select the complexity parameter associated with minimum error Regtree_optimal <- prune(classification_tree, cp=cp) plot(Regtree_optimal, uniform=T, main="pruned regression Tree for trait") text(Regtree_optimal, use.n=TRUE, xpd=TRUE)
