This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get part of a sequence from a .bam file, considering the reference genome position.

Hi everybody.

I am trying to get part of a sequence from a .bam file, considering the reference genome position. Like, which sequence from the .bam file returns at position 50,000 to 50,400 nucleotides... Can someone help me? I need all in details, since I never work with this.

assembly
• 3,445 views • link
updated by Petr Ponomarenko Petr Ponomarenko • written by valopes • valopes

Thank you very much. Can I make you more questions? So the ony thing that I need to change in this command is "ref.fasta" for my genome.fata and "in1.bam" for my.bam? And of course the coordinates...

• 0 views • link
written by valopes • valopes

I did like this and I got this message

[bam_parse_region] fail to determine the sequence name. [mpileup] malformatted region or wrong seqname for US-18.bam

• 0 views • link
written by valopes • valopes

You will have to run following command to create the indexes first.

samtools faidx ref.fasta
• 0 views • link
written by Sej Modha Sej Modha

1 answer

samtools documentation (i.e. http://www.htslib.org/doc/samtools.html) recomends to use sorted bam files for many tools like "view", "flagstat", "bedcov" and other. Also sorted bam files are a bit smaller and other programs work faster with them. So I usually sort and index bam files. Then I use samtools view to see the region of interest:

samtools sort -T tmp.in1 -o in1.sorted.bam in1.bam
samtools index in1.sorted.bam
samtools view in1.sorted.bam chr3:50000-50400

Make sure your prefix for the temporary files (option -T in samtools sort) is different for different bam files if you sort multiple files in parallel otherwise your sorted files will contain corrupted data.

• 0 views • link
written by Petr Ponomarenko Petr Ponomarenko

Log in to answer this question.

More posts like this