1. 處理數(shù)據(jù)

####初次質(zhì)控
plink --bfile mvp --maf 0.05 --geno 0.1 --mind 0.1 --recode vcf --out mvp_qc1
####Beagle填充
java -Xmx2g -jar beagle.28Sep18.793.jar nthreads=4 gt=mvp_qc1.vcf out=imput.mvp.geno 
####對(duì)填充后文件進(jìn)行再質(zhì)控
plink --vcf imput.mvp.geno.vcf --maf 0.05 --geno 0.1 --mind 0.1 --make-bed --out mvp_qc2

2. 運(yùn)行R包rMVP

#library("devtools") 
#devtools::install_github("xiaolei-lab/rMVP")
install.packages("rMVP")
library(rMVP)

# Step1: 數(shù)據(jù)格式轉(zhuǎn)換
MVP.Data(filePhe = "mvp.phe", fileBed = "mvp_qc2", out = "demo")

# Step2: 讀取轉(zhuǎn)換好的數(shù)據(jù)
pheno <- read.table("demo.phe", header = TRUE)       # 表型文件中第一列為個(gè)體 ID，第二列是要分析的性狀
geno  <- attach.big.matrix("demo.geno.desc")         # rMVP 使用的基因型文件是 .geno.desc 和 .geno.bin 兩個(gè)文件一組，前者儲(chǔ)存的是元數(shù)據(jù)，后者是以二進(jìn)制格式儲(chǔ)存的基因型數(shù)據(jù)。
map   <- read.table("demo.geno.map", header = TRUE)  # map 文件中存儲(chǔ)的是標(biāo)記的位置信息

# Step3: 進(jìn)行分析
imMVP <- MVP(
    phe=pheno,    #輸入表型文件
    geno=geno,    #輸入基因型文件
    map=map,      #輸入map文件
    nPC.GLM=5,    #該參數(shù)設(shè)置的是GLM模型中加入的PC數(shù)量，此處設(shè)置為前5個(gè)PC
    nPC.MLM=3,    #該參數(shù)設(shè)置的是MLM模型中加入的PC數(shù)量，此處設(shè)置為前3個(gè)PC
    nPC.FarmCPU=3,#該參數(shù)設(shè)置的是FarmCPU模型中加入的PC數(shù)量，此處設(shè)置為前3個(gè)PC
    threshold=0.05,#該參數(shù)設(shè)置的是顯著水平線的閾值，采用bonfferoni校正，threshold = 0.05/n, n為marker數(shù)量
    method=c("GLM", "MLM", "FarmCPU") # 模型可以設(shè)置多個(gè)同時(shí)進(jìn)行分析。
)

# 多表型情況下，可如下運(yùn)行
for (i in 2:ncol(pheno)){
    imMVP <- MVP(
        phe=pheno,    #輸入表型文件
        geno=geno,    #輸入基因型文件
        map=map,      #輸入map文件
        nPC.GLM=5,    #該參數(shù)設(shè)置的是GLM模型中加入的PC數(shù)量，此處設(shè)置為前5個(gè)PC
        nPC.MLM=3,    #該參數(shù)設(shè)置的是MLM模型中加入的PC數(shù)量，此處設(shè)置為前3個(gè)PC
        nPC.FarmCPU=3,#該參數(shù)設(shè)置的是FarmCPU模型中加入的PC數(shù)量，此處設(shè)置為前3個(gè)PC
        threshold=0.05,#該參數(shù)設(shè)置的是顯著水平線的閾值，采用bonfferoni校正，threshold = 0.05/n, n為marker數(shù)量
        method=c("GLM", "MLM", "FarmCPU") # 模型可以設(shè)置多個(gè)同時(shí)進(jìn)行分析。
    )
}