bam stats F1F2/F1R2/R1F2/R1R2/F2F1/F2R1/R2F1/R2R1

Question

[14134125465346445]

Hi,

I would like to know the statistics of the different read pair orientations for a bam file aligned with bwa mem.

I tried bamtools stats, but it does not split the counts into all possible orientations: F1F2/F1R2/R1F2/R1R2/F2F1/F2R1/R2F1/R2R1

Is there a piece of software that will give me the stats of each orientation possible?

Thx

Answer 1

[Pierre Lindenbaum]

using bioalcidae: https://github.com/lindenb/jvarkit/wiki/BioAlcidae

var count={};

while(iter.hasNext())
    {
    var rec = iter.next();
    if(!rec.getReadPairedFlag()) continue;
    if(rec.getReadUnmappedFlag()) continue;
    if(rec.getMateUnmappedFlag()) continue;
    /* I'm not sure you want this */
    if(rec.getAlignmentStart() >= rec.getMateAlignmentStart()) continue;

    if(rec.getReadNegativeStrandFlag())
        {
        key="R";
        }
    else
        {
        key="F";
        }

    key+= (rec.getFirstOfPairFlag()?"1":"2");

    if(rec.getMateNegativeStrandFlag())
        {
        key+="R";
        }
    else
        {
        key+="F";
        }

    key+=(rec.getFirstOfPairFlag()?"2":"1");

    if(!(key in count))
        {
        count[key]=1;
        } 
    else
        {
        count[key]++;
        }
    }

for(var key in count)
    {
    out.println(key+" "+count[key]);
    }

usage:

 java -jar src/jvarkit/dist/bioalcidae.jar -f file.js in.bam

or

 samtools view -F (something) -bu in.bam | java -jar src/jvarkit/dist/bioalcidae.jar -f file.js -F BAM

Answer 2

[John]

I hate competing with Pierre because I always lose :P but SeQC can do this so I ought to plug it.

Download all the code from https://github.com/JohnLonginotto/SeQC You want to use the TYPE module, which will give a number for all the possible orientations a single or paired end read can align, as a number from 1 to 20. To decode this number, use this chart: http://i.imgur.com/DCMrtwe.jpg

If the image quality isn't good enough or you want help running SeQC, i'll be sitting in ac.gt/chat for the next 4 hours to help