在上篇帖子中介绍了如何用qtl这个R包进行单QTL定位，但有很多复杂性状都是由多位点共同控制的。有些位点之间可能存在连锁，即加性效应或上位性效应，本篇文章将介绍如何通过二维双QTL基因组扫描，对复杂数量性状多位点之间的QTL连锁或相互作用进行分析。
多QTL联合分析有以下三个优点：
- 考虑大效应的QTL，从而减少残差，更好地识别中等效应的QTL；
- 通过比较双QTL模型与单QTL模型的拟合，可以较好地实现连锁QTL的分离；
- QTL间的上位性只能通过明确考虑多个QTL的模型来评估。
0.安装与准备工作
R包安装与数据格式在上一篇帖子中有详细介绍：
1.QTL定位：Rqtl—— Single-QTL analysis - 简书 (jianshu.com)
官方说明书：
https://rqtl.org/tutorials/rqtltour2.pdf
1.读取数据
> 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()函数抖动基因图谱中的标记位置，这样就不会有两个标记相互重叠。
jittermap()
> jit<-jittermap(data)
2. 计算条件QTL基因型概率
根据可用的标记数据，利用 calc.genoprob()函数计算条件QTL基因型概率，这里使用疏松的网格（step=2）来进行计算，提高计算速度：
> calc <- calc.genoprob(jit, step=2)
3. 二维双QTL扫描
默认情况下，分析是通过EM算法的最大似然来执行的，(method="em")，verbose=FALSE阻止跟踪信息。
> out2d_em <- scantwo(calc, verbose=FALSE)
- 使用Haley-Knott回归方法进行2D扫描
> out2d_hk <- scantwo(calc, method="hk")
- 组合检验
最好是对观察到的数据使用组合检验，这样的结果将考虑表型分布、标记密度和缺失基因型数据，这一步需要耗费很长的时间。
> out2d_perm <- scantwo(calc, method="hk",n.perm=1000)
多性状双QTL扫描
因为数据量比较大，win下的R内存达不到，需要提交Linux，由于没有服务器的root权限，我自己技术有限，Linux下的Rstudio装不上，所以我把表型文件拆开，分别计算，分别保存为RData，再转回到win里进行后续分析。
这里我写了一个循环解决这个问题，其中1-20列为表型，21列为样本的ID，22列之后为基因型，如果和我一样有多个表型的，可以参照下面的循环更改。
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)
以说明书中的结果为例进行解释，详细内容如下：
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之间。
引用转载请注明出处，如有错误敬请指出。
