BSA虽然不是我最早接触的高通量数据分析项目（最早的是RNA-seq），但是却是我最早独立开展的数据分析项目， 我曾经专门写过一篇文章介绍如何使用SHOREMap做拟南芥的EMS诱变群体的BSA分析
在遗传定位上，相对于GWAS和binmap，BSA是一个比较省钱的策略，它只需要测两个亲本和后代中两个极端差异群体即可，但是它对实验设计，表型考察，样本挑选都有比较高的要求。如果你的表型差异并不是泾渭分明，那么还是不要用BSA比较合适。这方面知识，建议阅读徐云碧老师2016年发表在PBJ上的"Bulked sample analysis in genetics, genomics and crop improvement", 这篇教程侧重于实际分析，而非理论讲解。
用于讲解的文章题为"Identification of a cold-tolerant locus in rice (Oryza sativa L.) using bulked segregant analysis with a next-generation sequencing strategy", 发表在Rice上。 文章用的耐寒(Kongyu131)和不耐寒(Dongnong422)的两个亲本进行杂交，得到的F2后代利用单粒传法得到RIL(重组自交系)群体.之后用BWA+GATK识别SNP，计算ED和SNP-index作图。
这次实战的目标也就是得到文章关键的两张图：
环境准备
新建一个项目，用于存放本次分析的数据和结果
mkdir -p rice-bsa
cd rice-bsa
后续分析还需要用到如下软件:
- wget: 一般Linux系统会自带，用于下载数据
- seqkit: 多能的序列处理软件
- fastp: 数据质控
- bwa: 数据比对到参考基因组
- samtools: 处理sam/bam文件
- sambamba: 标记重复序列
- bcftools: 处理vcf/bcf，能用于变异检测
- R: 数据分析
我们需要分别下载参考基因组(IRGSP)数据和测序数据。
# rice-bsa目录下
mkdir -p ref && cd ref
wget https://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_genome.fasta.gz
gunzip IRGSP-1.0_genome.fasta.gz
# 建立索引
bwa index IRGSP-1.0_genome.fasta
之后 根据文章提供的编号，SRR6327815, SRR6327816, SRR6327817, SRR6327818 ，我们到ENA 查找对应的下载链接进行数据下载。
# rice-bsa目录下
mkdir -p data && cd data
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR632/005/SRR6327815/SRR6327815_{1,2}.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR632/006/SRR6327816/SRR6327816_{1,2}.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR632/007/SRR6327817/SRR6327817_{1,2}.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR632/008/SRR6327818/SRR6327818_{1,2}.fastq.gz
因为在NCBI上下载的是SRA格式，需要做数据转换，而从ENA上可以直接下载压缩过的fastq文件，所以我更偏好于ENA。
上游预处理
对于测序数据而言，为了能够获得后续用于计算SNP-index或者ED的SNP信息，我们需要将二代测序得到的数据回贴到参考基因组上，然后利用变异检测软件找到每个样本和参考基因组的区别。
数据质控
原始数据中可能有一部分序列质量不够高，会影响后续分析，比如说测序质量不够，或者说存在接头。通常我们会都会对原始数据进行一波过滤，这里用的是fastp，优点就是快。
# rice-bsa目录下
mkdir -p 01-clean-data
fastp -i data/SRR6327815_1.fastq.gz -I data/SRR6327815_2.fastq.gz -o 01-clean-data/SRR6327815_1.fastq.gz -O 01-clean-data/SRR6327815_2.fastq.gz
fastp -i data/SRR6327816_1.fastq.gz -I data/SRR6327816_2.fastq.gz -o 01-clean-data/SRR6327816_1.fastq.gz -O 01-clean-data/SRR6327816_2.fastq.gz
fastp -i data/SRR6327817_1.fastq.gz -I data/SRR6327817_2.fastq.gz -o 01-clean-data/SRR6327817_1.fastq.gz -O 01-clean-data/SRR6327817_2.fastq.gz
fastp -i data/SRR6327818_1.fastq.gz -I data/SRR6327818_2.fastq.gz -o 01-clean-data/SRR6327818_1.fastq.gz -O 01-clean-data/SRR6327818_2.fastq.gz
序列比对
之后将各个样本的序列回帖到参考基因组
# rice-bsa目录下
mkdir -p 02-read-align
NUMBER_THREADS=80
REFERENCE=ref/IRGSP-1.0_genome.fasta
for i in `ls 01-clean-data/`; do
sample=${i%%_*}
(bwa mem -M -R "@RG\\tID:${sample}\\tSM:${sample}\\tPL:ILLUMINA" \
-t $NUMBER_THREADS $REFERENCE 01-clean-data/${sample}_1.fastq.gz 01-clean-data/${sample}_2.fastq.gz || echo -n 'error' ) \
| samtools sort -@ 20 -o 02-read-align/${sample}_sort.bam -
samtools index -@ 20 02-read-align/${sample}_sort.bam
done
这里用到了稍微比较高级的shell操作
- 变量名替换
sample=${i%%_*} - 逻辑判断:
|| - 子shell:
() - for循环:
for do done
接着我门需要去除重复序列，这里的重复序列指的是拥有相同位置信息，且序列也一模一样的read，通常是由PCR扩增引起。过滤重复序列的目的是为了提高变异检测的准确性，如果一条read上出现的“变异”其实是在第一轮PCR扩增时引入的错误，那么后续的扩增只会让这个错误一直保留着，随后测序的时候这条许多又被多次测到，那么在后续的分析中由于多次出现，就有可能会变异检测软件当作真实变异。
这里用sambamba，因为它的速度比较快，且结果和picard一模一样。
# rice-bsa目录下
for i in `ls 02-read-align/*_sort.bam`; do
sample=${i%%_*}
sambamba markdup -r -t 10 ${sample}_sort.bam ${sample}_mkdup.bam
done
变异检测
变异检测最常见的就是bcftools, freebayes和GATK. 这里用的是BCFtools，主要原因还是它的速度比较快。
这里为了让他的速度更快，我用了--region参数分染色体并行，由于水稻有12条染色体，相当于提速了12倍
# rice-bsa目录下
mkdir -p 03-variants
ls -1 02-read-align/*_mkdup.bam > bam_list.txt
seqkit seq -n ref/IRGSP-1.0_genome.fasta | \
while read region
do
bcftools mpileup -f ref/IRGSP-1.0_genome.fasta \
--redo-BAQ --min-BQ 30 \
--per-sample-mF \
--annotate FORMAT/AD,FORMAT/DP \
--regions ${region} \
-Ou --bam-list bam_list.txt | \
bcftools call -mv -Ob -o 03-variants/${region}.bcf &
done
之后将拆分运行的结果合并到单个文件中
# rice-bsa目录下
mkdir -p 04-variant-filter
bcftools concat --naive -o 04-variant-filter/merged.bcf 03-variants/*.bcf
变异过滤
得到VCF文件还需要进行一些过滤，来提高变异的准确性. 这个需要根据具体的项目来进行, 但是有一些固定的指标可以用来对结果进行过滤，例如
- 测序深度
- 非参考基因组的高质量read数
- 是否和indel紧邻，通常和indel比较近的snp都不可靠
- 平均的比对质量,
举个例子:
# rice-bsa目录下
cd 04-variant-filter
bcftools filter -g3 -G10 -e'%QUAL<10 || (RPB<0.1 && %QUAL<15) || (AC<2 && %QUAL<15) || MQ < 30 || MQSB <=0.1' merged.bcf > filter.vcf
之后，我只选择了snp用于后续分析
# 04-variant-filter目录下
bcftools view -i 'TYPE="snp" & N_ALT =1 & STRLEN(ALT) = 1' filter.vcf > snps.vcf
最终留下了将近80w的SNP。这就是后续R语言下游分析的基础。
版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。
