目錄：
1.?Data Input
2. Calculating Genotype Probabilities
3.?Performing a genome scan
4.?Performing a permutation test
5.?
ok徙赢，進(jìn)入正文唱较。
1. Data Input
和 R/qtl 要求的數(shù)據(jù)格式不同，R/qtl2 要求每一個(gè)數(shù)據(jù)為一個(gè)單獨(dú)的文件扳炬，但可以最后總結(jié)到一個(gè) YAML 或 JSON 格式的文本里，一次性讀入堵幽，還是挺方便的捕透。
首先看一下這個(gè)?YAML 或 JSON 格式的文本文件菇爪，指定了群體類型，每一個(gè)文件的名字袋励，以及基因型信息等侥啤。
其中当叭，geno、pheno盖灸、gmap 三個(gè)文件是必須的蚁鳖，其它文件是可選的。
讀入使用 read_cross2 函數(shù)赁炎。
read_cross2(file, quiet=TRUE)
# file 指定 YAML 文件的路徑位置即可醉箕，quiet 指定是否輸出讀入過程。
以下分別是這 3 個(gè)文件的格式示例：
2. Calculating Genotype Probabilities
類似 qtl 包中的分析方法徙垫，QTL 分析首先需要計(jì)算已有標(biāo)記之間的基因型概率讥裤。
使用的函數(shù)為 calc_genoprob()?。不過在使用之前姻报，需要先在 markers 之間插入?pseudomarkers己英，以確定要計(jì)算基因型可能性的位置。在遺傳圖譜中插入 pseudomarkers 的函數(shù)為?insert_pseudomarkers() 吴旋。
舉例示范一下：
library(qtl2)
## 讀入包里自帶的示例文件
iron <- read_cross2(file = system.file("extdata", "iron.zip", package="qtl2") )
## 在遺傳圖譜內(nèi)插入?pseudomarkers
map <- insert_pseudomarkers(map=iron$gmap, step=1)
# step 指的是?pseudomarkers 和標(biāo)記之間的舉例為 1 cM
## 計(jì)算 pseudomarkers?的基因型概率
pr <- calc_genoprob(cross=iron, map=map, error_prob=0.002, cores=4)
# error_prob 指定假設(shè)的基因分型錯(cuò)誤概率损肛，默認(rèn)為0.0001；cores 是在集群上指定運(yùn)算使用的內(nèi)核數(shù)荣瑟，parallel::detectCores() 查看可用的內(nèi)核數(shù)
這個(gè)過程就到此為止了治拿，接下來通過幾張圖看一下每一步的效果。
summary(iron)
head(iron$gmap)
map <- insert_pseudomarkers(map=iron$gmap, step=1)
head(map, n=1)? ? ????????# 查看1號染色體上的標(biāo)記
pr <- calc_genoprob(cross=iron, map=map, error_prob=0.002)
(pr$`19`)[1:3,,"c19.loc4"]? ?????# 查看19號染色體上該標(biāo)記處前3個(gè)個(gè)體的基因型概率
此外笆焰，如果需要使用加性模型進(jìn)行基因組掃描忍啤，需要將基因型概率轉(zhuǎn)換為等位基因概率。
apr <- genoprob_to_alleleprob(probs=pr)
3.?Performing a genome scan
原理就不在這塊詳細(xì)說了仙辟，有興趣的可以去看?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
雖然可以看到部分染色體上存在峰，但是否顯著戴尸？對結(jié)果如何評估粟焊？
類似于 qtl 包中的操作，需要進(jìn)行排列檢驗(yàn)孙蒙，打亂基因型和表型项棠，去除二者之間的關(guān)聯(lián)，然后進(jìn)行計(jì)算全基因組的最大 LOD 評分作為閾值挎峦，從而確定結(jié)果是否是隨機(jī)出現(xiàn)的香追。需要用到的函數(shù)為 scan1perm() 。
operm <- scan1perm(genoprobs = pr, pheno = iron$pheno, Xcovar = Xcovar, n_perm = 1000)
#?n_perm 指定排列檢驗(yàn)的重復(fù)次數(shù)
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 染色體單獨(dú)計(jì)算閾值坦胶，可使用下面的命令：
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
這一步就要根據(jù)排列檢驗(yàn)得到的閾值透典，以及計(jì)算結(jié)果晴楔，確定目標(biāo)的 QTL 區(qū)間了。首先要確定 LOD peaks 的位置峭咒，再計(jì)算其可靠的區(qū)間范圍税弃。
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 確定貝葉斯置信區(qū)間，0.95 指的是有 95% 的可能性落在該區(qū)間
#?expand2markers 指的是是否擴(kuò)展至相鄰標(biāo)記
該函數(shù)也可以確定一條染色體上的多個(gè)峰凑队，不過需要 peakdrop?指定兩個(gè)相鄰峰值之間的最低值必須下降的量则果。
find_peaks(scan1_output = out, map = map, threshold = thr, peakdrop = 1.8, prob = 0.95, expand2markers = FALSE)
