This is a test version of Biostars. For the public version, visit https://www.biostars.org.
how to count the number of reads align per viruses in my sample

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!

rnaseq samtools illumina
• 1,729 views • link
updated by GenoMax GenoMax • written by viviane20 • 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

• 0 views • link
written by viviane20 • viviane20

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.

• 0 views • link
written by GenoMax GenoMax

2 answers

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
• 0 views • link
written by GenoMax GenoMax

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

• 0 views • link
written by viviane20 • viviane20

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

• 0 views • link
written by GenoMax GenoMax 
$ samtools view

(...)
  -c, --count                Print only the count of matching records
(....)
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this