count number of references with at least one hit
1
0
Entering edit mode
5.0 years ago
vmicrobio ▴ 290

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 • 1.4k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

on X references Y are mapped by at least one read

ok, I understand now

ADD REPLY
2
Entering edit mode
5.0 years ago

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
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

and wc -l output file

ADD REPLY

Login before adding your answer.

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