library(qtl)								#ライブラリー起動
testmap <- read.cross("csvr", file="IN-RIL.csv", estimate.map=FALSE)	#データ読み込み
testmap <- convert2riself(testmap)					# RI selfとしてデータを処理する


#QTL解析

testmap <- calc.genoprob(testmap, step=1) 	#1cM毎にQTL検索
#interval mapping(IM)を行う。1番目の形質を解析する
out.simhk <- scanone(testmap, pheno.col=1, method="hk") 

#並替え検定により有意水準(pval)を推定する。
operm.hk <- scanone (testmap, n.perm=100, method="hk", pheno.col=1) 
summary(out.simhk, perms=operm.hk, alpha=0.05, pvalues=TRUE)


#composite interval mapping(CIM)を行う
out.cimhk <- cim(testmap, n.marcovar=5, pheno.col=1, method="hk") 
#並べ替え検定を行う
operm.chk <- cim(testmap,n.perm=100, pheno.col=1,n.marcovar=5,method="hk")
#operm.chk <- cim(testmap,n.perm=1000, pheno.col=1,n.marcovar=5,method="hk")


#5%水準で有意なQTLピークを書き出すには、
summary(out.cimhk, perms=operm.chk, alpha=0.05, pvalues=TRUE)


#IMとCIMの結果を表示する
plot(out.simhk, out.cimhk, col=c("blue", "red"),show.marker.names=TRUE)

plot(out.simhk, out.cimhk, col=c("blue", "red"),show.marker.names=TRUE,chr=c(6,8))#有意なQTLが検出された染色体6,8のみを表示する



#twoQTLscan
testmap <- calc.genoprob(testmap, step=1)	#1cM毎にQTLを検索する
out2.hk <- scantwo(testmap, method="hk")	# two QTLのスキャン
summary(out2.hk)	#染色体対毎の最大LODピークを表示


operm2 <- scantwo(testmap, method="hk", n.perm=10, pheno.col=1)
#operm2 <- scantwo(testmap, method="hk", n.perm=1000)  #大量のメモリが必要
summary(out2.hk, perms=operm2, alpha=0.05, pvalues=TRUE)



#QTL fitting
#imputation n.drawsはシミュレーション回数、stepは解像度(cM)
testmap <- sim.geno(testmap, step=1, n.draws=128, err=0.001)
  
chr <- c(6,8)
pos <- c(9,56)
qtl <- makeqtl(testmap, chr, pos)
out.fq <- fitqtl(testmap, qtl=qtl,formula = y~Q1+Q2)
summary(out.fq)