Question: How get intervals of low coverage from a genome sequence?
2
gravatar for fhsantanna
5.1 years ago by
fhsantanna470
Brazil
fhsantanna470 wrote:

I have a file with the following format, that relates a determined position of the genome to coverage data:

Sequence position coverage

NC_006510    1    19
NC_006510    2    22
NC_006510    3    44
NC_006510    4    1
NC_006510    5    1
NC_006510    6    0
NC_006510    7    0

... (low coverage positions, under 6)

NC_006510    1000    0
NC_006510    1001   66
NC_006510    1002    66
 

I would like to get intervals larger than 800 pb with low coverage along a sequence. For example, a desirable output would be like that:

NC_006510 4..1000 #In this position interval, the sequence has low coverage (<6).

Do you know a clever way to do it?

Thanks in advance.

 

interval script low coverage • 1.5k views
ADD COMMENTlink modified 5.1 years ago by geek_y10k • written 5.1 years ago by fhsantanna470

What do you mean by 800pb ?

ADD REPLYlink written 5.1 years ago by geek_y10k
6
gravatar for geek_y
5.1 years ago by
geek_y10k
Barcelona
geek_y10k wrote:

If you are looking at 800bp or more stretches with coverage less than or equal to 6X, some thing like following would be useful.

bedtools genomecov -ibam input_sorted.bam -g reference.fa.fai -bg | awk '{ if ($4<=6) print }' | bedtools merge -i - | awk '{ if (($3-$2)>=800 ) print }' | less

Anyway this is not an exact answer to your question, but different approach.

ADD COMMENTlink modified 5.1 years ago • written 5.1 years ago by geek_y10k

Very clever, thank you very much!

ADD REPLYlink written 5.1 years ago by fhsantanna470

-bga option is better because computes regions of coverage "0". Also in "bedtools merge" the parameter -i has to be written.

ADD REPLYlink written 5.1 years ago by fhsantanna470

updated.         

ADD REPLYlink written 5.1 years ago by geek_y10k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2220 users visited in the last hour