We have released Goby 1.9.8.2. This version offers the vcf-subset and vcf-compare replacements tools I mentioned in my earlier VCF post.

The release also packs an option to call indels with Goby. We use the method of Krawitz et al (Bioinformatics 2010) to find equivalent indel regions (EIR). This approach can re-conciliate distinct indel observations into canonical  indel boundaries (an EIR). The genotype and compare-groups formats of the discover-sequence-variants mode will output EIRs at a frequency that sum over all the possible indel variations observed at the site that can be explained by that EIR. Of course, there is quite more to the Goby indel calling approach than the Krawitz method. For instance, the approach is integrated with the fast algorithm for local realignment around indels, so that indels that open when realigning end of reads contribute to the frequency of an EIR.

Programmers will find that Goby represents  observed indels at a site in a very similar way to base genotypes. Reading a base or indel frequency at a position in a sample is done with the same API (see the SampleCountInfo class). This makes it easy to support indels in different output formats.

The vcf-compare replacement (new in this release) can keep random samples of positions that differ between input files according to each category of differences it tallies (e.g., missed one allele RA vs RR, missed two alleles AA vs RR, genotypes differ C/T vs A/T where R=G). This is quite useful in inspecting positions in a genome viewer to try and understand differences between calls made by two approaches.

More details about this release are in the ChangeLog.

Stumbled on PLINK/SEQ while looking for a tool to estimate Ti/Tv from VCF files. At the moment, PLINK/SEQ seems limited to some older version of VCF, so it does not quite work with the files Goby generate (4.1), but I am interested in the mention of a binary file alternative to VCF, which could speed up some of the work we do (see section about project creation).

PLINK/SEQ genetics library

Tutorial: working with 1000 Genomes Pilot 3 VCFs. By way of introducing some of the features and approaches of PLINK/Seq, this page provides a tutorial that uses PSEQ and the R interface to PLINK/Seq …

Here's a nice blog-review about studies that discovered variations causing diseases with next-generation sequencing.

In addition to the content, I like the comment of the author that there is no time to write a review paper given the speed of the field development. Obviously what he means (since he has written the material already), is that publishing in a journal has very high overheads for the author(s) that simply slow down communication of information. In this case, I would argue that the comments have served a similar role as peer-review, prompting the author to add links to the papers, adding to the information or modulating the claims of novelty (see comment about RET).

We have recently started a large-scale evaluation of the genotype calling features of Goby and GobyWeb. To this end, we decided to obtain exome data from the 1000 genome project, and compare the genotypes called by Goby when all processing is done with GobyWeb (alignment and genotype calls). Since the 1000g project has way more data than we need for this evaluation, we picked two exome samples semi-randomly. Both are paired-end and one has length 76bp, while the other is 90bp long.

Exome data realignments

Realigning the reads to the 1000g reference was no trouble, we simply converted the bam files distributed by the 1000g project to the compact-reads format and uploaded this to GobyWeb. The rest is pretty much automated and was done in a matter of hours.

Extracting a few samples from the 1000g VCF files

The 1000g genome project distributes many versions of the genotype calls, in the VCF format. Locating the version that was produced against the 1000g reference (based on hg19) that we have installed in GobyWeb was a bit tricky since there is really no summary of all the versions. Thanks to Juan Rodriguez-Flores (in Jason Mesei’s lab at the ICB), for recommending this version:

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/working/20111111_old_phase1_release_files/

(we had initially gone directly to the very latest release, being wary of the “old” keyword in this version, but that the most recent version turned out to have been aligned against some ancestral reference reconstructed from primates, as far as we could tell and would not work for this validation).

As one can see form this directory, genotype calls are given in the VCF format, and split in one file per chromosome (excluding the X and Y chromosome and MT). The files are in the tens of GB even though they are compressed. The files are large because they contain genotypes and annotations for hundreds of samples studied in the 1000g genome.

