QTL-seq系列 | 使用easyQTLseq进行QTL-seq分析

作者：老王
链接：https://laowang2023.cn/2023/11/07/20231107-easyQTLseq/
来源：请叫我老王（个人博客）
本文采用 CC BY 4.0 许可协议 进行共享。转载请注明作者姓名和链接

QTL-seq是一种结合混池分离分析（Bulked segregant analysis，BSA）和高通量测序的方法^[1]，这种方法对于质量性状 QTL 或数量性状主效 QTL 定位具有省时省力省钱的特点。我们在前面已经介绍过QTL-seq分析的原理以及如何计算QTL-seq中的的置信区间。这里我们介绍如何使用R包 easyQTLseq 进行 QTL-seq 分析。

easyQTLseq 包的一个优势是可以灵活处理有无亲本信息、有一个还是两个亲本信息以及亲本是否有参考基因等各种情况，并且可以导出各个计算过程的结果以及生成高质量的图片。

安装easyQTLseq

easyQTLseq 包在 GitHub 上，使用 devtools 包安装。

# install devtools install.packages("devtools") # install easyQTLseq devtools::install_github("laowang1992/easyQTLseq")

准备输入数据

使用 GATK VariantsToTable

对亲本和混池的测序数据使用 GATK best practice pipeline 得到 VCF 文件，其中包含了各个样本的基因型（Genotype，GT）、基因型质量（Genotype Quality，GQ）、等位基因的深度（Allelic Depth，AD）以及一些其他的信息。为了加快数据读取和处理速度，使用 GATK 的 VariantsToTable 功能提取 CHROM、POS、REF、ALT 以及各个样本的 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

如果说你只有一个 VCF 文件，而且没有装 GATK，那如何做格式转换呢？easyQTLseq 包提供了一个函数 vcf2table() 用于格式转换。可以先用 vcfR 包中的 read.vcfR() 函数读取 VCF，然后用 vcf2table() 转换，函数输出可以直接用于 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)

数据来源

这里我们使用的数据来源于一篇甘蓝型油菜花色研究的文章^[2]（注意：样本命名和原文不同）。这篇文章使用两个 F₂ 群体将花色 QTL 定位到 A07 染色体的 21.96-30.14Mb 区间上，最终经过精细定位得到候选基因 BnaA07.PAP2 (BnaA07G0287000ZS) 的位置是 A07:26,522,430-26,524,093。

使用方法

加载包并导入数据

使用 VariantsToTable 准备好数据文件并加载 easyQTLseq 包后就可以导入数据了。

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位点

读取数据后，使用 select_sample_and_SNP() 函数对数据进行初步处理。选择亲本和混池并指定对应的材料信息，并指定群体类型和混池样本数量等信息，其他参数参考帮助信息。这个函数最终会返回一个包含所有信息的 QTLseq S3 对象。

