Samtools Header From A Bacterial Genome Ba File

Question

[Sakti]

Hello dear Bio masters,

I'm analyzing some reads from a bacterial genome (E. coli), and I'm trying to get the pileup to calculate coverage for a specific region from a bam file. However, whenever I check the header of my file it has the form:

samtools view -H coli.sorted.bam
@SQ    SN:gi|49175990|ref|NC_000913.2|    LN:4639675
@PG    ID:bwa    PN:bwa    VN:0.5.9-r16

And when I try to retrieve the information from my region (1670-2550) of interest using the "chromosome name" after the SN

samtools view -u coli.sorted.bam gi|49175990|ref|NC_000913.2|:1670-2550 | samtools pileup - > pileup.txt

I get:

-bash: 49175990: command not found
-bash: ref: command not found
-bash: NC_000913.2: command not found
-bash: :190-255: command not found
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[bam_pileup] fail to read the header: non-exisiting file or wrong format.
[main_samview] fail to get the reference name. Continue anyway

So it is not accepting the gi|49175990|ref|NC_000913.2| as a name for search. If I omit it I get:

samtools view -u coli.sorted.bam 1670-2550 | samtools pileup - > pileup.txt
[main_samview] fail to get the reference name. Continue anyway.

And my file is empty. What am I doing wrong????? Any help will be highly appreciated.

Thanks!!!

Answer 1

[Docroberson]

I'm not sure if the problem is with samtools. The bash errors come from bash seeing the vertical pipes and trying to pipe the output of one command into another.

Not sure if this actually will fix it, but try quoting it instead:

samtools view -u coli.sorted.bam "gi|49175990|ref|NC_000913.2|:1,670-2,550" | samtools pileup - > pileup.txt

That should differentiate the name with pipes in it from the actual pipe to samtools pileup.

[Sakti]

Thank you, yes, that worked :)

Answer 2

[Drio]

In your first one, your shell is trying to expand the command prior to execute it. Also you are not passing the reference genome in the pileup. Try:

 samtools view -u -h coli.sorted.bam "gi|49175990|ref|NC_000913.2|:1670-2550" | samtools pileup -vcf ./ref.fa -

where ref.fa is the fasta file for your reference genome.

You are also missing the reference genome in your second one. Try:

 samtools view -u -h coli.sorted.bam 1670-2550 | samtools pileup -vcf ./ref.fa -

[Sakti]

Yes, that worked :)

[Drio]

If it worked, could you please mark my reply as question answered ?

[Sakti]

Done :) Thanks!