Is there a tool that will give detailed information on mapped reads beyond basic numbers? I'm looking for breakdown of distances of supplementary read mappings of long reads versus primary. Like how many within 200 bp, 500bp, 1000bp 3000bp etc
•
written
by
rob234king
using bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html and the SA:Z tag:
final Counter<String> contigs = new Counter<>();
final Counter<String> distances = new Counter<>();
stream().
filter(R->!(R.isSecondaryOrSupplementary() || R.getReadUnmappedFlag())).
forEach(R->{
for(SAMRecord R2: htsjdk.samtools.SAMUtils.getOtherCanonicalAlignments(R)) {
if(R2.contigsMatch(R)) {
int distance = Math.abs((int)(R2.getStart() - R.getStart())/1000)*1000;
distances.incr(String.valueOf(distance));
}
else
{
contigs.incr(R.getContig()+"->"+R2.getContig());
}
}
});
println("#contigs");
for(String c:contigs.keySet()) {
println(c+"\t"+contigs.count(c));
}
println("#distances");
for(String d:distances.keySet()) {
println(d+"\t"+distances.count(d));
}
invoke:
$ java -jar dist/bioalcidaejdk.jar -f input.code src/test/resources/retrocopy01.bwa.bam
#contigs
chr4->chr5 1
chr4->chr6 1
chr4->chr13 4
chr4->chr10 1
chr4->chr7 1
chr4->chrX 2
chr4->chr11 1
chr4->chrY 1
chr4->chr20 1
chr4->chr18 1
chr4->chr1 1
chr4->chr2 1
#distances
0 22
6000 13
5000 11
12000 13
1000 25
10000 8
2000 33
Thanks. it works, I feel like I should be able to unpick and work it via awk but why fight something that does the job
Just checking in the distances output 0 = reads from 0-1000 or just 0 1000 = reads from 1000-2000 or 1-1000?
0 = reads from 0-1000
thanks, how easy is it using the same tool to count how many alignments a read has, so from 1-10 alignments and how many alignments reads have? e.g. 1 alignment only - 1000 reads, 2 alignments - 2000 reads. I think that is different to what I have calculated here but applicable to understanding the context of those numbers.
Log in to answer this question.
