生物信息tsRNA-seqRNA-seq

BWA比对及Samtools提取目标序列

2022-05-12  本文已影响0人  吕强强学生信

今天想看一下自己的序列里面会不会有某细菌基因组存在,主要使用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

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

上一篇下一篇

猜你喜欢

热点阅读