目录：

1. Data Input

2. Calculating Genotype Probabilities

3. Performing a genome scan

4. Performing a permutation test

5. 

---

ok，进入正文。

1. Data Input

和 R/qtl 要求的数据格式不同，R/qtl2 要求每一个数据为一个单独的文件，但可以最后总结到一个 YAML 或 JSON 格式的文本里，一次性读入，还是挺方便的。

首先看一下这个 YAML 或 JSON 格式的文本文件，指定了群体类型，每一个文件的名字，以及基因型信息等。

其中，geno、pheno、gmap 三个文件是必须的，其它文件是可选的。

读入使用 read_cross2 函数。

read_cross2(file, quiet=TRUE)

# file 指定 YAML 文件的路径位置即可，quiet 指定是否输出读入过程。

以下分别是这 3 个文件的格式示例：

geno pheno gmap yaml

2. Calculating Genotype Probabilities

类似 qtl 包中的分析方法，QTL 分析首先需要计算已有标记之间的基因型概率。

使用的函数为 calc_genoprob() 。不过在使用之前，需要先在 markers 之间插入 pseudomarkers，以确定要计算基因型可能性的位置。在遗传图谱中插入 pseudomarkers 的函数为 insert_pseudomarkers() 。

举例示范一下：

library(qtl2)

## 读入包里自带的示例文件

iron <- read_cross2(file = system.file("extdata", "iron.zip", package="qtl2") )

## 在遗传图谱内插入 pseudomarkers

map <- insert_pseudomarkers(map=iron$gmap, step=1)

# step 指的是 pseudomarkers 和标记之间的举例为 1 cM

## 计算 pseudomarkers 的基因型概率

pr <- calc_genoprob(cross=iron, map=map, error_prob=0.002, cores=4)

# error_prob 指定假设的基因分型错误概率，默认为0.0001；cores 是在集群上指定运算使用的内核数，parallel::detectCores() 查看可用的内核数

这个过程就到此为止了，接下来通过几张图看一下每一步的效果。

summary(iron)

head(iron$gmap)

map <- insert_pseudomarkers(map=iron$gmap, step=1)

head(map, n=1)            # 查看1号染色体上的标记

pr <- calc_genoprob(cross=iron, map=map, error_prob=0.002)

(pr$`19`)[1:3,,"c19.loc4"]       # 查看19号染色体上该标记处前3个个体的基因型概率

此外，如果需要使用加性模型进行基因组扫描，需要将基因型概率转换为等位基因概率。

apr <- genoprob_to_alleleprob(probs=pr)

3. Performing a genome scan

原理就不在这块详细说了，有兴趣的可以去看 R/qtl 定位分析（三）Single-QTL analysis 。

以 Haley-Knott regression 方法为例：

Xcovar <- get_x_covar(iron) 

out <- scan1(genoprobs = pr, pheno = iron$pheno, Xcovar=Xcovar, cores=4)

输出是 LOD 得分矩阵，位置×表型。

head(out, n=10)

plot_scan1(out, map = map, lodcolumn = "liver")

4. Performing a permutation test

虽然可以看到部分染色体上存在峰，但是否显著？对结果如何评估？

类似于 qtl 包中的操作，需要进行排列检验，打乱基因型和表型，去除二者之间的关联，然后进行计算全基因组的最大 LOD 评分作为阈值，从而确定结果是否是随机出现的。需要用到的函数为 scan1perm() 。

operm <- scan1perm(genoprobs = pr, pheno = iron$pheno, Xcovar = Xcovar, n_perm = 1000)

# n_perm 指定排列检验的重复次数

hist(operm[,'liver'], breaks = 50, xlab = "LOD", main = "LOD scores for liver scan with threshold in red")abline(v = summary(operm)[,'liver'], col = 'red', lwd = 2)

summary(operm, alpha=c(0.1, 0.05))

如果要对 X 染色体单独计算阈值，可使用下面的命令：

operm2 <- scan1perm(pr, iron$pheno, Xcovar=Xcovar, n_perm=1000, perm_Xsp=TRUE, chr_lengths=chr_lengths(map))

summary(operm2, alpha=c(0.2, 0.05))

5. Finding LOD peaks

这一步就要根据排列检验得到的阈值，以及计算结果，确定目标的 QTL 区间了。首先要确定 LOD peaks 的位置，再计算其可靠的区间范围。

operm <- scan1perm(pr, iron$pheno, Xcovar=Xcovar, n_perm=1000)

thr <- summary(operm)

find_peaks(scan1_output = out, map = map, threshold = thr, prob = 0.95, expand2markers = FALSE)

# prob 确定贝叶斯置信区间，0.95 指的是有 95% 的可能性落在该区间

# expand2markers 指的是是否扩展至相邻标记

该函数也可以确定一条染色体上的多个峰，不过需要 peakdrop 指定两个相邻峰值之间的最低值必须下降的量。

find_peaks(scan1_output = out, map = map, threshold = thr, peakdrop = 1.8, prob = 0.95, expand2markers = FALSE)

---