以水稻为例教你如何使用BSA方法进行遗传定位（上篇）

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`; dosample=${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-filterbcftools 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语言下游分析的基础。