bam stats F1F2/F1R2/R1F2/R1R2/F2F1/F2R1/R2F1/R2R1
2
2
Entering edit mode
7.6 years ago

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

bam • 6.2k views
ADD COMMENT
2
Entering edit mode
7.6 years ago

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
ADD COMMENT
1
Entering edit mode
7.6 years ago
John 13k

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

ADD COMMENT

Login before adding your answer.

Traffic: 840 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