在上篇帖子中介紹了如何用qtl這個R包進行單QTL定位，但有很多復雜性狀都是由多位點共同控制的。有些位點之間可能存在連鎖，即加性效應或上位性效應，本篇文章將介紹如何通過二維雙QTL基因組掃描，對復雜數(shù)量性狀多位點之間的QTL連鎖或相互作用進行分析。
多QTL聯(lián)合分析有以下三個優(yōu)點：
- 考慮大效應的QTL，從而減少殘差，更好地識別中等效應的QTL；
- 通過比較雙QTL模型與單QTL模型的擬合，可以較好地實現(xiàn)連鎖QTL的分離；
- QTL間的上位性只能通過明確考慮多個QTL的模型來評估。
0.安裝與準備工作
R包安裝與數(shù)據(jù)格式在上一篇帖子中有詳細介紹：
1.QTL定位：Rqtl—— Single-QTL analysis - 簡書 (jianshu.com)
官方說明書：
https://rqtl.org/tutorials/rqtltour2.pdf
1.讀取數(shù)據(jù)
> library(qtl)
> data <- read.cross("csv", ".", "gen_phe.csv")
Warning messages:
1: In read.cross.csv(dir, file, na.strings, genotypes, estimate.map, :
The following unexpected genotype codes were treated as missing.
|--|
2: In summary.cross(cross) :
Some markers at the same position on chr 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27; use jittermap().
這里看到warning信息，第一條是“--”基因型識別為缺失值，第二條是有一些marker在相同的位置上。因此需要通過jittermap()函數(shù)抖動基因圖譜中的標記位置，這樣就不會有兩個標記相互重疊。
jittermap()
> jit<-jittermap(data)
2. 計算條件QTL基因型概率
根據(jù)可用的標記數(shù)據(jù)，利用 calc.genoprob()函數(shù)計算條件QTL基因型概率，這里使用疏松的網(wǎng)格（step=2）來進行計算，提高計算速度：
> calc <- calc.genoprob(jit, step=2)
3. 二維雙QTL掃描
默認情況下，分析是通過EM算法的最大似然來執(zhí)行的，(method="em")，verbose=FALSE阻止跟蹤信息。
> out2d_em <- scantwo(calc, verbose=FALSE)
- 使用Haley-Knott回歸方法進行2D掃描
> out2d_hk <- scantwo(calc, method="hk")
- 組合檢驗
最好是對觀察到的數(shù)據(jù)使用組合檢驗，這樣的結果將考慮表型分布、標記密度和缺失基因型數(shù)據(jù)，這一步需要耗費很長的時間。
> out2d_perm <- scantwo(calc, method="hk",n.perm=1000)
多性狀雙QTL掃描
因為數(shù)據(jù)量比較大，win下的R內(nèi)存達不到，需要提交Linux，由于沒有服務器的root權限，我自己技術有限，Linux下的Rstudio裝不上，所以我把表型文件拆開，分別計算，分別保存為RData，再轉(zhuǎn)回到win里進行后續(xù)分析。
這里我寫了一個循環(huán)解決這個問題，其中1-20列為表型，21列為樣本的ID，22列之后為基因型，如果和我一樣有多個表型的，可以參照下面的循環(huán)更改。
library('qtl')
raw<-read.csv("gen_phe.csv")
for (i in 1:20) {
T<-raw[,c(i,22:9788)]
T[is.na(T)]<-""
file_name <- paste0("gen_phe_", i, ".csv")
write.csv(T,file = file_name,row.names = F,quote =FALSE)
data <- read.cross("csv", ".", file_name)
jit<-jittermap(data)
calc <- calc.genoprob(jit, step=2)
out2d_hk <- scantwo(calc,method="hk")
file_name <- paste0("out2d_hk_", i, ".RData")
save(out2d_hk, file = file_name)
out2d_perm <- scantwo(calc,method="hk",n.perm=100)
file_name <- paste0("out2d_perm_", i, ".RData")
save(out2d_perm, file = file_name)
}
4. 總結P值
如果是win中進行的計算，直接運行：
> summary(out2d_hk, perms=out2d_perm, alpha=0.2, pvalues=TRUE)
如果是Linux中進行的計算，首先要讀入計算得到的.RData，這里以第一個性狀計算得到的out2d_hk_1.RData和out2d_perm_1.RData為例：
> load("out2d_hk_1.Rdata", temp_env <- new.env())
> out2d_hk_1 <- as.list(temp_env)
> out2d_hk_1<-out2d_hk_1[[1]]
> load("out2d_perm_1.Rdata", temp_env <- new.env())
> out2d_perm_1 <- as.list(temp_env)
> out2d_perm_1<-out2d_perm_1[[1]]
> summary(out2d_hk_1, perms=out2d_perm_1, alpha=0.1, pvalues=TRUE)
以說明書中的結果為例進行解釋，詳細內(nèi)容如下：
https://rqtl.org/tutorials/new_summary_scantwo.pdf
lod.full，lod.fv1和lod.int對應pos1f和pos2f；
lod.add和lod.av1對應pos1a和pos2a。
結果中共提供五個LOD值，分別是：
lod.full： 全模型full model (two QTL plus interaction)的效應值。
lod.fv1：比較j、k染色體上QTL的全模型和單QTL模型的LOD值。允許上位性的可能性。
lod.int：比較j、k染色體上QTL的全模型和加性模型的LOD值。證明j、k兩條染色體上是否存在QTL間的相互作用。
lod.add：加性模型的LOD值。
lod.av1：比較j、k染色體上QTL的加性模型和單QTL模型的LOD值。假設沒有上位性。
官方說明中建議忽略lod.int，而對其余4個閾值使用共同的顯著性水平(α = 5 or 10%)。
因此，在示例結果中，共存在兩對QTL存在相互作用，分別在Chr1、4和Chr6、15之間。
引用轉(zhuǎn)載請注明出處，如有錯誤敬請指出。
