BWA比对及Samtools提取目标序列

今天想看一下自己的序列里面会不会有某细菌基因组存在，主要使用BWA和Samtools：

bwa主要用于将低差异度的短序列与参考基因组进行比对。主要包含三种比对算法：backtrack、SW和MEM，第一种只支持短序列比对（<100bp），后两种支持长序列比对(70bp~1M)，并支持分割比对（split alignment）。MEM算法是最新的也是官方推荐的。

BWA-MEM 是一种新的比对算法，用于将测序 reads 或者组装后 contigs 比对至大型参考基因组，例如人参考基因组。它会自动选择局部比对和 end-to-end 比对模式，支持PE reads 比对和嵌合体比对。该算法对测序错误有良好的稳定性，适用的reads 长度范围较广，从70bp至几Mb。

bwa的工作原理

所有的比对工具均基于相同的原则：

1. 从参考基因组建立一个索引

2. 将FASTA和FASTQ文件中的序列同索引进行比对

一．ncbi上下载了wol细菌的100条序列，作为ref。

二．BWA比对

1.先是构建索引：

./bwa index /share/home/myname/wol/wol100.fas -p wol
-p索引文件前缀名
-a bwtsw :参考基因组大于2G的时候添加该参数

生成的索引文件：wol100.fas.amb、wol100.fas.ann、wol100.fas.bwt  wol100.fas.pac、wol100.fas.sa

2.bwa比对及samtools转为bam文件,并根据比对情况进行提取

bwa比对生成的为sam（sequence Alignment mapping)文件，将SAM转换为二进制对应的BAM格式。二进制格式对于计算机程序来说更容易使用。要将SAM转换为BAM，使用samtools view命令。

在医院服务器上用转录组的数据成功运行，学校服务器上报错，见后面。

bwa mem wol100.fas sample_R1.fq.gz sample_R2.fq.gz | /path/to/samtools view -S -bF 4 - > sample.bam

三．samtools根据比对情况提取：

#提取比对到参考序列上的比对结果
samtools view -bF 4 test.bam > test.F.bam
 
#提取paired reads中两条reads都比对到参考序列上的比对结果；
#提取paired reads中两条reads都比对到参考序列上的比对结果，只需要把两个4+8的值12作为过滤参数即可
samtools view -bF 12 test.bam > test.F12.bam
 
#提取没有比对到参考序列上的比对结果
samtools view -bf 4 test.bam > test.f.bam

四．Flag参数：

flag是一种描述read比对情况的标记，对于双端比对的数据，生成的BAM文件中，R1端序列和R2端序列的标识符是一样的，之前一直不知道如何根据bam文件区分哪条序列是R1端，哪条序列是R2端，代表R1端和R2端的信息都存储在flag中，即bam文件的第二列；

在bam文件格式中定义了各种flag代表的意思，一种12种，可以搭配使用。

-f 正确比对 only include reads with all  of the FLAGs in INT present [0]

-F NOT正确比对 only include reads with none of the FLAGS in INT present [0]

查询flag值含义：samtools flag 4

更多Flag信息见：http://www.htslib.org/doc/samtools-flags.html

五.提取特定位置上的比对结果

# 提取bam文件中比对到chr1上的比对结果，并保存到sam文件格式
samtools view test.bam chr1 > test.chr1.sam

# 线粒体
samtools view test.bam chrM > test.chrM.sam

# 提取chr1上能比对到30k到100k区域的比对结果
samtools view test.bam chr1:30000-100000 > test.chr1_30k-100k.sam

六．将sam文件、bam文件、fastq文件之间转换

## bam文件转fastq文件

samtools fastq -N -1 sample_P1.fq -2 sample_P2.fq sample.bam

## sam文件转bam文件：

samtools view -bS test.sam > test.bam

## bam文件转sam文件

附：samtools功能蛮强大的，功能很多，都可以单独写一篇了，参考的里面有非常详细的记录。

七．错误

错误1.

[M::bwa_idx_load_from_disk] read 0 ALT contigs

[M::process] read 67652 sequences (10000283 bp)…

[main_samview] fail to read the header from “-”.

中划线“-”是前一句命令的，也就是bwa出了错误。

发现错误2

错误2.单独运行bwa程序：

/opt/gridview//pbs/dispatcher/mom_priv/jobs/21873.node1.SC: line 10: 20117 Bus error  (core dumped)

一时找不到原因，但换了服务器运行就没有问题，暂时记录一下。

参考：

bwa官网：http://bio-bwa.sourceforge.net/

samtools命令详解：https://blog.csdn.net/qq_27390023/article/details/121164168