QTL分析是進行基因精細定位和克隆的基礎,今天小編教大家使用R包" qtl "進行QTL分析。
在開始分析前,我們需要準備兩個輸入檔案:基因型和表型檔案。
基因型檔案:
表型檔案:
基因型和表型檔案均儲存為逗号分隔的csv檔案。
準備好兩個輸入檔案後,我們就可以開始分析啦!
## 安裝R包
install.packages("qtl")
## 加載R包
library("qtl")
## 導入基因型和表型資料
sug <- read.cross("csvs", ".", "gen.csv", "phe.csv")
## 檢視輸入檔案相關資訊
summary(sug)
複制
此外,還有一些函數可以統計對應的資訊。
## 樣本數
nind(sug)
## 染色體數
nchr(sug)
## 标記數
totmar(sug)
## 每個染色體上的标記數
nmar(sug)
## 表型數
nphe(sug)
除了文字資訊外,我們還可以用圖來展示這些資訊。
複制
plot(sug)
複制
這三張圖分别展示了缺失的基因型資料,遺傳圖譜和表型資料分布。
也可以單獨展示這三張圖。
## 展示缺失基因型資料(黑色為缺失的基因型)
plotMissing(sug)
複制
## 繪制遺傳圖譜
plotMap(sug)
複制
## 繪制表型分布直方圖
plotPheno(sug, pheno.col=1)
複制
## 計算基因型機率
sug <- calc.genoprob(sug, step=1)
## 使用預設方法進行single-QTL全基因組掃描
out.em <- scanone(sug)
## 檢視掃描結果
summary(out.em)
## 挑選LOD > 3的結果
summary(out.em, threshold=3)
## 展示結果
plot(out.em)
複制
我們可以看到7号和15号染色體上各有一個顯著的峰。
## 使用Haley-Knott回歸方法進行全基因組掃描
out.hk <- scanone(sug, method="hk")
## 使用Multiple imputation法進行全基因組掃描
sug <- sim.geno(sug, step=1, n.draws=64)
out.imp <- scanone(sug, method="imp")
## 比較三種方法結果的差異
plot(out.em, out.hk, out.imp, col=c("blue", "red", "green"))
複制
我們可以看到,三種方法的結果并沒有明顯差異。
## 進行1000次Permutation test
operm <- scanone(sug, method="hk", n.perm=1000)
## 獲得顯著性門檻值
summary(operm, alpha=c(0.05, 0.2))
複制
## 從掃描結果中挑選顯著的位點
summary(out.hk, perms=operm, alpha=0.2, pvalues=TRUE)
複制
接下來,我們需要估計QTL區間。因為我們通過LOD值過濾後的QTL位點位于7号和15号染色體上,是以我們首先對7号染色體上的QTL區間的進行估計。
## 獲得7号染色體1.5倍LOD區間和95%貝葉斯區間
lodint(out.hk, chr=7, drop=1.5)
bayesint(out.hk, chr=7, prob=0.95)
複制
第一行和第三行是區間的範圍,第二行是預測QTL的位置。
## 獲得區間兩側最近的标記
lodint(out.hk, chr=7, expandtomarkers=TRUE)
bayesint(out.hk, chr=7, expandtomarkers=TRUE)
複制
## 獲得離QTL最近的标記
max(out.hk)
mar <- find.marker(sug, chr=7, pos=47.7)
## 統計不同基因型個體的表型
plotPXG(sug, marker=mar)
複制
紅色的點表示基因型進行過填補的個體。
## 統計不同基因型個體的效應
effectplot(sug, mname1=mar)
複制
Tips:我們常說的LOD值=log10 (L1/L0) ,L1指該位點含QTL的機率,L0指該位點不含QTL的機率。LOD值為3表示該位點含QLT的機率是不含QTL機率的1000倍。
參考資料:
https://rqtl.org/tutorials/rqtltour2.pdf