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
2 answers
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
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
Log in to answer this question.