拿到质控后的数据就可以开始进行数据的比对了。（这里复习一下）需要确定好参考基因组，参考基因组选择标准：质量高、完整、和某一亲本相近（是亲本之一的基因组当然更好了）
比对软件
由于转录组存在可变剪切，所以和DNA重测序数据还是有区别的。简单来说就是DNA重测序数据包含基因组所有的碱基信息（内含子外显子都有），但是RNA数据只有表达数据（外显子，同时存在可变剪切）。所以比对时要注意这一点。DNA数据一般使用BWA，RNA除了BWA还可以使用Hisat2等。
BWA
下载安装
wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
tar jxf bwa-0.7.12.tar.bz2
cd bwa-0.7.12 && make
#编译完成后将软件写入环境变量中
echo 'PATH=$PATH:~/softwarepath/bwa-0.7.12/setup' >> ~/.bashrc && source ~/.bashrc
#测试是否安装好
bwa
Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.12
Contact: Heng Li <lh3@sanger.ac.uk>
Usage: bwa <command> [options]
Command: index index sequences in the FASTA format
mem BWA-MEM algorithm
fastmap identify super-maximal exact matches
pemerge merge overlapping paired ends (EXPERIMENTAL)
aln gapped/ungapped alignment
samse generate alignment (single ended)
sampe generate alignment (paired ended)
bwasw BWA-SW for long queries
fa2pac convert FASTA to PAC format
pac2bwt generate BWT from PAC
pac2bwtgen alternative algorithm for generating BWT
bwtupdate update .bwt to the new format
bwt2sa generate SA from BWT and Occ
Note: To use BWA, you need to first index the genome with `bwa index'.
There are three alignment algorithms in BWA: `mem', `bwasw', and
`aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
first. Please `man ./bwa.1' for the manual.
使用
bwa包含三种算法：
BWA-backtrack: 比对Illumina 序列，read length： ≤100bp。
BWA-SW: 比对 long read，也支持剪切比对。read length：70bp-1Mbp。
BWA-MEM: 比对 long read，也支持剪切比对。read length：70bp-1Mbp。（该算法更新，准确性和速度都更优）
一般使用步骤：
#ref.fa即为参考基因组，genome为建立索引库的前缀，不加-p genome，默认建立索引以ref.fa为前缀
bwa index ref.fa -p genome
#因为BSA分析中的重测序数据大多是150 PE，这里就直接介绍mem算法的用法，别的算法不再详细介绍
bwa mem [-aCHMpP] [-t nThreads] [-k minSeedLen] [-w bandWidth] [-d zDropoff] \
[-r seedSplitRatio] [-c maxOcc] [-A matchScore] [-B mmPenalty] [-O gapOpenPen] \
[-E gapExtPen] [-L clipPen] [-U unpairPen] [-R RGline] [-v verboseLevel] \
db.prefix reads.fq [mates.fq]
常用参数注释：
-t INT 线程数，默认是1。
-M 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
-p 若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R STR 完整的read group的头部，可以用 '\t' 作为分隔符，sam文件会识别。需要指定RG表头是因为同一样品可包括多个不同lane、不同文库的测序结果，此外不同样品比对结果合并时，也需要通过RG进行标记区分。
-T INT 当比对的分值比 INT 小时，不输出该比对结果，这个参数只影响输出的结果，不影响比对的过程。
-a 将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。
SE（单）：
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam
PE（双）：
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam
#实例（双端），sam文件即为我们的比对结果文件。
##不指定样品名
bwa mem ref.fa read1.fq read2.fq > out.sam
##指定样品名，sample即为样品名
bwa mem -M -R "@RG\tID:sample\tLB:sample\tSM:sample" ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log
当我们想要提高比对的速率的时候，可以直接使用shell的基本命令split，或者写python脚本对测序数据进行分割再并行比对（每四行为一组），最后再将比对文件合并起来。
注意：1、因为涉及到后面的变异检测，所以当基因组存在染色体长度大于500M时，一定要对染色体拆分成小于500M的序列再进行比对，在变异检测完毕后再进行合并（拆分基因组的脚本我会放在流程脚本里，这里不再详细记录）；2、一定不能是基因组拆分进行比对，因为数据是无序的，且涉及到多重比对的问题！！！
比对结果格式解读
@HD VN:1.5 SO:coordinate
@SQ SN:Chr01 LN:91897250
@RG ID:sample LB:sample SM:sample
@PG ID:bwa CL:~/softwarepath/bwa-0.7.12/setup/bwa mem -M -R "@RG\tID:sample\tLB:sample\tSM:sample" ref.fa read1.fastq read2.fastq > out.sam 2> out.sam.log PN:bwa VN:0.7.17-r1188
E150015117L1C019R0062820889 433 Chr01 186 0 72M78H ChrD02_S6 89876638 0 AAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCTA ?A+<80-?C5-@A5BC<7>A<@CC-2A9CCB67@2BCC/;8:@CCC=66@CCA<>=BCCCAB68B@C<?<=C NM:i:0 MD:Z:72 MC:Z:150M AS:i:72 XS:i:72 RG:Z:CY619 SA:Z:ChrD05_S5,59252485,-,51S99M,0,1;
E150015079L1C036R0342398185 161 Chr01 361 2 6S14M1I86M1I14M4D28M ChrD02_S1 14286 0 CCTAAAAAACCCTAAACCCTAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCAAAAAACCCTAAAAAAACCCTAAACCCTAAACCCTAAACCCTAAACCATAAACCCCTAAACCCTAAAAAACCCTAAACCCTAAACCCTAAACCCT B-AB<CBA/65@=5B<.C96C4:B<#AB2=8);18B7/;=??C86>)@AB<8C4<CCA2<C69@>=?2A=A76=6/867/C=4>A--=(2B=A1;@@/9.#@'4C>).@?92<C421AAC(B2@99-C6#=&'<<9@4C2'B.B$B/CC% NM:i:7 MD:Z:95C18^CCCT28 MC:Z:147M3S AS:i:113 XS:i:109 RG:Z:CY619
E150015079L1C003R0260125142 81 Chr01 368 0 38M1D112M ChrD03 10858 0 AAACCGTAAACCCTAAAGCCTAAATCCTAAACCCTAAAGCAAAAAACCCTAAAAAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAAAGCTAAACCCTAAACGCTAAACCCTAAACC @2C;15CCCC26<BCAB$?BCBC>/<.@@CB9C.C?CC#CBBA6BC:C6ACCBCC>CC=@CCCA7A>CACAC6CCACC?1?CCC=B9CCCCC.2-CCC<;C4CCCC;5?CCCCC-CCCAC&(CCC;A*C5CAC7C'9BBCCC?7CCC>8B NM:i:8 MD:Z:5C11C6C13^C0C81C0C13C14 MC:Z:3M1D147M AS:i:108 XS:i:106 RG:Z:CY619
#@开头的列都是表头信息行，不能少
HD：（排序类型）
头部区第一行：VN是格式版本；SO表示比对排序的类型，有unknown（default），unsorted，queryname和coordinate几种。
SQ：染色体号\t染色体长度
RG：样本信息
PG: 比对命令行，可以直接由此查到比对方式，参考基因组，测序数据，样本名等信息
#比对内容行每列内容注释：
第一列：QNAME。比对的序列名称，类似于E150015117L1C019R0062820889。
第二列：FLAG Bwise FLAG，比对类型：paring，strand，mate strand等。（这个值需要进行计算）
第三列：RENAME 比对上的参考序列名，类似于Chr01
第四列：POS 1-Based的比对起始物理位置
第五列：MAPQ 比对质量
第六列：CIGAR Extended CIGAR string（操作符：MIDNPSH=X）比对结果信息，包含数字和字符。
第七列：MRNM 相匹配的另外一条序列，比对上的参考序列名，“=”：比对到一个位置；“*”：无匹配；染色体号：比对多个位置。
第八列：MPOS 1-Based leftmost Mate Position 与该reads对应的mate pair reads的比对位置
第九列：ISIZE 插入片段长度
第十列：SEQ 和参考序列在同一个链上比对的序列（若比对结果在负义链上，则序列是其反向重复序列，反向互补序列）
第十一列：QUAL 比对序列的质量（ASCII-33=Phred base quality）reads碱基质量值
可选列：以TAG：TYPE：VALUE的形式提供额外的信息。
比对文件格式处理
samtools
下载安装
# 建议直接conda安装
conda create -n samtools
conda activate samtools
conda install -c bioconda samtools -y
#测试是否安装成功
samtools
Program: samtools (Tools for alignments in the SAM format)
Version: 1.3.1 (using htslib 1.3.1)
Usage: samtools <command> [options]
Commands:
-- Indexing
dict create a sequence dictionary file
faidx index/extract FASTA
index index alignment
-- Editing
calmd recalculate MD/NM tags and '=' bases
fixmate fix mate information
reheader replace BAM header
rmdup remove PCR duplicates
targetcut cut fosmid regions (for fosmid pool only)
addreplacerg adds or replaces RG tags
-- File operations
collate shuffle and group alignments by name
cat concatenate BAMs
merge merge sorted alignments
mpileup multi-way pileup
sort sort alignment file
split splits a file by read group
quickcheck quickly check if SAM/BAM/CRAM file appears intact
fastq converts a BAM to a FASTQ
fasta converts a BAM to a FASTA
-- Statistics
bedcov read depth per BED region
depth compute the depth
flagstat simple stats
idxstats BAM index stats
phase phase heterozygotes
stats generate stats (former bamcheck)
-- Viewing
flags explain BAM flags
tview text alignment viewer
view SAM<->BAM<->CRAM conversion
depad convert padded BAM to unpadded BAM
使用
这里主要介绍如何将sam文件处理成二进制bam文件，建立索引，排序，去重，以及统计比对信息
#sam文件查看
less -S out.sam
#bam文件查看
samtools view -hS out.sort.rmdup.bam |less -S
##以上二者区别就是是否为二进制文件
#sam转bam
samtools view -hb out.sam >out.bam
#拆分并行比对后会涉及到bam合并的问题,merge也可以合并
samtools cat 1.bam 2.bam ... -o out.bam
#bam排序,先建索引（生成*.bam.bai文件，当染色体长度过长时则会出现索引无法建立的情况）再排序
samtools index out.bam
samtools sort out.bam -o out.sort.bam
#去重，s:SE,S:PE,默认双端
samtools rmdup -S out.sort.bam out.sort.rmdup.bam
#统计比对信息
samtools flagstat out.sort.rmdup.bam > flagstat.txt
#flagstat.txt文件内容
7524 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7523 + 0 mapped (99.99% : N/A)
7524 + 0 paired in sequencing
3760 + 0 read1
3764 + 0 read2
7510 + 0 properly paired (99.81% : N/A)
7522 + 0 with itself and mate mapped
1 + 0 singletons (0.01% : N/A)
12 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
第一行：QC pass的reads的数量为7514，未通过QC的reads数量为0，意味着一共有7524*2条reads；
第四行：重复reads的数量，QC pass和failed，这里都是0条
*第五行：比对到参考基因组上的reads数量，这里是共7523条比对上了参考基因组，比对率为99.99%；
第六行：paired reads数据数量，这里是7524对；
第七行：read1比对上的数量；
第八行：read2比对上的数量；
*第九行：双端都正确匹配到参考序列的reads数量，这里是7510，比对率为99.81%；
第十行：双端都匹配到了参考序列上的reads数量(但是并不一定比对到同一条染色体上)；
*第十一行：双端中只有一条与参考序列相匹配的数量；
*第十二行：双端比对到不同染色体的数量；
第十三行：双端比对到不同染色体的且比对质量值大于5的数量。
## *为主要数据关注行
#覆盖率统计
/work1/Software/samtools-1.4/setup/bin/samtools depth -a out.sort.rmdup.bam > out_depth.txt 2>out_depth.statlog
#out_depth.txt文件内容
#染色体号 物理位置 覆盖碱基数
Chr01 1 0
Chr01 2 10
Chr01 3 0
Chr01 4 0
## a为有碱基覆盖的位点总长度
a=`awk '$3!=0' out_depth.txt|wc -l `
## b是基因组总长度
samtools faidx ref.fa
b=`awk '{sum+=$2}END{print sum} ref.fai`
## a/b 即为覆盖率了
往期回顾
参考资料
https://blog.csdn.net/genome_denovo/article/details/78712972
