This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How get intervals of low coverage from a genome sequence?

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 low-coverage script
• 2,984 views • link
updated by Ram Ram • written by fhsantanna fhsantanna

What do you mean by 800pb ?

• 1 views • link
written by GouthamAtla GouthamAtla

1 answer

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.

• 1 views • link
updated by Ram Ram • written by GouthamAtla GouthamAtla

Very clever, thank you very much!

• 1 views • link
written by fhsantanna fhsantanna

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

• 1 views • link
written by fhsantanna fhsantanna

Log in to answer this question.

More posts like this
  • Generating artificial NGS reads from a FASTA file specifying read length, coverage, and start seque…
    written by Hamish • Hamish
    I am looking to generate artificial single-end Illumina sequencing data from a FASTA containing a number of short sequences representing unique molecular identifiers. However, I've …
  • Criteria for filtering contigs after spades assembly
    written by fhsantanna fhsantanna
    I have assembled a bacterial genome using spades. The coverage among the contigs is very heterogeneous, while some have more than 30 fold of coverage, …
  • How can I calculate coverage from the base frequencies file?
    written by SaltedPork SaltedPork
    Hi What I have is a tab delimited file with the frequency of each base at each position in the sequence. How can I take …
  • How to plot mutation hotspot according to mutation counts
    written by wangshx • wangshx
    There is a `data.frame` including the `number of mutation` & `start` position and `end` position of mutation interval, just like > head(chr1) start end mutation_count …
  • findOverlapsOfPeaks() of package ChIPpeakAnno in R
    written by salamandra salamandra
    When running the following code with ChIPpeakAnno package in R: gr3 = GRanges(seqnames = c("chr1"), strand = c("+","+","+"), ranges = IRanges(start = c(1,3,5), width = …
  • searching reads with a certain sequence in fastq file
    written by sumithrasank75 sumithrasank75
    I have a file say master.fastq which looks like : ``` @M00990:202:000000000-ADM27:1:1101:21678:1536 2:N:0:291 CCTTTTACCGACCCGCTCTTTCTCTCCTACGCTTATTTCCGTCTACCCTTCTCTTCACTCGCTATTTCTATTCTTAAAACTATCTTAATGTTCTGCCTTTGCTCTTTTCTTTTTTCTATAACCTCTCTACAGCCAACTCACCCATCTCCTTCCCTGCTACGCTATTCCTCTGTTAGTTTTTCTTCATCATACTTTTCTCATCTCACACTACCTTTGCACTTCTTCCTTTCCACGTCCCCTTTCTCCTACC + -----,,6;6,@+8++6+,,,<C5@,,,,,,+8,,;,,<,,,,,,,;,,;,,,<C,,,;,,8,,8CC,C,,,5,C,,,,99,+,+4,,,3,,9,,6,@<,,,,,,,9,,,,,,,4),,0**,,5),,5**,)59*0),,*)5,)))*,9,3,0++))5D:)))))5;+;+0)*;+*6++++),******3**50,6,+++**,0*,,31)88*0*1*5*1)0*:*7>C;3,035:0)))8).*2**.*:) @M00990:202:000000000-ADM27:1:1101:22685:1539 2:N:0:291 CTTATCACCGACTCTCTCCTTCTCTTCCAAGTTTATTTCCGACTCCCCTTATCTTCACTTGCTATTTCTATTCTTAAAACTATCTCGACCTTTCACCTTTCCCTCTTTCCTTCTTTTCTCTCCTTCTACACTCCCACCCACTCTTACTTCTTTCTTGTCACCGTTTCCATATTATACTTTCTTCTCTTACATAATTTTCTTCCTGCAAACTATTTAAGCAATCTCTTTCTTTCACCCCTTTTATCTCGCC + -----,,<;67@+B,,6,,,,;C5@,,,,,66<,,6,,<,,,66,,;,66,,;;C,,,;,4<,,<CC6C,;,,,;,,6,:,6?A9=,,+++2,5,,,9<?D,,,,,:C?,@,@,,5,?*3*,9**,,0**,*)93))4+))19*0**,,52,56+**5*03*)3)))42+2***5*+=3+,*4*2****,**2*,3,+++*0,50,,**5*****0****)5)***0**,,***3*)0)))3***0*))) ``` I want only …
  • Can these fusion predictions by tophat be trusted
    written by Dataminer Dataminer
    Hi, I have got these results from tophat-fusion: ``` chr15-chr17 74317267 38506035 ff 1 13 1 7 43 33 10.000000 @ 13 26 38 51 …
  • plink analysis including covariates and gxe
    written by ddzhangzz ddzhangzz
    <p>I would conduct a GWAS using plink with covariates and gxe terms. Here is part of phenotypic data:</p> <pre> FID IID luperml2wk18 luperml2chg rotaarm sex …
  • Interval Format Variant Files Convertion Problems And Igv And Annotation Application
    written by Tonyzeng • Tonyzeng
    <p>I produced .sam formal variants files by Mpileup and then converted them into .bam. however, after several steps for filtering variants, I got interval format …
  • Aligning Miseq With Deletions >10 Bases
    written by Leszek Leszek
    <p>We expect deletions (10-1000 bases) in the genome of interest, therefore, we run paired MiSeq 300bp+200bp. Can you recommend me the best aligner for such …