x <- select_sample_and_SNP(data = data, highP = "qY", # 高表型亲本名称 lowP = "R3", # 低表型亲本名称 highB = "Y", # 极端高表型混池名称 lowB = "R", # 极端低表型混池名称 popType = "F2", # 群体类型，F2或者RIL。如果你的是DH群体也需要填RIL bulkSize = c(30, 30) # 混池大小，也就是高表型池和低表型池中个体的数量 ) 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"

除此之外，这个函数还可以处理各种不同的情况，例如上述的分离群体双亲都存在的情况、只有高表型亲本或低表型亲本存在的情况、没有亲本信息的情况、或者在 call SNP 时使用的是其中一个亲本的参考基因组。

# 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))

统计覆盖深度分布密度

在选择样本和 SNP 位点后，我们统计了每个样本所有 SNP 位点覆盖深度的分布，这样可以为后面根据覆盖深度过滤 SNP 位点提供参考。

depth_statistics(x = x, outPrefix = "outprefix" # "outprefix"是输出文件的前缀 )


覆盖深度分布密度

根据覆盖深度过滤SNP位点

低覆盖深度 SNP 的可靠性和准确性较低，而极高覆盖深度的 SNP 位点可能来自重复序列。这些 SNP 在进行后续分析前应该排除。

# default minimum coverage depth is 6, default maximum coverage depth is `average+3*sd`. x_filter <- filterDP(x = x)

统计SNP位点分布

在经过上述处理后，剩下的 SNP 可以用于 QTL-seq 分析，在分析之前应当检查一下 SNP 位点沿染色体的分布情况。这个分析所使用的参考基因组有很多没有挂载到染色体的 scaffold 片段，这里使用 targetChr 指定只分析19条染色体上的分布，并且这个基因组的染色体命名是按照 scaffoldA01、scaffoldA02 ... 这种格式，因此使用 chrLabel 重新指定染色体名称为 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是一个向量。你所用的基因组中可能有很多scaffold或contig，而这个向量中包含了需要展示的染色体名称，如果不指定这个参数，那么默认展示数据集中所有的染色体 chrLabel = chrLabel # 指定展示的染色体的名称，默认是数据集中染色体的名称，和"targetChr"参数对应。例如这个例子中数据集中第一个染色体名称是"scaffoldA01"，通过这个参数可以修改为"A01" )


SNP沿染色体分布

导出等位基因覆盖深度数据

如果想使用其他软件或方法来进行 QTL-seq 分析，就可以使用 export_dp() 函数导出等位基因覆盖深度信息，这个函数将生成名为 <outprefix>.Depth_information.txt|csv 的文件。

export_dp(x = x_filter, outPrefix = "outprefix")

滑窗统计计算 ΔSNP-index 和 ED 值

为了减小噪音使结果更加平滑，这里只用滑窗统计的方法来计算 SNP index、ΔSNP index、欧几里得距离（Euclidean Distance, ED）。在这一步，如果存在至少一个亲本的信息或者 call SNP 使用的是其中一个亲本的参考基因组，则 ΔSNP index 和 ED 值可以同时计算，如果没有亲本信息，则只计算 ED 值。同时这一步也会自动在当前目录导出滑窗统计的结果。

x_filter <- calc_index_etc(x = x_filter, outPrefix = "outprefix", winSize = 2000000, # 滑窗统计的窗口大小 winStep = 200000 # 滑窗统计的步长 )

导出图片

前面使用滑窗统计计算出结果，这一步则将上述计算结果按照染色体位置绘制图片并导出。

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)


delta SNP index

获取QTL信息

现在图表结果都有了，那么 QTL 位点的具体位置是哪里呢？如何得到QTL所在染色体的放大图片呢？这里可以使用 getQTL_and_exportFigure()来获得这些信息和结果。这个函数会将QTL在染色体上的位置区间、峰值所在位置、峰值 ΔSNP index 等信息导入名为 <outprefix>.99|95CI.csv 的表格中，并导出目标染色体的图片。当然和前面一样如果没有亲本信息，那么这一步就可以省略了。

getQTL_and_exportFigure(x = x_filter, outPrefix = "outprefix", minN = 20)


A07 染色体

到这里整个 QTL-seq 分析就结束了，在99%的置信区间下鉴定到一个 QTL 位于 A07:22,280,000-29,840,000，与文章结果基本相同，峰值所在位置为 25,740,000 bp，与候选基因 BnaA07.PAP2 的位置较为吻合。运行过程所产生的各种图表结果会自动导出，并且以 outPrefix 前缀命名。这个包的作者起名叫 easyQTLseq，主要是它可以自动处理有无亲本、亲本是否有参考基因组的各种情况，并且可以自动导出各种图表结果。但是很显然整个分析过程还是需要很多步骤才能完成，并且使用者至少需要了解R语言编程才能使用，并没有达到作者所说的 easy 的程度，因此后续可以结合 shiny 为这个扩展包写个图形界面，将各个参数设置好以后点击 run 可以直接运行所有步骤并生成所有结果，这种程度才算是真正的 easyQTLseq 。

引用

如果使用 easyQTLseq 进行 QTL-seq 分析，在撰写文章时可以在 Method 部分注明：

The QTL-seq analysis was performed using R package easyQTLseq (https://github.com/laowang1992/easyQTLseq.git).

参考文献

1. 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 ↩

2. 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 ↩