R/qtl是一個可擴(kuò)展的暮刃、交互的、用于QTL定位的R包爆土，集數(shù)據(jù)分析和繪圖于一體椭懊。特點是使復(fù)雜的QTL定位方法被廣泛使用，使用戶專注于建模而不是計算步势。
官網(wǎng)：
https://rqtl.org/
在網(wǎng)頁下部還有一塊氧猬，是其他常用的QTL定位軟件，都可以了解一下坏瘩。
官方說明書（簡易版）：
https://rqtl.org/tutorials/rqtltour2.pdf
還有詳細(xì)的說明書盅抚，專門對Rqtl的原理進(jìn)行了介紹，遺傳學(xué)知識非常豐富倔矾，有時間的話強(qiáng)烈建議把這本書看完妄均，絕對受益匪淺。
A Guide to QTL Mapping with R/qtl
這篇帖子就只介紹一下最基礎(chǔ)的用法哪自，想了解更多丰包，還是自己去看說明書哈！
0. 準(zhǔn)備工作
安裝：
> install.packages("qtl")
調(diào)用：
> library(qtl)
1. 數(shù)據(jù)格式：
官方各格式示例數(shù)據(jù)：
https://rqtl.org/sampledata/
基因型和表型數(shù)據(jù)整理在一張以逗號分隔的csv文件中壤巷，以樣本ID為界限邑彪，ID前為表型值，ID后為基因型胧华。
前幾列為表型值寄症，一種表型一列；
sex：性別
ID：sample名稱矩动，注意表型有巧，基因型，ID都要對應(yīng)
后幾列為基因型悲没。
其中表型和基因型是必須的篮迎，其余可有可無。
基因型格式：
第一行：maker名
第二行：染色體編號
第三行：maker所在的位置
第四行之后為基因型
這篇帖子中用到的示例數(shù)據(jù)是：
https://rqtl.org/sampledata/listeria.csv
2.讀取數(shù)據(jù)
> listeria.a <- read.cross("csv", ".", "listeria.csv")
3. 數(shù)據(jù)完整性檢查
> summary(listeria.a)
F2 intercross
No. individuals: 120
No. phenotypes: 3
Percent phenotyped: 96.7 100 100
No. chromosomes: 20
Autosomes: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
X chr: X
Total markers: 133
No. markers: 13 6 6 4 13 13 6 6 7 5 6 6 12 4 8 4 4 4 4 2
Percent genotyped: 88.5
Genotypes (%):
Autosomes: AA:25.8 AB:48.9 BB:24.4 not BB:0.0 not AA:0.9
X chromosome: AA:51.7 AB:48.3
從summary中可以看到檀训，共有120個樣本柑潦，19條常染色體享言，一條性染色體峻凫。
有許多簡單的函數(shù)可以提取匯總信息：
> nind(listeria.a) #樣本數(shù)
[1] 120
> nchr(listeria.a) #染色體數(shù)
[1] 20
> totmar(listeria.a) #標(biāo)記數(shù)
[1] 133
> nmar(listeria.a) #每條染色體上的標(biāo)記數(shù)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 X
13 6 6 4 13 13 6 6 7 5 6 6 12 4 8 4 4 4 4 2
> nphe(listeria.a) #表型數(shù)
[1] 3
4. 匯總信息畫圖
> pdf("summary_plot.pdf",width=8,height=10)
> plot(listeria.a)
> dev.off()
上述圖的各個部分可以分別按如下方式得到：
# 缺失基因型數(shù)據(jù)
> plotMissing(listeria.a)
# 繪制遺傳圖譜
> plotMap(listeria.a)
# 繪制表型數(shù)據(jù)直方圖
> plotPheno(listeria.a, pheno.col=1)
5. Single-QTL analysis
首先，利用calc.genoprob函數(shù)計算基因型的概率览露，step參數(shù)設(shè)置步長荧琼，單位是cM，步長確定了后期QTL定位的密度。
> jit<-jittermap(listeria.a)
> sug <- calc.genoprob(jit, step=2)
通過scanone函數(shù)進(jìn)行QTL掃描
> out.em <- scanone(sug)
查看每條染色體上的最大LOD值
> summary(out.em)
chr pos lod
c1.loc81 1 81.00 2.1057
c2.loc36 2 36.00 1.0528
D3M147 3 63.19 1.8446
D4M251 4 68.10 1.1744
c5.loc28 5 28.00 6.7131
D6M15 6 59.37 3.3279
D7M105 7 60.11 0.5715
D8M94 8 0.00 0.6940
D9M328 9 4.22 1.0563
D10M42_ 10 40.71 0.5819
D11M333 11 64.34 0.2801
c12.loc44 12 44.00 2.1008
D13M147 13 26.16 5.8292
c14.loc11 14 11.00 0.0819
c15.loc23 15 23.00 3.1928
c16.loc37 16 37.00 1.2660
c17.loc16 17 16.00 0.6180
D18M186 18 20.90 0.8158
D19M117 19 16.36 0.4822
DXM186 X 0.00 0.7381
設(shè)置閾值命锄，查看大于設(shè)定閾值的染色體堰乔，這里設(shè)置LOD>3
LOD=log10 (L1/L0) ，L1指該位點含QTL的概率脐恩，L0指該位點不含QTL的概率镐侯。LOD值為3表示該位點含QLT的概率是不含QTL概率的1000倍。
> summary(out.em, threshold=3)
chr pos lod
c5.loc28 5 28.0 6.71
D6M15 6 59.4 3.33
D13M147 13 26.2 5.83
c15.loc23 15 23.0 3.19
繪制結(jié)果圖
> plot(out.em)
使用Haley-Knott回歸方法進(jìn)行掃描
> out.hk <- scanone(sug, method="hk")
同時繪制默認(rèn)方法驶冒，和Haley-Knott回歸方法的LOD值圖
> plot(out.em, out.hk, col=c("blue", "red"))
使用Multiple imputation法進(jìn)行掃描
> sug2 <- sim.geno(sug, step=1, n.draws=64)
> out.imp <- scanone(sug2, method="imp")
三種方法苟翻，只繪制Chr7,15的LOD值
> plot(out.em, out.hk, out.imp, col=c("blue", "red", "green"), chr=c(7,15))
6. 組合檢驗 Permutation tests
> operm <- scanone(sug, method="hk", n.perm=1000)
顯著性閾值可以通過summary函數(shù)得到
> summary(operm)
LOD thresholds (1000 permutations)
lod
5% 3.62
10% 3.17
> summary(operm, alpha=c(0.05, 0.2))
LOD thresholds (1000 permutations)
lod
5% 3.62
20% 2.77
組合檢驗結(jié)果可以與掃描結(jié)果一起使用，自動計算顯著性閾值和p值
> summary(out.hk, perms=operm, alpha=0.2, pvalues=TRUE)
chr pos lod pval
c5.loc28 5 28.0 6.68 0.000
D6M15 6 59.4 3.33 0.081
D13M147 13 26.2 5.83 0.001
c15.loc23 15 23.0 3.20 0.095
7. QTL區(qū)間估計
通過前面LOD值閾值設(shè)定骗污，確定QTL位點位于5崇猫、6、13號和15號染色體上需忿，這里只對5號染色體上的QTL區(qū)間的進(jìn)行估計诅炉。
獲得7號染色體1.5倍LOD區(qū)間和95%貝葉斯區(qū)間
> lodint(out.hk, chr=5, drop=1.5)
chr pos lod
c5.loc16 5 16 5.145233
c5.loc28 5 28 6.682493
c5.loc38 5 38 4.832705
> bayesint(out.hk, chr=5, prob=0.95)
chr pos lod
c5.loc19 5 19 5.575862
c5.loc28 5 28 6.682493
c5.loc35 5 35 5.723181
其中第一行為區(qū)間開始位置 ，第三行為區(qū)間結(jié)束位置屋厘，中間一行是QTL的預(yù)估位置涕烧。
獲得區(qū)間兩側(cè)最近的標(biāo)記
> lodint(out.hk, chr=5, expandtomarkers=TRUE)
chr pos lod
D5M232 5 6.10396 1.947197
c5.loc28 5 28.00000 6.682493
D5M338 5 38.06807 4.805719
> bayesint(out.hk, chr=5, expandtomarkers=TRUE)
chr pos lod
D5M232 5 6.10396 1.947197
c5.loc28 5 28.00000 6.682493
D5M338 5 38.06807 4.805719
8. QTL效應(yīng)計算
這里以上面檢測到的Chr5，position=28.00的位置為例：
基因型分布
> mar <- find.marker(sug, chr=5, pos=28.0)
> plotPXG(sug, marker=mar)
統(tǒng)計不同基因型效應(yīng)
> effectplot(sug, mname1=mar)
等同于
> effectplot(sug, mname1="5@28.0")
5@28.0即為Chr5汗洒，position=28.0澈魄。
統(tǒng)計兩個位點間的聯(lián)合效應(yīng)分析
> effectplot(sug, mname1="5@28.0", mname2="6@59.37")
9. 其他表型的QTL定位
通過pheno.col函數(shù)設(shè)置分析的表型是第幾個：
> out.hr <- scanone(sug, pheno.col=2, method="hk")
一次性對所有表型進(jìn)行掃描，比如第一列到第四列：
> out.all <- scanone(sug, pheno.col=1:4, method="hk")
> summary(out.all, threshold=3)
引用轉(zhuǎn)載請注明出處仲翎，如有錯誤敬請指出痹扇。
