Expect Library - Bam-Readcount

Question

[tayebwajb]

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!

[ernfrid]

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?

[tayebwajb]

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

Answer 1

[ernfrid]

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

[tayebwajb]

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?

[ernfrid]

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

[tayebwajb]

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

[ernfrid]

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

[xrao]

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!

[sadiq.saleem.09]

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

[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.

[ernfrid]

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

[MAPK]

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.