Expect Library - Bam-Readcount
1
0
Entering edit mode
8.5 years ago
tayebwajb ▴ 110

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!

• 4.7k views
ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode
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
ADD REPLY
1
Entering edit mode
8.5 years ago
ernfrid ▴ 400

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

ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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!

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 807 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6