Sort/create list of locations with highest # of reads?

[ronton]

Hello,

Hopefully there is a simple way to do this that I am unaware of.

I have fastq and bam files (100bp paired end Illumina). In the alignment files there are lots of areas with one or a few reads spread throughout the genome.

What I would like to do is either sort results to visualize, for example in Integrative Genomics Viewer, or generate a list that will show in descending order the locations that have the greatest number of reads or the greatest coverage. This would most likely be within a window, say 50bp.

For example:

chr21:33,031,597 has 58 reads
chr21:33,031,200 has 56 reads
chr21:33,025,000 has 42 reads

and so on.

Any advice or guidance on different ways I might achieve this is very appreciated.

Thank you

[Richard Finney]

Check out samtools ...

$ samtools depth

Usage: samtools depth [options] in1.bam [in2.bam [...]]
Options:
-b <bed> list of positions or regions
-f <list> list of input BAM filenames, one per line [null]
-l <int> minQLen
-q <int> base quality threshold
-Q <int> mapping quality threshold
-r <chr:from-to> region

Example:
$ samtools depth -r chr1:1-100000 1221792.bam | head
chr1 10033 1
chr1 10034 1
chr1 10035 1
chr1 10036 1
chr1 10037 1
chr1 10038 1
chr1 10039 1
chr1 10040 1
chr1 10041 1
chr1 10042 1


Use "awk" to filter on column 3 and "sort" on column 3 and "tail" it.
Like this ...
samtools depth 1221792.bam | awk '{if ($3>300) print $0}' | sort -k3n | tail

[ronton]

Thanks so much Richard,

So if I use:

samtools depth 1221792.bam | awk '{if ($3>50) print $0}' | sort -k3n | tail

for example, this will give me a list of all locations with more than 50 reads, correct?

But no 'window,' or in other words, it will likely print the individual locations for each block of reads, rather than one result for each bin or chunk of 50 bases or something like that. Meaning, at chr1 10038 +/- 25bp, there are max 50 reads.

Certainly a start, thank you for the help.

[swbarnes2]

samtools view -F 4 sorted.bam | cut -f 3,4 | uniq -c | sort -nr | head -n 20 > most_common.txt

That'll give you the top 20 most common start sites for your reads.