Software to find gaps in coverage

[Scotch]

I'm trying to make a list of gaps in coverage in Illumina sequence data? IGV displays the number of reads covering each base, but is there a way to output this count data to a text file? This will save me having to manually scroll through the data in IGV and annotate by hand. Many thanks in advance.

[simonandrews]

You could do this in SeqMonk. The process would be:

1. Import your raw data
2. Do a contig probe generation with a depth of 1 read
3. Do a read count quantitation
4. Go back to probe generation and do an interstitial probe generation (so you put probes between the existing probes, which will be all the gaps in your assembly) - make sure not to exclude probes with no data!
5. Do a read count quantitation (even through they will all be 0)
6. Create an annotated probe report and you'll have a list of all of the gaps in your assembly
7. Alternatively you can convert your probe set into an annotation track so you can see where the gaps are when performing further analyses

[msincan]

I have a python program that can do that. Please contact me at [email protected] if you are interested in using it. It can do sums, averages and determine streches of coverages based on thresholds using aligned bam files as input.

[quinlana]

I'm trying to make a list of gaps in coverage in Illumina sequence data? IGV displays the number of reads covering each base, but is there a way to output this count data to a text file? This will save me having to manually scroll through the data in IGV and annotate by hand. Many thanks in advance.

genomeCoverageBed in the BEDTools package will do this for you.

The "-d" option will report depth at every base (perhaps excessive). The -bga option will create a BEDGRAPH file documenting intervals with common depth. You could use awk on the output to report solely those intervals with 0 or <X coverage. Something like this:

genomeCoverageBed -ibam aln.possorted.bam -bga | awk '$4 == 0' > intervals.with.0.cov.bedg

Check out the BEDTools manual or see the examples here:

http://groups.google.com/group/bedtools-discuss/web/community-usage-examples

[Adjuvant]

I also go the Bedtools route, but I wrote a quick perlscript to output some coverage statistics that can then be imported into Excel. I've been working with alignments against multiple fragments, but it works for single sequence alignments as well.

Here's my method:

Calculating alignment coverage

genomecoveragebed –ibam alignment_sorted.bam –g reference.fa.fai –bga > alignment.coverage

** Generates a file listing read coverage at each base for each genome. The "reference.fa.fai" file is the one created by samtools faidx.

perl path/to/get_fragment_coverage.pl –c alignment.coverage –o alignment_coverage.txt –d 1

** Generates a text file with the number of uncovered bases per genome or fragment and the total percent coverage of the entire genome or fragment. The –d flag is optional and designates how deep the coverage at a base should be for it to be counted as “covered.” The default is that if even a single read covers a base, it is counted as covered.

By the way, if you want you can pipe genomecoveragebed and get_fragment_coverage.pl together if you want, like so:

genomecoveragebed –ibam alignment_sorted.bam –g reference.fa.fai –bga | perl path/to/get_fragment_coverage.pl –c STDIN –o alignment_coverage.txt –d 1

And here's my script. Remember, I wrote it myself and it works pretty well for me, but feel free to adapt it for your purposes as needed.

[msincan]

I just tried bedtools genomeCoverage and it gave an error and exit. The error is Chromosome chr_random10 found in your bed file but not in your genome file . So my hg18.fai only has regular chromosomes and it is these type of issues that led me to write my own program. I agree that bedtools is very fast and efficient but it doesn't fit every requirement easily at least for me. Do you know any workaround to this (i.e a flag to tell bedtools to ignore such errors) other than creating another fai with random chromosome positions.

[quinlana]

I'm not sure why you wouldn't be using the same .fai as the header in your BAM, but in this case you could just make a new "genome" file from your BAM header.

samtools view -H [BAM] | grep "@SQ" | sed -e "s/SN://" -e "s/LN://" | cut -f 2,3 > chroms.txt

[msincan]

I'm not sure why you wouldn't be using the same .fai as the header in your BAM, but in this case you could just make a new "genome" file from your BAM header.

samtools view -H [BAM] | grep "@SQ" | sed -e "s/SN://" -e "s/LN://" | cut -f 2,3 > chroms.txt

It works, thanks, but how can I get coverage for only those regular (non-random) chromosomes?

[bpetersen]

At the risk of this being a stupid question: What exactly are these "random" chromosomes? I stumbled over this when converting a soap alignment to bam. I think alignments to these random chromosome s were just skipped in the conversion and I would like to know if I'm now missing important data in my bamfile.

[Bruins]

They are unassembled supercontigs. They can be found in hg19 (at least, in 1kg.37) but not in hg18.

Can someone confirm if I got this right?

Anyway, we do use them to align to, because there may be reads that should be aligned there instead of at the regular chromosomes. If the supercontigs are missing, those reads may be forced to align in a place where they don't belong. In short: reduce the number of false positives.

Sorry for being slightly off topic.

[quinlana]

It works, thanks, but how can I get coverage for only those regular (non-random) chromosomes?

Just add a grep -v "_random" after the output of genomeCoverageBed to filter any intervals arising from the "random" chroms.

[colindaven]

I couldn't get the following command to work:

genomeCoverageBed -ibam aln.possorted.bam -bga | awk '$4 == 0' > intervals.with.0.cov.bedg

The following worked fine:

samtools view -b mybam.bam | genomeCoverageBed -ibam stdin -g PA14genome.txt -bga | awk '$4 == 0' > intervals.with.0.cov.bedg