I have two PCR amplicons that have been multiplexed and sequenced using the nanopore minion.
I have aligned the fastq reads using minimap2 with a reference file containing both amplicon sequences and generated a bam file that I have viewed using IGV.
I am looking for a way to generate some simple summary statistics.
In particular, is there a way to extract the total number of fastq reads aligning to each amplicon reference from the bam file?
The quick way to get the number of alignments on each reference is
samtools idxstats my_bam.bam
Number of reads on each reference is column 3. Although, as has been pointed out, this will give you the total number of alignments per reference, not the total number of reads (each read might give rise to more than one alignment). That said I do tend to us this as generally I'm after a rough approximation, rather than an accurate number.
In theory, only one alignment for each read should be marked as primary, so the following should give you what you need quickly and at low memory usage:
samtools view -bF 2304 my_bam.bam > primary_only.bam
samtools index primary_only.bam
samtools idxstats primary_only.bam
This one-liner below will work better for long reads than samtools flagstat in that it only counts the primary alignment for each read and samtools flagstat doesn't seem to calculate some stats for long reads. I have never seen samtools flagstat output stats on a per-reference basis, but am curious if so!
This answer filters out secondary and supplementary alignments for your reads (-F 2304) that have some alignment to both amplicon reference and just keep the best one. This might give a more accurate idea of how many reads of each amplicon are in the library.
samtools view -F 2304 myfile.bam | awk -F $'\t' '{a[$1, $3]++} END{for (i in a) {split (i, sep, SUBSEP); print sep[1], sep[2], a[i]}}' | uniq | awk '{print($2)}' | uniq -c | sort -k1 -nr
This can be done using
samtools flagstat your_bam_file
