how to count the number of reads align per viruses in my sample

Question

[viviane20]

Hello,

I would like to count the number of reads aligning per viruses in my sample based on the output of minimap2 in a bam file.

Does anyone how to proceed?

By using the command: samtools view -f 2 ech22ReadsViralSorted.bam, I am able to see the alignments, but I am unsure about how to process the file.

Attached is the link to download the sample ech22ReadsViralSorted.bam

Thank you!

[viviane20]

All virus sequences are in a single file, which means that in the minimap2 output, I need to know the aligned viruses and the number of reads aligned for each virus

[GenoMax]

All virus sequences are in a single file

If the reference is a multi-fasta file in following format

> virus_1
Sequence
> virus_2
sequence
....
>virus_N
sequence

then run the command I posted below. It will show you all reference sequence names and reads aligned to each.

Answer 1

[GenoMax]

count the number of reads aligning per viruses

Assuming the virii were in the multi-fasta reference as independent entries you should use

samtools idxstats ech22ReadsViralSorted.bam

[viviane20]

yes the command line samtools idxstats ech22ReadsViralSorted.bam worked. Thanks you

[GenoMax]

Please accept the answer (green checkmark) to provide closure to this thread.

Answer 2

[Pierre Lindenbaum]

$ samtools view

(...)
  -c, --count                Print only the count of matching records
(....)