This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Expect Library - Bam-Readcount

I have installed bam-readcount using git clone --recursive as per the instructions on github.

I am now using the command below to obtain read-count metrics 

bam-readcount -q 1 -b 20 -f reference.fasta -l varScan.variants BAM_FILE >varScan.variants.readcounts

However, I get an empty output with a message

Expect library: SX215 in BAM

My bam file has the @RG field (tab-delimited) as follows 

@RG     ID:SampleID     PL:illumina     PU:L006     LB:SX215        SM:SampleID

This was added using Picard. Could anyone be knowing what the problem here is? Thanks!

• 7,931 views • link
updated by ernfrid ernfrid • written by tayebwajb tayebwajb

The "Expect library: $name" message is just a status message. It does not indicate something is wrong with the BAM. The fact you were able to call variants with VarScan suggests that you should be able to also get readcounts, so that is weird. What does your varScan.variants file look like?

• 0 views • link
updated by Ram Ram • written by ernfrid ernfrid 
chrom    position    ref    var    normal_reads1    normal_reads2    normal_var_freq    normal_gt    tumor_reads1    tumor_reads2    tumor_var_freq    tumor_gt    somatic_status    variant_p_value    somatic_p_value    tumor_reads1_plus    tumor_reads1_minus    tumor_reads2_plus    tumor_reads2_minus    normal_reads1_plus    normal_reads1_minus    normal_reads2_plus    normal_reads2_minus
chr1    4843531    G    A    11    0    0%    G    4    5    55.56%    R    Somatic    1    0.008126935    0    4    2    3    0    11    0    0
• 0 views • link
updated by Ram Ram • written by tayebwajb tayebwajb

1 answer

bam-readcount expects that you will provide it a file that is formatted as follows:

chr   start_pos   end_pos

There should be no header. Each field should be tab separated and coordinates should be 1-based.

For this site your file should look like:

 chr1   4843531   4843531

where each field is separated by a tab

• 406 views • link
written by ernfrid ernfrid

Thanks, that worked. Also, I get a warning that "Couldn't find single-end mapping quality. Check to see if the SM tag is in BAM". Does this mean that in this case bam-readcount is only using reads mapped as pairs?

• 0 views • link
updated by Ram Ram • written by tayebwajb tayebwajb

All it means is that the single-end mapping quality field of the output will not be populated.

• 47 views • link
written by ernfrid ernfrid

Does bam-readcount work with indels? From my understanding, it should so long as you specify start and stop, right? Although this is not straightforward from the output of VarScan2. Found this thread that says bam-readcount only works with point mutations

• 0 views • link
updated by Ram Ram • written by tayebwajb tayebwajb

bam-readcount will work for indels. Make sure to give it a little bit of padding around the coordinates (+/- 2bp or so) to attempt to ensure you've picked them up

• 0 views • link
written by ernfrid ernfrid

Thank you for this useful information. But how to supply the positions for indel? So for an insertion of "AAAGGG" at chr1 11175116, I will just specify it as chr1 11175116-2 11175116+2? I just noticed an option of bam-readcount, which is -i (generate indel centric readcounts. Reads containing insertions will not be included in per-base counts). Not sure if it is related.

Thank you in advance!

• 0 views • link
updated by Ram Ram • written by xrao • xrao

@ernfrid and @tayebwajb How to do that? I am facing a similar problem and get an empty output. Thanks!

• 0 views • link
written by sadiq.saleem.09 • sadiq.saleem.09

I mean regarding the first suggestion, how to remove the header. As I remember without header BAM file is supposedly useless. Thank you once again.

• 0 views • link
written by sadiq.saleem.09 • sadiq.saleem.09

Only in the file of positions should there be no header. You definitely want to keep the header of your BAM file.

• 0 views • link
written by ernfrid ernfrid

Do you need reference.fasta here? I need to get pile up reads for ref and alt alleles using chr and position as input in bam files.

• 0 views • link
updated by Ram Ram • written by MAPK MAPK

Log in to answer this question.

More posts like this