2.QTL定位：Rqtl —— Two-QTL scans

2022-01-09 本文已影响0人 Wei_Sun

在上篇帖子中介绍了如何用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之间。

引用转载请注明出处，如有错误敬请指出。