将经过过滤后得到的高质量数据比对到参考基因组上，并进行排序和去重复等处理，用于后续的变异检测。

软件准备

bwa
samtools
picard

文件准备

参考基因组：genome.fasta
测序数据文件：sample1_1.fq.gz、sample1_2.fq.gz、sample2_1.fq.gz、sample2_2.fq.gz

构建基因组index

重测序很多软件不能直接读取基因组，需要对基因组构建index才能使用。

samtools faidx genome.fasta #基因组的统计信息

bwa index  genome.fasta #bwa比对需要的索引文件

picard CreateSequenceDictionary R=genome.fasta #生成字典序文件

比对

采用bwa (Li,H.,and R. Durbin等，2009) mem程序将经过过滤后得到的高质量数据比对到参考基因组上，比对的参数均按照bwamem的默认参数。采用samtools软件对sam文件进行排序并转换为bam文件。测序数据的生成过程涉及到文库的扩增和簇的形成，这两个步骤容易产生一些 Duplicates（即PCR Duplicates和Optical Duplicates），这些Duplicates 不能作为变异检测的证据。采用Picard软件包中的“MarkDuplicates” 去除Duplicates，即如果多个Paired Reads比对后具有相同的染色体坐标，则仅保留分值最高的Paired Reads。

bwa mem -t 2 \ # -t 线程数
 -R '@RG\tID:{sample}\tSM:{sample}\tPL:illumina' \ # 添加表头信息
 genome.fasta \  参考基因组index路径
 S1_1.fq.gz   \  fq1数据
S1_2.fq.gz    \  fq2数据
2>S1.bwa.log |  \ 写到日志文件
samtools sort -@ 2 -m 1G -o S1.sort.bam -   \ 排序

##生成文件：S1.sort.bam

#去除重复
picard -Xmx4g  \ #设置内存
MarkDuplicates I=S1.sort.bam \ #输入文件
O=S1.sort.rmdup.bam  \ #输出文件
CREATE_INDEX=true  \ #创建bam index
REMOVE_DUPLICATES=true M=S1.marked_dup_metrics.txt

#比对率统计
samtools  flagstat  S1.sort.bam > S1.sort.bam.flagstat

#覆盖率统计
samtools  coverage  S1.sort.bam

欢迎关注Bioinfor 生信云微信公众号！