A: how to prepare the input file for the ReadDepth - http://www.biostars.org/p...
Feb 8, 2013
from
bamToBed would get you part of the way there, but there are two caveats: 1) You need per-chromosome bed files 2) readDepth actually only looks at the first and second columns. You can save yourself some disk by omitting the third column. Use the following command on your sorted bam file instead:
for $chr in $(seq 1 22) X Y;do samtools view myfile.bam $chr:1-999999999 | cut -f 3,4 >$chr.bed;done
If your bam uses the "chr" prefix ("chr1" instead of "1"), use this instead:
for $chr in $(seq 1 22) X Y;do samtools view myfile.bam chr$chr:1-999999999 | cut -f 3,4 >chr$chr.bed;done
Keep in mind that this will use a lot of disk space! Just make sure the files match your annotations. (either use "chr" or don't, but be consistent).
Sidenote: If I had to write this package over again, it would read bam/sam natively, but a lot of this work was done before bam/sam was in widespread use. I'm working on a new followup tool that should be out this summer.
- Chris Miller