BWA比对及Samtools提取mapped reads

筛选出能比对到参考序列的reads,形成新的fastq文件

文件:全基因组二代测序双端fastq文件,参考序列fasta文件

目的:筛选出能比对到参考序列的reads,形成新的fastq文件

使用软件:BWA,Samtools


一、构建索引

bwa index ref.fasta
可选参数:
-p 索引文件前缀名
-a bwtsw :参考基因组大于2G的时候添加该参数

生成的索引文件:ref.fasta.amb、ref.fasta.ann、ref.fasta.bwt、ref.fasta.pac、ref.fasta.sa


二、bwa比对及samtools转为bam文件

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

bwa mem ref.fas sample_R1.fq.gz sample_R2.fq.gz | samtools view -Sb - > sample.bam


三.samtools根据比对情况提取

#提取比对到参考序列上的比对结果
samtools view -bF 4 sample.bam > sample_mapped.bam
#提取paired reads中两条reads都比对到参考序列上的比对结果,只需要把两个4+8的值12作为过滤参数即可;
samtools view -bF 12 sample.bam > sample_all.mapped.bam
#提取没有比对到参考序列上的比对结果
samtools view -bf 4 sample.bam > sample_unmapped.bam


四、bam文件转化为fastq文件

将输入的bam文件转化为fastq文件,并将结果保存至指定的输出-1-2-0-0-s中.

samtools fastq [options] in.bam

其对序列的分类依据是序列末尾的READ1 flag和READ2 flag,flag有3类:

1:Only READ1 is set.

2:Only READ2 is set.

O:Either both READ1 and READ2 are set;or neither is set.


对于PE测序reads,同时指定-1 R1.fq-2 R2.fq-s singleton.fq 时,samtools会将flag1和flag2配对的序列分别输出到-1-2指定的文件,对于无法匹配的flag 1/2,全部的flag 0都会保存到-s指定的文件中。

${ samtools} fastq -n \
-1 R1.fq.gz \
-2 R2.fq.gz \
-s unmapped_singleton.fq.gz \
sample_mapped.bam


参考链接:https://zhuanlan.zhihu.com/p/514094203

  • 发表于 2022-12-20 14:27
  • 阅读 ( 5391 )
  • 分类:软件工具

0 条评论

请先 登录 后评论
星莓
星莓

生物信息工程师

58 篇文章

作家榜 »

  1. omicsgene 698 文章
  2. 安生水 347 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 82 文章
  6. 红橙子 78 文章
  7. rzx 74 文章
  8. CORNERSTONE 72 文章