今天想看一下自己的序列里面會不會有某細(xì)菌基因組存在，主要使用BWA和Samtools：

bwa主要用于將低差異度的短序列與參考基因組進(jìn)行比對。主要包含三種比對算法：backtrack、SW和MEM，第一種只支持短序列比對（<100bp），后兩種支持長序列比對(70bp~1M)，并支持分割比對（split alignment）。MEM算法是最新的也是官方推薦的。

BWA-MEM 是一種新的比對算法，用于將測序 reads 或者組裝后 contigs 比對至大型參考基因組，例如人參考基因組。它會自動選擇局部比對和 end-to-end 比對模式，支持PE reads 比對和嵌合體比對。該算法對測序錯(cuò)誤有良好的穩(wěn)定性，適用的reads 長度范圍較廣，從70bp至幾Mb。

bwa的工作原理

所有的比對工具均基于相同的原則：

1. 從參考基因組建立一個(gè)索引

2. 將FASTA和FASTQ文件中的序列同索引進(jìn)行比對

一．ncbi上下載了wol細(xì)菌的100條序列，作為ref。

二．BWA比對

1.先是構(gòu)建索引：

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

2.bwa比對及samtools轉(zhuǎn)為bam文件,并根據(jù)比對情況進(jìn)行提取

bwa比對生成的為sam（sequence Alignment mapping)文件，將SAM轉(zhuǎn)換為二進(jìn)制對應(yīng)的BAM格式。二進(jìn)制格式對于計(jì)算機(jī)程序來說更容易使用。要將SAM轉(zhuǎn)換為BAM，使用samtools view命令。

在醫(yī)院服務(wù)器上用轉(zhuǎn)錄組的數(shù)據(jù)成功運(yùn)行，學(xué)校服務(wù)器上報(bào)錯(cuò)，見后面。

三．samtools根據(jù)比對情況提取：

四．Flag參數(shù)：

flag是一種描述read比對情況的標(biāo)記，對于雙端比對的數(shù)據(jù)，生成的BAM文件中，R1端序列和R2端序列的標(biāo)識符是一樣的，之前一直不知道如何根據(jù)bam文件區(qū)分哪條序列是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

五.提取特定位置上的比對結(jié)果

六．將sam文件、bam文件、fastq文件之間轉(zhuǎn)換

附：samtools功能蠻強(qiáng)大的，功能很多，都可以單獨(dú)寫一篇了，參考的里面有非常詳細(xì)的記錄。

七．錯(cuò)誤

錯(cuò)誤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出了錯(cuò)誤。

發(fā)現(xiàn)錯(cuò)誤2

錯(cuò)誤2.單獨(dú)運(yùn)行bwa程序：

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

一時(shí)找不到原因，但換了服務(wù)器運(yùn)行就沒有問題，暫時(shí)記錄一下。

參考：

bwa官網(wǎng)：http://bio-bwa.sourceforge.net/

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

本文使用 文章同步助手 同步