Trying vcf-compare against just one of these files convinced us the comparison would be too slow against the complete files. We decided to extract just the two samples we selected for validation to yield smaller files that could be compared more efficiently.

Fortunately—we initially thought—VCF-tools provides a program called vcf-subset. Let’s just run this program with the names of the two samples we need to extract, on each of these chromosome, then concat the result. It turns out that vcf-subset is incredibly slow for the work it needs to perform. To be more specific, on a fast server, after a day of processing, we had not finished extracting the two samples from the chromosome 1 file. Upon closer inspection, the perl process was running at 100% CPU, but did not appear to make much process through this file (as judged from the speed at which results were added to the output). At this point, rather than throwing more CPUs at the problem and go on vacation (Chrismas time), we decided to go green (consider that inefficient programs are just as detrimental to our environment than other wasteful ways to burn oil).

Since Goby provides an efficient VCF parser, we reasoned we could write a more efficient way to extract the data we needed without too much trouble. To this end, we added a vcf-subset mode in Goby (an early implementation of this mode made it to 1.9.8.1, but we suggest getting the source code directly from our subversion server until we push 1.9.8.2 since the mode has improved a lot since that first release). Re-implementing the mode indeed provided a performance boost, but also offered an opportunity to add new options. One new feature that we added quickly was to process a number of files in parallel. We are now able to subset the 1000g VCF files in a few hours on a multi-threaded server.

Why is vcf-subset so slow?

This is obviously much better, but one has to wonder what is taking so long? After all, the input files are only a few Gigabytes, and we don’t need to do anything complicated, just extract a subset of information. It turns out that the design of the VCF format makes the task very computationally demanding (much more so than it would need to be). First of all, VCF is a text-based format, which by definition is slow to parse. To complicate matters, the format of the file can vary from line to line (look at the specification of the FORMAT and sample column). In my opinion, this is a very poor design decision. Consider  whether you have yet to encountered a VCF file that use the feature (e.g., that has different fields in the FORMAT field on different lines)? The need is not common, yet every program must be written to support this “flexibility” and there is a clear computing cost for supporting the feature. This is a typical red flag that should have told the designers of the added flexibility was not worth the compute cost.

Another complication introduced by VCF is that the type of delimiters varies by field (the INFO fields are delimited by ;, FORMAT fields use :, while other fields are delimited by tabs). All these “features” may seem to provide flexibility, but they combine to create significant inefficiencies. The format looks like as if it was designed so that humans can read it, yet is now used to store gigabytes of data. All this suggests to me that the committee/members of the mailing list who have been  responsible for the design of VCF could have paid more attention to the actual uses of the format and should have weighted  the impact of adding “nice to have, but not that frequent” features against the computational cost of these features. Then again, most people don’t care if a program or format is inefficient, at least until the computational cost makes some needed tasks impractical. It seems that VCF may be ripe for a redesign, it clearly meets requirements, but it is so complicated and inefficient that it would make sense to replace it with a leaner and efficient alternative.

I would not be surprised if the techniques we used in Goby formats could yield a more compressed VCF format alternative that could be subset in minutes rather than hours. Whether the community is feeling enough pain to consider adopting an alternative is a different question. What do you think?

We have released critical performance enhancements and bug fixes in Goby 1.9.8.1. A detailed list of changes can be found in the Change Log. All users are encouraged to install this latest distribution. Highlight include much better performance merging alignments with large >100MB .tmh files (as required when calling variations across many samples) and several fixes for subtle bugs that extensive testing has recently uncovered.

We have started to post GobyWeb training videos. Video demonstrate key features of GobyWeb and are designed to help users who cannot attend our regular training sessions at the Weill Cornell Medical College. Additional videos will be posted as we finalize the first GobyWeb release.

I have recently been interested in reanalyzing datasets from the study of Ajay et al, Genome Research (‘Accurate and comprehensive sequencing of personal genomes’). The data were deposited to the ENA, and the accession code given in the paper (ERP000765), so all seemed good.

