SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
BWA and gaps efoss Bioinformatics 5 04-15-2013 07:42 PM
How to find coverage for RNA-seq, solid data? pinki999 Bioinformatics 2 09-28-2011 04:59 AM
sequence coverage profile/gaps for various sequencing machines mjpablo23 Illumina/Solexa 1 04-08-2011 01:33 AM
How many gaps??? saima Bioinformatics 6 04-28-2010 06:33 PM
Bowtie with gaps jjk Bioinformatics 2 12-14-2009 08:33 PM

Reply
 
Thread Tools
Old 01-29-2011, 01:34 PM   #1
Scotch
Junior Member
 
Location: Maryland

Join Date: Aug 2010
Posts: 6
Default Software to find gaps in coverage

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.
Scotch is offline   Reply With Quote
Old 01-29-2011, 11:12 PM   #2
simonandrews
Simon Andrews
 
Location: Babraham Inst, Cambridge, UK

Join Date: May 2009
Posts: 871
Default

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
simonandrews is offline   Reply With Quote
Old 02-09-2011, 05:27 AM   #3
msincan
Member
 
Location: Maryland

Join Date: Dec 2009
Posts: 19
Default

I have a python program that can do that. Please contact me at sincanm@mail.nih.gov 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.
msincan is offline   Reply With Quote
Old 02-09-2011, 05:51 AM   #4
quinlana
Senior Member
 
Location: Charlottesville

Join Date: Sep 2008
Posts: 119
Default

Quote:
Originally Posted by Scotch View Post
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:

Code:
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/bedto...usage-examples
quinlana is offline   Reply With Quote
Old 02-09-2011, 06:32 AM   #5
Adjuvant
Member
 
Location: Chicago, IL

Join Date: Sep 2010
Posts: 13
Default

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

Code:
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.

Code:
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:
Code:
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.
Attached Files
File Type: pl get_fragment_coverage.pl (7.2 KB, 43 views)
Adjuvant is offline   Reply With Quote
Old 02-09-2011, 11:54 AM   #6
msincan
Member
 
Location: Maryland

Join Date: Dec 2009
Posts: 19
Default

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.
msincan is offline   Reply With Quote
Old 02-09-2011, 11:59 AM   #7
quinlana
Senior Member
 
Location: Charlottesville

Join Date: Sep 2008
Posts: 119
Default

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.

Code:
samtools view -H [BAM] | grep "@SQ" | sed -e "s/SN://" -e "s/LN://" | cut -f 2,3 > chroms.txt
quinlana is offline   Reply With Quote
Old 02-09-2011, 12:06 PM   #8
msincan
Member
 
Location: Maryland

Join Date: Dec 2009
Posts: 19
Default

Quote:
Originally Posted by quinlana View Post
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.

Code:
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?

Last edited by msincan; 02-09-2011 at 12:23 PM.
msincan is offline   Reply With Quote
Old 02-09-2011, 11:15 PM   #9
bpetersen
Member
 
Location: Germany

Join Date: Mar 2010
Posts: 20
Default

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.
bpetersen is offline   Reply With Quote
Old 02-09-2011, 11:50 PM   #10
Bruins
Member
 
Location: Groningen

Join Date: Feb 2010
Posts: 78
Default

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.
Bruins is offline   Reply With Quote
Old 02-10-2011, 04:48 AM   #11
quinlana
Senior Member
 
Location: Charlottesville

Join Date: Sep 2008
Posts: 119
Default

Quote:
Originally Posted by msincan View Post
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.
quinlana is offline   Reply With Quote
Old 03-28-2011, 04:54 AM   #12
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 415
Default

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
colindaven is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:37 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO