Counting number of short-read cover positions

[d_t_nguyen]

Hi everyone,
I used this command to count the number of short-read cover a position:

samtools view data.bam 20:1000-1000 -c

If I run this command several time, it would be very slow.
Is there any way to count the number of short-read cover millions of positions.
Thank you in advance!

[Rocketknight]

GATK's DepthOfCoverage can do this. Use the -L argument to specify which positions you want coverage information for. http://www.broadinstitute.org/gatk/g...fCoverage.html

[asiangg]

Hi, if you purpose is to get the coverage and visualize them, an easy to use and powerful tool is ngs.plot. It can accept a BED file or use its pre-compiled database and calculate the coverage for any bam file. It then visualize them either as an averaged profile plot or a heatmap.

[shi]

Hi d_t_nguyen,

You may try the featureCounts program included in the Subread package (http://subread.sourceforge.net/).

Here is the steps to use featureCounts to count your reads. Firstly, create a tab-delimited annotation file like this:

GeneID Chr Start End Strand
1 chr20 1000 1000 +
2 chr20 2000 2000 +
…

Then, issue the following command to count reads:

featureCounts -F SAF -f -O -a your_annotation.txt -i your_sam_file.sam -o results.txt

The counting results will be saved into file 'results.txt'.

Alternatively, you may provide a GTF format annotation for summarization. You may also use -T argument to run on multiple threads. Type 'featureCounts' for more information, or have a look at the userguide (http://bioinf.wehi.edu.au/subread-pa...UsersGuide.pdf)

Best wishes,

Wei

[d_t_nguyen]

Thank you all. I've just realized that pileup format of SAMTools have all I need. I used this command and I had all I want:
samtools mpileup individual.bam > read.txt