文件:全基因组二代测序双端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比对生成的为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 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文件,并将结果保存至指定的输出-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
如果觉得我的文章对您有用,请随意打赏。你的支持将鼓励我继续创作!