前言
我在犹豫是否补充一篇怎样拆分染色体进行比对和变异检测（主要是为了提高比对和变异检测效率）的帖子。开始没有规划这篇内容的原因有两点：一是内容比较简单，觉得可能用处不大；二是个人精力也是十分有限（也就是懒）。
后来想想写这个系列的初衷主要的是为了给一些刚入门需要做BSA分析的同学提供一些比较好的参考资料，貌似补充这样的内容又有些必要了。不管怎样，这个系列的主要目标就是大家跟着这个系列做一遍至少能自主完成一遍BSA分析，个人觉得中间可能会用到的常用脚本放在实用脚本分享的内容里了，有需要可以自取。希望能为大家减少一些学习过程中的磕绊吧！
如果觉得很需要拆分染色体进行比对和变异检测这篇介绍的话，可以后台私信或者评论留言，有需要我再写吧。插播一句，如果有些别的想看的系列，也可以留言，会将诉求多的系列靠前安排。后续还会继续推出GWAS，基因组组装注释，转录组，泛基因组，机器学习等多系列的分析流程，敬请期待吧。
正文
我们从上一篇推文的分析结果中拿到了VCF文件，在进一步分析的时候要对数据进行过滤整合。由于QTLseqr只支持两个样本的数据输入，当把亲本1（P1），亲本2（P2），子代1（S1），子代2 (S2) 四组样本数据合并成一个VCF的进行分析的时候就会报错。而我们还需要亲本数据进行位点的筛选，因此，我们可以先按照亲本纯合且差异的标准对位点进行过滤得到merge.filter.vcf（过滤脚本会在实用脚本分享里分享），再从中将子代两个池的样本SNP位点数据提取出来得到用于分析的input.snp.vcf。这里要用到软件VCFtools。
vcftools安装
# 普通下载
git clone https://github.com/vcftools/vcftools.git
tar zxvf ./vcftools-0.1.16.tar.gz
cd vcftools-0.1.16/
./configure --prefix=`pwd`
make && make install
#导入环境变量
echo "PATH=$PATH:~/softwarepath/vcftools-0.1.16/bin" >>~/.bashrc
source ~/.bashrc
# 懒人安装
conda create -n vcftools
conda activate vcftools
conda install -c bioconda vcftools=0.1.16 -y
#测试是否安装成功
vcftools
vcftools使用
这里主要用到vcftools中提取多个子样本变异位点信息以及分离SNP、InDel位点的功能。
#提取子样本SNP变异位点，该软件也支持压缩的vcf输入，--gzvcf vcf.gz即可。
vcftools --vcf merge.filter.vcf \
--recode --recode-INFO-all --stdout \
--indv S1
--indv S2 --remove-indels > c
##参数解读
--recode 在结果文件中显示等位基因频率，不能少，否则没有结果文件
--recode-INFO-all 保留原vcf文件INFO列的所有内容
--stdout 使文件可以重定向，可以配合管道符压缩文件
--indv 保留样本名
（这里还可以使用--remove-indv P1 --remove-indv P2，是同样的效果）
--remove-indels 保留SNP位点
（--keep-only-indels 保留InDel位点）
## 扩展
在重测序项目中，当样品较多时，也可将保留样品的样本名保存到文件Keep.txt中，格式如下：
S1
S2
然后使用--keep Keep.txt 保留样品
或者将保留样品的样本名保存到文件remove.txt中，格式如下：
P1
P2
然后使用--remove remove.txt 删除样品
QTLseqr
QTLseqr是基于R语言开发的一款用于BSA分析的R包。包含G's value 和Delta SNP-index两种方法。
详细信息，可以参考官网的介绍（https://github.com/bmansfeld/QTLseqR/），官网还有英文版说明文档，对算法都有具体的介绍。
安装
install.packages("devtools") #安装devtools
library(devtools)
install_github("bmansfeld/QTLseqr") #使用devtools安装QTLseqr
library(QTLseqr)
#还需要搭配ggplot使用
install.packages("ggplot2")
library(ggplot2)
#library后无报错信息，则说明安装成功
使用
QTLseqr包的输入文件不支持vcf格式的文件，所以需要对VCF文件格式进行转换。我们这里又用到了之前介绍的软件GATK。
本文使用“Mapping of Quantitative Trait Loci Underlying Cold Tolerance in Rice Seedlings via High- Throughput Sequencing of Pooled Extremes”文章中的数据作为测试数据。
gatk VariantsToTable -R ref.fa -V input.snp.vcf \
-F CHROM -F POS -F REF -F ALT \
-GF AD -GF DP -GF GQ -GF PL \
-O input.table
#参数介绍，所需是分析必须包括的信息，推荐可加可不加
ref.fa 参考基因组
input.snp.vcf vcf文件
所需的VCF字段（-F）是chrom（染色体）、POS（位置）。
所需的基因型字段（-GF）为AD（等位基因深度）、DP（深度）。
推荐的字段是REF（参考等位基因）和ALT（替代等位基因)
推荐的基因型字段是PL（phred比例可能性)和GQ(基因型质量）。
得到整理后的VCF，格式如下：
本文使用“Mapping of Quantitative Trait Loci Underlying Cold Tolerance in Rice Seedlings via High- Throughput Sequencing of Pooled Extremes”文章中的数据作为测试数据。数据获得方式如下：
#download and load data package (~50Mb)
devtools::install_github("bmansfeld/Yang2013data")
library("Yang2013data")
#Import the data
rawData <- system.file(
"extdata",
"Yang_et_al_2013.table",
package = "Yang2013data",
mustWork = TRUE)
library(QTLseqr)
library(ggplot2)
# 如果使用自己的数据，可以直接
rawData <- "input.table"
#示例数据中子代数据的名称是SRR834931，SRR834927，这个可以根据实际需求去更改，但是要与table表头命名一致，最好不要以N或数字开头。
HighBulk <- "SRR834931" #高值混池名
LowBulk <- "SRR834927" #低值混池名
#染色体号列表，这里是Chr1, Chr2 . . . Chr12，实际分析以ref.fa的编号为主。
Chroms <- paste0(rep("Chr", 12), 1:12)
#导入数据
df <- importFromGATK(file = rawData,
highBulk = HighBulk,
lowBulk = LowBulk,
chromList = Chroms)
#查看数据
head(df)
## CHROM POS REF ALT AD_REF.LOW AD_ALT.LOW DP.LOW GQ.LOW PL.LOW
## 1 Chr1 31071 A G 34 36 70 99 897,0,855
## 2 Chr1 31478 C T 34 52 86 99 1363,0,844
## 3 Chr1 33667 A G 20 48 68 99 1331,0,438
## 4 Chr1 34057 C T 38 40 78 99 1059,0,996
## 5 Chr1 35239 A C 25 36 61 99 987,0,630
## 6 Chr1 38389 T C 36 42 78 99 1066,0,906
## SNPindex.LOW AD_REF.HIGH AD_ALT.HIGH DP.HIGH GQ.HIGH PL.HIGH
## 1 0.5142857 26 22 48 99 522,0,698
## 2 0.6046512 40 34 74 99 848,0,1099
## 3 0.7058824 24 29 53 99 765,0,599
## 4 0.5128205 29 26 55 99 673,0,772
## 5 0.5901639 40 60 100 99 1632,0,1015
## 6 0.5384615 42 40 82 99 984,0,1105
## SNPindex.HIGH REF_FRQ deltaSNP
## 1 0.4583333 0.5084746 -0.055952381
## 2 0.4594595 0.4625000 -0.145191703
## 3 0.5471698 0.3636364 -0.158712542
## 4 0.4727273 0.5037594 -0.040093240
## 5 0.6000000 0.4037267 0.009836066
## 6 0.4878049 0.4875000 -0.050656660
# 结果信息
• CHROM - The chromosome this SNP is in
• POS - The position on the chromosome in nt
• REF - The reference allele at that position
• ALT - The alternate allele
• DP.HIGH - 高值池reads深度
• AD_REF.HIGH - 高值池中与参考基因组一致的基因型深度
• AD_ALT.HIGH - 高值池中与变异一致的基因型深度
• GQ.HIGH - 高值池在这个位点的基因型质量值
• SNPindex.HIGH - 高值池SNP-index，计算方式：AD_ALT.HIGH/DP.HIGH
• *.lOW - 和高值池一致
• REF_FRQ - The reference allele frequency as defined above
• deltaSNP - ∆(SNP-index) ，计算方式：SNPindex.HIGH - SNPindex.LOW
# 在上面的结果就可以看到，每个位点的index以及deltaSNPindex已经被计算出来了。导出计算结果,主要方便自己查看位点信息
write.table(df,"delta_index_caculate_result.xls",sep="\t",quote = F)
除了对亲本纯合且差异的位点做挑选以外，我们还需要结合数据具体的深度，等位基因频率等特征值的分布情况进行过滤。因此，我们首先对这些值的分布情况进行可视化。
深度分布图
# 深度分布
depth_plot < ggplot(data = df) +
geom_histogram(aes(x = DP.HIGH + DP.LOW)) +
xlim(0,1000) #结合实际情况进行调整
#输出图片
pdf(file="SNP_depth.pdf")
depth_plot
dev.off()
深度分布图
如上图，我们可以看到数据深度大部分分布在0-375之间,所以这个值定在375-400就差不多了。
等位基因频率分布图
# 等位基因频率
alle_plot <- ggplot(data = df) +
geom_histogram(aes(x = REF_FRQ))
#输出图片
pdf(file="ref_alle_frequency.pdf")
alle_plot
dev.off()
等位基因频率分布图
如上图，同理，我们可以看得到大部分的位点等位基因频率都是在0.2~0.8之间，有少部分异常的位点，后续可以进行过滤。
高低池SNP-index分布图
接下来我们可以查看一下高值、低值池的SNP-index的分布图，以此来检查数据是否存在异常。
#以高值池SNP-index分布，低值池类似
high_index_plot <- ggplot(data = df) +
geom_histogram(aes(x = SNPindex.HIGH))
#输出图片
pdf(file="high_index.pdf")
high_index_plot
dev.off()
在F2群体中，大多数SNP应该近似正态分布分布在0.5左右，极端值（0，1）是我们正期望看到的位点所在。如果不符合这个分布规律，需要进行问题排查，回去检查样本是否正确，数据是否有问题等等。
数据过滤
df_filt <-filterSNPs(
SNPset = df,
refAlleleFreq = 0.20, #等位基因频率过滤，0.2表示0.2~0.8之间
minTotalDepth = 100, #最小过滤深度
maxTotalDepth = 400, #最大过滤深度
minSampleDepth = 40, #单个样本过滤深度
minGQ = 99, ### genotype quality 过滤
verbose = TRUE
)
# 输出日志信息
## Filtering by reference allele frequency: 0.2 <= REF_FRQ <= 0.8
## ...Filtered 59754 SNPs
## Filtering by total sample read depth: Total DP >= 100
## ...Filtered 378814 SNPs
## Filtering by total sample read depth: Total DP <= 400
## ...Filtered 744 SNPs
## Filtering by per sample read depth: DP >= 40
## ...Filtered 4729 SNPs
## Filtering by Genotype Quality: GQ >= 99
## ...Filtered 6677 SNPs
## Filtering by difference between bulks <= 100
## ...Filtered 1418 SNPs
## Original SNP number: 969487, Filtered: 452136, Remaining: 517351
！！！最后一行就可以看到总的SNP个数为969487，过滤SNP数目为：452136，还剩下SNP数目为517351。
过滤完就可以进行我们最为关键的一步了。下面分析分为两种方法，一个是G'，一个是∆(SNP-index)。
runGprimeAnalysis 函数
G' analysis
# 滑窗计算
df_out <- runGprimeAnalysis(df_filt,
windowSize = 1e6, #窗口大小，根据基因组大小进行适当调整
outlierFilter = "deltaSNP", #两种过滤方式，Hampel和deltaSNP
filterThreshold = 0.1 #阈值线，相当于99%的置信区间
)
# 数据格式
head(df_filt)
## CHROM POS REF ALT AD_REF.LOW AD_ALT.LOW DP.LOW GQ.LOW PL.LOW
## 1 Chr1 31071 A G 34 36 70 99 897,0,855
## 2 Chr1 31478 C T 34 52 86 99 1363,0,844
## 3 Chr1 33667 A G 20 48 68 99 1331,0,438
## 4 Chr1 34057 C T 38 40 78 99 1059,0,996
## 5 Chr1 35239 A C 25 36 61 99 987,0,630
## 6 Chr1 38389 T C 36 42 78 99 1066,0,906
## SNPindex.LOW AD_REF.HIGH AD_ALT.HIGH DP.HIGH GQ.HIGH PL.HIGH
## 1 0.5142857 26 22 48 99 522,0,698
## 2 0.6046512 40 34 74 99 848,0,1099
## 3 0.7058824 24 29 53 99 765,0,599
## 4 0.5128205 29 26 55 99 673,0,772
## 5 0.5901639 40 60 100 99 1632,0,1015
## 6 0.5384615 42 40 82 99 984,0,1105
## SNPindex.HIGH REF_FRQ deltaSNP nSNPs tricubeDeltaSNP minDP
## 1 0.4583333 0.5084746 -0.055952381 876 -0.07376705 48
## 2 0.4594595 0.4625000 -0.145191703 876 -0.07376955 74
## 3 0.5471698 0.3636364 -0.158712542 878 -0.07378303 53
## 4 0.4727273 0.5037594 -0.040093240 878 -0.07378543 55
## 5 0.6000000 0.4037267 0.009836066 878 -0.07379271 61
## 6 0.4878049 0.4875000 -0.050656660 878 -0.07381211 78
## tricubeDP CI_95 CI_99 G Gprime pvalue
## 1 70 -0.1714286 -0.2285714 0.35697092 1.922409 0.4658957
## 2 70 -0.1714286 -0.2285714 3.38161778 1.922473 0.4658766
## 3 70 -0.1714286 -0.2285714 3.23692853 1.922813 0.4657738
## 4 70 -0.1714286 -0.2285714 0.20748641 1.922873 0.4657555
## 5 70 -0.1714286 -0.2285714 0.01521766 1.923057 0.4657000
## 6 70 -0.1714286 -0.2285714 0.41076961 1.923547 0.4655521
## negLog10Pval qvalue
## 1 0.3317113 0.6681496
## 2 0.3317291 0.6681338
## 3 0.3318249 0.6680525
12
## 4 0.3318420 0.6680330
## 5 0.3318938 0.6679837
## 6 0.3320317 0.6678859
# 格式解读
• nSNPs - the number of SNPs bracketing the focal SNP within the set sliding window
• tricubeDeltaSNP - the tricube-smoothed ∆(SNP-index)
• G - the G value for the SNP
• Gprime - the tricube-smoothed G value
• pvalue - the p-value calculated by non-parametric estimation
• negLog10Pval - the −log10(p-value)
• qvalue - Benjamini-Hochberg adjusted p-values
#绘制结果图
g_plot <- plotQTLStats(SNPset = df_filt, var = "Gprime", plotThreshold = TRUE, q = 0.01)
pdf("Gvalue.pdf")
g_plot
dev.off()
#可以随机挑选染色体区域进行绘制
QTLplots <- plotQTLStats(SNPset = df_out,
var = "negLog10Pval", #−log10(p-value)
plotThreshold = TRUE, #绘制阈值线
q = 0.01, #FDR
subset = c("Chr1", "Chr8")
)
# 获得候选区间
gregion <- getQTLTable(SNPset = df_filt, method = "Gprime",alpha = 0.01, export = FALSE)
write.table(gregion,file="g_candidate_region.xls",sep="\t",quote=F)
# 提取区间显著位点
g_results <- getQTLTable(SNPset = df_filt, method = "Gprime",alpha = 0.01, export = FALSE)
write.table(g_results[[1]],file="G_0.01.xls",sep="\t",quote=F)
结果图如下：
Gvalue.pdf
runQTLseqAnalysis 函数
SNP-index
# 滑窗计算
df_out <- runQTLseqAnalysis(df_filt,
windowSize = 1e6, #窗口长度，根据基因组大小进行适当调整
popStruc = "F2", #群体类型，F2或者RIL
bulkSize = c(385, 430), #混池单株数目，前者是高值池，后者是低值池
replications = 10000, # bootstrip次数，一般10000次足够了。
intervals = c(95, 99) # 置信区间。
)
#输出日志文件
## Counting SNPs in each window...
## Calculating tricube smoothed delta SNP index...
## Returning the following two sided confidence intervals: 95, 99
## Variable 'depth' not defined, using min and max depth from data: 40-197
## Assuming bulks selected from F2 population, with 385 and 430 individuals per bulk.
## Simulating 10000 SNPs with reads at each depth: 40-197
## Keeping SNPs with >= 0.3 SNP-index in both simulated bulks
## Joining, by = "tricubeDP"
# 绘制delta SNP-index图
delta_SNP_index <- plotQTLStats(SNPset = df_out, var = "deltaSNP", plotIntervals = TRUE)
pdf("delta_snp_index.pdf")
delta_SNP_index
dev.off()
#可以随机选择一条染色体出图
delta_SNP_index <- plotQTLStats(SNPset = df_out,
var = "deltaSNP",
plotIntervals = TRUE，
subset = c("Chr1", "Chr8") #添加要单独展示的染色体编号即可
)
# 获得候选区间
qtlregion <- getQTLTable(SNPset = df_filt, method = "QTLseq",alpha = 0.01, export = FALSE)
write.table(qtlregion,file="index_candidate_region.xls",sep="\t",quote=F)
# 提取区间显著位点
QTL <- getSigRegions(SNPset = df_filt, alpha = 0.01) #
write.table(index_results[[1]],file="index_0.01.xls",sep="\t",quote=F)
结果图：
delta_snp_index.pdf
其他
QTLseqr还提供了绘制标记分布的函数plotQTLStats
p1 <- plotQTLStats(SNPset = df_filt, var = "nSNPs")
# 和之前一样绘制出图即可
参考资料
https://github.com/bmansfeld/QTLseqR/
往期回顾
整理不易，给个好评再走呗！
