Better memory usage

We’ve made two changes to DISCOVAR de novo:

1. Memory usage during BAM file reading is now lower.

2. We’ve added a memory check to see how much memory appears to be available, and then throttle memory usage to that level. Sometimes (and particularly on scheduled systems like SGE) the amount of available memory is reduced. Please let us know if you observe aberrant behavior. The feature can be turned off by setting MEMORY_CHECK=False on the DiscoverExp command line.

Reference support added

We’ve added support for use of a reference sequence with DISCOVAR de novo. If you provide a reference sequence, then your assembly will be created de novo, then aligned to the reference sequence. Then when you view the assembly using NhoodInfo, the view will be marked to show reference coordinates.

To use this new feature, supply the command-line option REFHEAD=g, where you have files g.fasta and g.names. The file g.names should be a text file having the same number of lines as g.fasta has records, with each line being a name for the corresponding record (e.g. chr3). These names are displayed by NhoodInfo, so you want them to be reasonably short to avoid crowding the display.

Fractional read support added

We’ve added support for use of part of a read set. For example,

READS="frac:0.6 :: x.bam"

will cause 60% of the read pairs in x.bam to be chosen at random and assembled. This is useful to understand the effect of lower coverage and in cases where one has too much data to assemble on a given machine. The frac option can be combined with the sample option e.g.

READS="sample:T :: t.bam + sample:N,frac=0.5 :: n.bam"

to use all the reads in t.bam but only half of the reads in n.bam.

Memory limit added

DISCOVAR de novo now has an argument MAX_MEM_GB that can be used to limit memory usage to roughly the given amount. This can be useful on very large shared-memory systems.

New DISCOVAR de novo stats

We’ve added some new assembly statistics to DISCOVAR de novo. These are in the file stats in a.final and are mirrored in standard output. These along with the file frags.dist.png are often diagnostic.

Highest coverage paths now used in scaffolds

DISCOVAR de novo produces several output files, including a file of scaffolds a.lines.fasta in which a single path through a genomic locus is shown, even when multiple paths are possible (for one of several reasons, including polymorphism). (See “Edges, lines and scaffolds“.) This ‘flattened’ representation of the assembly loses information but has the advantage that it is FASTA and so can be processed by standard tools. With revision 51386, we now pick the paths used to be those having highest coverage. This is completely arbitrary in cases of bona fide polymorphism, but is helpful in cases where an assembly bubble occurs because of sequencing difficulty, making it uncertain which bubble branch is correct. In such cases, and in cases of ‘minor alleles’ in bacterial cultures, choosing the highest coverage branch makes sense.

Clarification of DISCOVAR input requirements

DISCOVAR and DISCOVAR de novo take as input read pairs from fragments of size 400-500 bp, with some larger and some smaller. The blog and manual contained references to fragments of size 700 bp, which were outdated, and have now been removed. Note that the protocol yields a wide size distribution, including some large fragments.

Understanding DISCOVAR output

A DISCOVAR de novo assembly is a graph. A typical assembly consists almost entirely of linear stretches, typically like this

simple

which we call ‘lines’, and providing a rich data type that captures polymorphism and other important features. Further, with some loss of information, these lines may be ‘flattened’ into standard contigs. We have added a tutorial explaining how these data types are available as part of the DISCOVAR output. We are also interested in hearing your thoughts regarding the utility of these output types and others that might be useful to you.

HiSeq 2500 data quality

We have been asked what our DISCOVAR input data looks like, and the best way to answer this question is with some examples. We don’t claim that these data are necessarily representative, but they do illustrate what we are able to generate here at the Broad Institute.