I prefer to start with FastQ files, but these were not available for this study. I located BAM files at ENA under the Submitted files tab.

Goby provides a tool to extract reads from BAM files, so I did not expect major difficulties. For single-end reads, I would typically do:

goby 1g sam-extract-reads input.bam -o output.compact-reads

and proceed to analyze the compact-reads files with GobyWeb. The first problem I experienced was that the study used a paired-end design. This is a great idea when you are sequencing a complete genome, but we did not have an option in sam-extract-reads to load two BAM files and put the read and its mate together into one compact reads file. We do have this option in the tool fasta-to-compact, but since the files are tens of Gigabytes each, I did not feel like creating several huge intermediary files, or waiting for the conversions to run ((from BAM to compact-reads, then to fastq) x 2, then to compact again).

Using the options of fasta-to-compact as a template, I added the same capability to sam-extract-reads (this option will be available in our next release).

Here’s what you can now do:

goby 1g sam-extract-reads SRA_GAIIX_RUN_1.bam SRA_GAIIX_RUN_2.bam --paired-end -o output.compact-reads

This would cause the two input files to be loaded and each read combined with its pair and written to output.compact-reads (Goby stores primary and paired read in same file so that we can be sure these sequences never get out of order).

In implementing these options, I was careful to check that the id of the reads matched for the primary read and its paired read. As soon as I started this new tool on the datasets I obtained from ENA, I got an error indicating that the id of the first read in SRA_GAIIX_RUN_1.bam did not match the id of first read in SRA_GAIIX_RUN_2.bam.

I guess I should have expected these BAM files to be sorted in chromosomal order. Fortunately, I saw that samtools sort has an option to resort a BAM file by read id. My plan was to resort each file by id, then use sam-extract-reads.

Here’s the samtools command to resort by id:

samtools sort -n SRA_GAIIX_RUN_1.bam SRA_GAIIX_RUN_namesorted_1.bam

Since the GAIIX file was 23GB, I started the sort overnight. Here’s the message I got this morning:

[bam_sort_core] merging from 316 files…

open: Too many open files

[bam_merge_core] fail to open file SRA_GAIIX_RUN-namesort_1.bam.0252.bam

This does not look good. Let’s look at what happened here. The program samtools sort handles alignments that are larger than can fit in memory by splitting the input in small pieces, each large enough to fit in memory, then sorts the piece, and finally writes the result to disk. The last step of the process is to sort all pieces to write the final sorted output. The problem I encountered was that the default size for each piece resulted in so many files that the merge step failed (most operating systems have a limit on the number of files a process can open at any one time).

The work-around would consist in increasing the default -m parameter to a few GB and run on a machine with sufficient memory. This will reduce the number of files that need to be kept open in the merge step. However, this means that the whole process is now limited in the amount of memory one can get on a server. Consider that server memory is not growing nearly as fast as sequencing capacity and one can envision a point when the dataset is too large to process on any server. This is unlikely to happen tomorrow however, and the biggest problem is that one has to go through several steps (that do not run in parallel and hence are quite slow on such large files) when trying to reanalyze large paired-end datasets provided in BAM files.

It would have been so much simpler if the reads had been deposited in a format that keeps read pairs together. Goby compact reads does this nicely by keeping everything in one file, FASTQ does the trick with two files. Considering this experience, BAM files sorted by genomic location should probably not be the first choice when considering how to distribute full-genome paired-end NGS read data.

We have released Goby 1.9.7.3 [changelog]. New features include the ability to enter a short-cut for a mode name. You can call goby 1g ftc, where ftc is a short-cut for fasta-to-compact. See goby 1g ‐‐help for the list of all short-cuts. Another useful feature is the ability to convert sorted BAM files to sorted Goby alignments with sam-to-compact (short-cut: stc). This feature was suggested by Bradford Powell. A few bug fixes have been implemented, especially for VCF output. Finally, a bug in tally-reads was fixed that caused the tool to fail when processing reads with different lengths (thanks Adrian Platt for the bug report).