作者：老王
鏈接：https://laowang2023.cn/2023/11/07/20231107-easyQTLseq/
來(lái)源：請(qǐng)叫我老王（個(gè)人博客）
本文采用 CC BY 4.0 許可協(xié)議 進(jìn)行共享。轉(zhuǎn)載請(qǐng)注明作者姓名和鏈接
QTL-seq是一種結(jié)合混池分離分析（Bulked segregant analysis，BSA）和高通量測(cè)序的方法[1]，這種方法對(duì)于質(zhì)量性狀 QTL 或數(shù)量性狀主效 QTL 定位具有省時(shí)省力省錢(qián)的特點(diǎn)。我們?cè)谇懊嬉呀?jīng)介紹過(guò)QTL-seq分析的原理以及如何計(jì)算QTL-seq中的的置信區(qū)間。這里我們介紹如何使用R包 easyQTLseq 進(jìn)行 QTL-seq 分析。
easyQTLseq 包的一個(gè)優(yōu)勢(shì)是可以靈活處理有無(wú)親本信息、有一個(gè)還是兩個(gè)親本信息以及親本是否有參考基因等各種情況，并且可以導(dǎo)出各個(gè)計(jì)算過(guò)程的結(jié)果以及生成高質(zhì)量的圖片。
安裝easyQTLseq
easyQTLseq 包在 GitHub 上，使用 devtools 包安裝。
# install devtools
install.packages("devtools")
# install easyQTLseq
devtools::install_github("laowang1992/easyQTLseq")
準(zhǔn)備輸入數(shù)據(jù)
使用 GATK VariantsToTable
對(duì)親本和混池的測(cè)序數(shù)據(jù)使用 GATK best practice pipeline 得到 VCF 文件，其中包含了各個(gè)樣本的基因型（Genotype，GT）、基因型質(zhì)量（Genotype Quality，GQ）、等位基因的深度（Allelic Depth，AD）以及一些其他的信息。為了加快數(shù)據(jù)讀取和處理速度，使用 GATK 的 VariantsToTable 功能提取 CHROM、POS、REF、ALT 以及各個(gè)樣本的 GT、GQ 和 AD 信息。
java -Xmx30g -jar ${GATK} \
-R ${genome} -T VariantsToTable \
-F CHROM -F POS -F REF -F ALT -GF GT -GF AD -GF GQ \
-V ${BSA}.filter.SNPs.vcf.gz -o ${BSA}.filter.SNPs.table
使用 vcf2table
如果說(shuō)你只有一個(gè) VCF 文件，而且沒(méi)有裝 GATK，那如何做格式轉(zhuǎn)換呢？easyQTLseq 包提供了一個(gè)函數(shù) vcf2table() 用于格式轉(zhuǎn)換。可以先用 vcfR 包中的 read.vcfR() 函數(shù)讀取 VCF，然后用 vcf2table() 轉(zhuǎn)換，函數(shù)輸出可以直接用于 select_sample_and_SNP()。
library(vcfR)
library(easyQTLseq)
file_path <- system.file("extdata", "A07.SNPs.vcf.gz", package = "easyQTLseq")
x <- read.vcfR(file = file_path)
data <- vcf2table(x = x)
數(shù)據(jù)來(lái)源
這里我們使用的數(shù)據(jù)來(lái)源于一篇甘藍(lán)型油菜花色研究的文章[2]（注意：樣本命名和原文不同）。這篇文章使用兩個(gè) F2 群體將花色 QTL 定位到 A07 染色體的 21.96-30.14Mb 區(qū)間上，最終經(jīng)過(guò)精細(xì)定位得到候選基因 BnaA07.PAP2 (BnaA07G0287000ZS) 的位置是 A07:26,522,430-26,524,093。
使用方法
加載包并導(dǎo)入數(shù)據(jù)
使用 VariantsToTable 準(zhǔn)備好數(shù)據(jù)文件并加載 easyQTLseq 包后就可以導(dǎo)入數(shù)據(jù)了。
library(easyQTLseq)
file_path <- "BSA.filter.SNPs.table"
# readr::read_tsv() has a faster speed than read.table() when reading a file.
data <- readr::read_tsv(file = file_path)
## Rows: 2803830 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (11): CHROM, REF, ALT, R.GT, R.AD, R3.GT, R3.AD, Y.GT, Y.AD, qY.GT, qY.AD
## dbl (5): POS, R.GQ, R3.GQ, Y.GQ, qY.GQ
##
## ? Use `spec()` to retrieve the full column specification for this data.
## ? Specify the column types or set `show_col_types = FALSE` to quiet this message.
指定樣本信息并篩選SNP位點(diǎn)
讀取數(shù)據(jù)后，使用 select_sample_and_SNP() 函數(shù)對(duì)數(shù)據(jù)進(jìn)行初步處理。選擇親本和混池并指定對(duì)應(yīng)的材料信息，并指定群體類(lèi)型和混池樣本數(shù)量等信息，其他參數(shù)參考幫助信息。這個(gè)函數(shù)最終會(huì)返回一個(gè)包含所有信息的 QTLseq S3 對(duì)象。
x <- select_sample_and_SNP(data = data,
highP = "qY", # 高表型親本名稱(chēng)
lowP = "R3", # 低表型親本名稱(chēng)
highB = "Y", # 極端高表型混池名稱(chēng)
lowB = "R", # 極端低表型混池名稱(chēng)
popType = "F2", # 群體類(lèi)型，F(xiàn)2或者RIL。如果你的是DH群體也需要填RIL
bulkSize = c(30, 30) # 混池大小，也就是高表型池和低表型池中個(gè)體的數(shù)量
)
x
## $data
## # A tibble: 722,424 × 12
## CHROM POS REF ALT HP.DP LP.DP HB.HP.AD HB.LP.AD HB.DP LB.HP.AD
## <chr> <dbl> <chr> <chr> <int> <int> <int> <int> <int> <int>
## 1 scaffoldA01 275959 A C 17 24 6 6 12 16
## 2 scaffoldA01 320766 G A 19 11 6 6 12 10
## 3 scaffoldA01 361230 A T 18 6 6 5 11 8
## 4 scaffoldA01 361486 G T 26 9 13 5 18 8
## 5 scaffoldA01 361875 A C 19 5 8 0 8 4
## 6 scaffoldA01 361884 T G 19 5 8 0 8 6
## 7 scaffoldA01 362128 A C 19 6 7 5 12 4
## 8 scaffoldA01 362153 T G 20 5 8 5 13 6
## 9 scaffoldA01 364170 C T 17 14 12 9 21 14
## 10 scaffoldA01 365454 T C 18 10 11 5 16 15
## # ? 722,414 more rows
## # ? 2 more variables: LB.LP.AD <int>, LB.DP <int>
##
## $highP
## [1] "qY"
##
## $lowP
## [1] "R3"
##
## $highB
## [1] "Y"
##
## $lowB
## [1] "R"
##
## $popType
## [1] "F2"
##
## $bulkSize
## [1] 30 30
##
## $slidwin
## data frame with 0 columns and 0 rows
##
## $chrLen
## # A tibble: 523 × 2
## CHROM Len
## <chr> <dbl>
## 1 scaffold0022 70801
## 2 scaffold0025 4039467
## 3 scaffold0026 8552783
## 4 scaffold0027 5993037
## 5 scaffold0032 11673
## 6 scaffold0034 36998
## 7 scaffold0036 262
## 8 scaffold0037 23809
## 9 scaffold0038 32054
## 10 scaffold0039 6403
## ? 513 more rows
## ? Use `print(n = ...)` to see more rows
##
## attr(,"class")
## [1] "QTLseq" "WithParent" "BothParent"
除此之外，這個(gè)函數(shù)還可以處理各種不同的情況，例如上述的分離群體雙親都存在的情況、只有高表型親本或低表型親本存在的情況、沒(méi)有親本信息的情況、或者在 call SNP 時(shí)使用的是其中一個(gè)親本的參考基因組。
# If only one parent is present, e.g. high parent.
x_onlyHP <- select_sample_and_SNP(data = data, highP = "qY", highB = "Y", lowB = "R", popType = "F2", bulkSize = c(30, 30))
# If no parent is present.
x_noParent <- select_sample_and_SNP(data = data, highB = "Y", lowB = "R", popType = "F2", bulkSize = c(30, 30))
# If no parent is present, but high parent has a reference genome, this reference genome is used for SNP calling. Then the `highP` parameter should be "REF".
x_HPisREF <- select_sample_and_SNP(data = data, highP = "REF", highB = "Y", lowB = "R", popType = "F2", bulkSize = c(30, 30))
統(tǒng)計(jì)覆蓋深度分布密度
在選擇樣本和 SNP 位點(diǎn)后，我們統(tǒng)計(jì)了每個(gè)樣本所有 SNP 位點(diǎn)覆蓋深度的分布，這樣可以為后面根據(jù)覆蓋深度過(guò)濾 SNP 位點(diǎn)提供參考。
depth_statistics(x = x,
outPrefix = "outprefix" # "outprefix"是輸出文件的前綴
)
根據(jù)覆蓋深度過(guò)濾SNP位點(diǎn)
低覆蓋深度 SNP 的可靠性和準(zhǔn)確性較低，而極高覆蓋深度的 SNP 位點(diǎn)可能來(lái)自重復(fù)序列。這些 SNP 在進(jìn)行后續(xù)分析前應(yīng)該排除。
# default minimum coverage depth is 6, default maximum coverage depth is `average+3*sd`.
x_filter <- filterDP(x = x)
統(tǒng)計(jì)SNP位點(diǎn)分布
在經(jīng)過(guò)上述處理后，剩下的 SNP 可以用于 QTL-seq 分析，在分析之前應(yīng)當(dāng)檢查一下 SNP 位點(diǎn)沿染色體的分布情況。這個(gè)分析所使用的參考基因組有很多沒(méi)有掛載到染色體的 scaffold 片段，這里使用 targetChr 指定只分析19條染色體上的分布，并且這個(gè)基因組的染色體命名是按照 scaffoldA01、scaffoldA02 ... 這種格式，因此使用 chrLabel 重新指定染色體名稱(chēng)為 A01、A02 ... 。
targetChr <- c(paste0("scaffoldA", sprintf("%02d", 1:10)), paste0("scaffoldC", sprintf("%02d", 1:9)))
chrLabel <- c(paste0("A", sprintf("%02d", 1:10)), paste0("C", sprintf("%02d", 1:9)))
SNP_distribution(x = x_filter, outPrefix = "outprefix",
targetChr = targetChr, # 這里的targetChr是一個(gè)向量。你所用的基因組中可能有很多scaffold或contig，而這個(gè)向量中包含了需要展示的染色體名稱(chēng)，如果不指定這個(gè)參數(shù)，那么默認(rèn)展示數(shù)據(jù)集中所有的染色體
chrLabel = chrLabel # 指定展示的染色體的名稱(chēng)，默認(rèn)是數(shù)據(jù)集中染色體的名稱(chēng)，和"targetChr"參數(shù)對(duì)應(yīng)。例如這個(gè)例子中數(shù)據(jù)集中第一個(gè)染色體名稱(chēng)是"scaffoldA01"，通過(guò)這個(gè)參數(shù)可以修改為"A01"
)
導(dǎo)出等位基因覆蓋深度數(shù)據(jù)
如果想使用其他軟件或方法來(lái)進(jìn)行 QTL-seq 分析，就可以使用 export_dp() 函數(shù)導(dǎo)出等位基因覆蓋深度信息，這個(gè)函數(shù)將生成名為 <outprefix>.Depth_information.txt|csv 的文件。
export_dp(x = x_filter, outPrefix = "outprefix")
滑窗統(tǒng)計(jì)計(jì)算 ΔSNP-index 和 ED 值
為了減小噪音使結(jié)果更加平滑，這里只用滑窗統(tǒng)計(jì)的方法來(lái)計(jì)算 SNP index、ΔSNP index、歐幾里得距離（Euclidean Distance, ED）。在這一步，如果存在至少一個(gè)親本的信息或者 call SNP 使用的是其中一個(gè)親本的參考基因組，則 ΔSNP index 和 ED 值可以同時(shí)計(jì)算，如果沒(méi)有親本信息，則只計(jì)算 ED 值。同時(shí)這一步也會(huì)自動(dòng)在當(dāng)前目錄導(dǎo)出滑窗統(tǒng)計(jì)的結(jié)果。
x_filter <- calc_index_etc(x = x_filter, outPrefix = "outprefix",
winSize = 2000000, # 滑窗統(tǒng)計(jì)的窗口大小
winStep = 200000 # 滑窗統(tǒng)計(jì)的步長(zhǎng)
)
導(dǎo)出圖片
前面使用滑窗統(tǒng)計(jì)計(jì)算出結(jié)果，這一步則將上述計(jì)算結(jié)果按照染色體位置繪制圖片并導(dǎo)出。
export_figure(x = x_filter,
outPrefix = "outprefix",
targetChr = targetChr, # Target chromosome to be drawn in figures, default is all chromosomes in the data.
chrLabel = chrLabel, # The label for chromosome shown in figures, default is chromosome names in the data.
minN = 20, # Too few SNPs in a window will result in noise, the windows containing SNPs less than minN will be omitted in figures.
width = 12, height = 4)
獲取QTL信息
現(xiàn)在圖表結(jié)果都有了，那么 QTL 位點(diǎn)的具體位置是哪里呢？如何得到QTL所在染色體的放大圖片呢？這里可以使用 getQTL_and_exportFigure()來(lái)獲得這些信息和結(jié)果。這個(gè)函數(shù)會(huì)將QTL在染色體上的位置區(qū)間、峰值所在位置、峰值 ΔSNP index 等信息導(dǎo)入名為 <outprefix>.99|95CI.csv 的表格中，并導(dǎo)出目標(biāo)染色體的圖片。當(dāng)然和前面一樣如果沒(méi)有親本信息，那么這一步就可以省略了。
getQTL_and_exportFigure(x = x_filter, outPrefix = "outprefix", minN = 20)
到這里整個(gè) QTL-seq 分析就結(jié)束了，在99%的置信區(qū)間下鑒定到一個(gè) QTL 位于 A07:22,280,000-29,840,000，與文章結(jié)果基本相同，峰值所在位置為 25,740,000 bp，與候選基因 BnaA07.PAP2 的位置較為吻合。運(yùn)行過(guò)程所產(chǎn)生的各種圖表結(jié)果會(huì)自動(dòng)導(dǎo)出，并且以 outPrefix 前綴命名。這個(gè)包的作者起名叫 easyQTLseq，主要是它可以自動(dòng)處理有無(wú)親本、親本是否有參考基因組的各種情況，并且可以自動(dòng)導(dǎo)出各種圖表結(jié)果。但是很顯然整個(gè)分析過(guò)程還是需要很多步驟才能完成，并且使用者至少需要了解R語(yǔ)言編程才能使用，并沒(méi)有達(dá)到作者所說(shuō)的 easy 的程度，因此后續(xù)可以結(jié)合 shiny 為這個(gè)擴(kuò)展包寫(xiě)個(gè)圖形界面，將各個(gè)參數(shù)設(shè)置好以后點(diǎn)擊 run 可以直接運(yùn)行所有步驟并生成所有結(jié)果，這種程度才算是真正的 easyQTLseq 。
引用
如果使用 easyQTLseq 進(jìn)行 QTL-seq 分析，在撰寫(xiě)文章時(shí)可以在 Method 部分注明：
The QTL-seq analysis was performed using R package easyQTLseq (https://github.com/laowang1992/easyQTLseq.git).
參考文獻(xiàn)
-
Takagi H, Abe A, Yoshida K, et al. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J. 2013;74(1):174-183. doi:10.1111/tpj.12105 ?
-
Ye S, Hua S, Ma T, et al. Genetic and multi-omics analyses reveal BnaA07.PAP2In-184-317 as the key gene conferring anthocyanin-based color in Brassica napus flowers. J Exp Bot. 2022;73(19):6630-6645. doi:10.1093/jxb/erac312 ?
