This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to collapse samtools depth generated bed file to regions with same coverage

I am creating per base depth files samtools depth. I want to collapse to a proper BED file for regions with similar depth:

Chr_01  140     10
Chr_01  141     10
Chr_01  142     10
Chr_01  143     30
Chr_01  144     30
Chr_01  145     30
Chr_01  146     30
Chr_01  147     20
Chr_01  148     20
Chr_01  149     20
Chr_01  150     20
Chr_01  151     20
Chr_01  152     10
Chr_01  153     10
Chr_01  154     10

Which should become (with trivial chosen names):

Chr_01  140     142    normal 10
Chr_01  143     146    high 30
Chr_01  147     151    medium  20
Chr_01  152     154    normal 10

Is this possible with bedtools or other existing tools?

samtools bedtools
• 4,104 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by BalderNGS • BalderNGS

2 answers 

$ cut -f 3 input.txt | sort| uniq | while read X; do awk -v X=$X '($3==X) { printf("%s\t%d\t%d\n",$1,$2,int($2)+1);}' input.txt | sort -k1,1 -k2,2n | bedtools merge -i - | sed "s/\$/\t${X}/" ; done

Chr_01  140 143 10
Chr_01  152 155 10
Chr_01  147 152 20
Chr_01  143 147 30

as it's BED, there is a +1 in the 'end' position.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks Pierre, I knew a single commandline would be out there. Thanks for the +1 on the end of the 'feature'. i missed that

• 0 views • link
written by BalderNGS • BalderNGS

Tao Liu (one of the authors of MACS) has a script for this in his Git. Just double check if his script assumes 1-based or 0-based coordinates. You might need to add +1 to the coordinates of the output. Anyway, usage is ./script.pl samtools_depth.txt > your.output

• 0 views • link
written by ATpoint ATpoint

Thanks, I will checkout which of the solutions is faster on large files.

• 0 views • link
written by BalderNGS • BalderNGS

Log in to answer this question.

More posts like this