This is a test version of Biostars. For the public version, visit https://www.biostars.org.
count number of references with at least one hit

Hello,

In a context of mapping reads with multireferences, I was wondering how to count number of references with at least one hit (and not number of reads matching with ref).

A way would be to make a custom script with files obtain from samtools idxstats after alignment. Is there a samtools command or a software able to do that?

Thanks a lot

alignment references
• 2,136 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by vmicrobio vmicrobio

not clear: how the output would be different from samtools idxstats ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I just would like to get the information : on X references Y are mapped by at least one read. Do you think I should go with idxstats output?

• 0 views • link
written by vmicrobio vmicrobio

on X references Y are mapped by at least one read

ok, I understand now

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

1 answer

using bioalcidaejdk http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

we test the reads don't have : MAQ=0 || contains 'SA' tag || contains 'XA' tag

java -jar dist/bioalcidaejdk.jar -e 'stream().filter(R->!(R.getReadUnmappedFlag() || R.isSecondaryOrSupplementary() || R.getDuplicateReadFlag() || R.getMappingQuality()<=0 || R.hasAttribute("XA") || R.hasAttribute("SA"))).map(R->R.getContig()).collect(Collectors.groupingBy(Function.identity(), Collectors.counting())).forEach((K,V)->println(K+"\t"+V));' src/test/resources/S1.bam

RF02    290
RF03    280
RF11    72
RF01    358
RF10    82
RF08    114
RF09    114
RF06    146
RF07    116
RF04    256
RF05    170
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I had to format a little bit the heads of multifasta file but it works. Thanks a lot!

• 0 views • link
written by vmicrobio vmicrobio

and wc -l output file

• 0 views • link
written by vmicrobio vmicrobio

Log in to answer this question.

More posts like this