今天想看一下自己的序列里面会不会有某细菌基因组存在，主要使用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都比对到参考序列上的比对结果，只需要把两个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
附：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
本文使用 文章同步助手 同步
