Error in running a Samtools mpileup script for getting base frequency information

Question

[jigarnt]

Hi All,

I intend to know the base frequencies at each positions in my BAM files. For the same, I am using samtools mpileup script, which can serve the purpose. While running the script, I am encountering an error which is not understandable by me. So I am reaching out to all the expert bioinformaticians for their assistance. I have posted the error below for your reference:

/Users/Desktop/Tools/samtools-1.3/samtools mpileup -f GCF_000184105.1_Geom_destructans_V1_genomic.fa -r supercont1.1-supercont1.1847 712206_realigned_reads.bam N_american_realigned_reads.bam Nova_Scotia_1_realigned_reads.bam

[E::mpileup] fail to parse region 'supercont1.1-supercont1.1847' with 712206_realigned_reads.bam

Thank you in advance.

Answer 1

[Devon Ryan]

supercont1.1-supercont1.1847 isn't a meaningful region designation. I assume you mean something like supercont1.1:1-1000.

[jigarnt]

Hi Devon,

Thank you for your response. I tried your suggestion of indicating the chromosome positions but it is giving me the same error:

/Users/lindakohn/Desktop/Tools/samtools-1.3/samtools mpileup -f GCF_000184105.1_Geom_destructans_V1_genomic.fa -r supercont1.1:1-1847 712206_realigned_reads.bam N_american_realigned_reads.bam Nova_Scotia_1_realigned_reads.bam

[E::mpileup] fail to parse region 'supercont1.1:1-1847' with 712206_realigned_reads.bam

Also, I want to add that my Reference Fasta is composed of 1847 contigs and they are described as: Supercont1.1 to Supercont1.1847. In that case, can there be any other errors?

[Devon Ryan]

If the contig is called Supercont1.1 in the fasta file, then you need to type that instead of supercont1.1.

[jigarnt]

Hi Devon,

In the file, it is actually supercont1.1. and not Supercont1.1. Would you mind having a look at the fasta file? It is available at the NCBI. It is an 30mb file so would not consume much time of yours. Organisms name is Geomyces destructans

[Devon Ryan]

In my copy (here) the contigs are named NW_012234580.1 and NW_012234581.1 and so on.

[jigarnt]

Hi Devon,

Many thanks for having a look at it. I tried with those numbers before, but got the same error. Is there any other way, to get the base pair information at each positions in your BAM files?

[Devon Ryan]

Please post your BAM file somewhere and I'll have a look.