This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Number Of Reads In A Sam File Which Are Assigned To Each Chromosome?

Given a sam file which was produced by using bwa to align a fastq of reads to a reference genome, how can I enumerate the number of reads which were assigned to each chromosome of the reference genome?

alignment sam

6 answers

Your question refers to SAM format, but in the event that you may also have a sorted/indexed BAM, you can simply do the following to get this very quickly.

samtools idxstats aln.bam | cut -f 1,3

The solution based on samtools idxstats aln.bam should be used with caution. This is because AFAIK the numbers reported by samtools idxstats (& flagstat) represent the number of alignments of reads that are mapped to chromosomes, not the (non-redundant) number of reads, as stated in the documentation. I stumbled across this by observing that the total numbers reported by samtools idxstats (& flagstat) exceeded the total number of reads in my input fastq files. See further discussion on this issue in the comments on this post: http://left.subtree.org/2012/04/13/counting-the-number-of-reads-in-a-bam-file/

Incidentally, I wrote a program that calculates read coverage from a sam file (number of reads mapped to each sequence, average depth, median depth, standard deviation, fraction of bases covered, etc).

pileup.sh in=mapped.sam out=coverage.txt

It has a toggle for including or ignoring secondary alignments; secondary=false or secondary=true (default).

When generating your SAM file make sure it has the header section before using pileup.sh.

samtools view -h file.bam > file.sam

Example pipe with filtering:

cat file.bam | samtools view -h -F 0x4 -F 0x100 | samtools idxstats -  |  cut -f 1,3

Brian, I'd be interested to read the source code of your program, but your link is broken. Could I ask you to provide a new link?

based on the first column only:

awk ' $1 !~ /@/ {print $1}' toy.sam | uniq -c

if you want to count only aligned reads you can also use awk to count only aligned reads.

Won't $1 denote the query name? $3 is the chrom

Here's how I count the number per sequence of primary alignments with proper read pairs:

Now there is a more convenient way with seqkit bam:

seqkit bam -C mybam.bam

This is the Samtools Solution + bash

Enjoy

PROCESSORS=10
BAM=blah.bam | blah.sam

samtools view --threads $PROCESSORS $BAM | cut -f 3 | sort | uniq -c | sort -nr | sed -e 's/^ *//;s/ /\t/' | awk 'OFS="\t" {print $2,$1}' | sort -n -k1,1 > Chr_freqs.txt

Log in to answer this question.