利用R/qtl進行MQM程式碼

# MQM方法
library(qtl)
#permutation所需要的package
library(snow)
fsi_LA1543<-read.cross("csv",dir="d:",file="RqtlLA1543.csv")
#確定data是常態分佈
mqmtestnormal(fsi_LA1543)
transfsi_LA1543 <- transformPheno(fsi_LA1543, transf=log)
mqmtestnormal(transfsi_LA1543)
#把marker之間可能的基因型也算進去，所以出現的QTL不以marker為主，以cM為主
mqmsqt_LA1543<- calc.genoprob(transfsi_LA1543)
#對genetype缺值進行最有可能基因型計算，在strategy中有三種方法可進行
mqm.out<-mqmaugment(mqmsqt_LA1543,strategy="impute",verbose=T)
#進行permutation
mqm_perm<-mqmpermutation(mqm.out,scanfunction=mqmscan,n.perm=5000,plot=T,method="permutation",verbose=T)
mqm_threshold=summary(mqm_perm,alpha=c(0.1,0.05,0.01))
result_perm<-mqmprocesspermutation(mqm_perm)
summary(result_perm)
#進入MQM的QTL偵測
result<-mqmscan(mqm.out)
summary(result,threshold=3)
plot(result)
for(i in 1:3){abline(h=mqm_threshold[i],lty=3)}
summary(result,perm=result_perm,alpha=0.1)
summary(result,perm=result_perm,alpha=0.05)