Simple tools for assembling and searching high-density picolitre pyrophosphate sequence data.

BACKGROUND: The advent of pyrophosphate sequencing makes large volumes of sequencing data available at a lower cost than previously possible. However, the short read lengths are difficult to assemble and the large dataset is difficult to handle. During the sequencing of a virus from the tsetse fly, Glossina pallidipes, we found the need for tools to search quickly a set of reads for near exact text matches.
METHODS: A set of tools is provided to search a large data set of pyrophosphate sequence reads under a "live" CD version of Linux on a standard PC that can be used by anyone without prior knowledge of Linux and without having to install a Linux setup on the computer. The tools permit short lengths of de novo assembly, checking of existing assembled sequences, selection and display of reads from the data set and gathering counts of sequences in the reads.
RESULTS: Demonstrations are given of the use of the tools to help with checking an assembly against the fragment data set; investigating homopolymer lengths, repeat regions and polymorphisms; and resolving inserted bases caused by incomplete chain extension.
CONCLUSION: The additional information contained in a pyrophosphate sequencing data set beyond a basic assembly is difficult to access due to a lack of tools. The set of simple tools presented here would allow anyone with basic computer skills and a standard PC to access this information.

Background

The introduction of micro-fabricated high-density picolitre reactor pyrophosphate sequencing [1,2] by the company 454 Life Sciences (454 Life Sciences Corp., 20 Commercial Street, Branford, Connecticut 06405, USA; hereafter referred to as 454 sequencing) makes available for the first time large quantities of sequence data at reasonable cost. The continual reduction in sequencing cost will encourage sequencing by small groups or individual researchers with modest computer resources and limited experience of bioinformatics tools.

The nature of the data from pyrophosphate sequencing is however both qualitatively and quantitatively different from that generated by Sanger sequencing [3] using fluorescent chain-terminating nucleotide analogues [4]. Instead of receiving a single consensus sequence with associated chromatogram (scf) file, this form of pyrophosphate sequencing generated for us 300 000 short sequence reads, around 100 bases long [5], assembled into several hundred contigs. However the normal system of checking dubious base calls on the chromatogram is not available, the qual and other associated files were not recognized by standard assembly tools (Phred/Phrap, DNAStar, NTiVector etc.) at the time we received the data (the ace file format has subsequently been made compatible with consed) and the very large number of short reads are difficult to assemble [6]. Assemblers designed for short reads are available [7-9] that utilize the de-Bruijn graph method [10] but they still suffer from shortcomings that result in the production of a large number of short contigs. Assembly is improved by the use of paired-end ditags [11], but this requires an additional sequencing run, thereby doubling the cost and does not necessarily resolve all the problems with repeat regions [12]. The additional cost though may be acceptable when compared to the cost of finishing by hand. The assemblers are also sensitive to read errors which generally occur at higher frequency in pyrophosphate sequencing than in Sanger sequencing [1] so the assembly protocol often incorporates an error correction step that corrects or removes sequences with low frequency [7,13]. This approach is acceptable when the source material is clonal and to produce a first consensus assembly, but when it is not clonal this error correction will mask polymorphisms in the source material. The pyrophosphate sequencing technique is still developing and both the average read length (now about 230 bases) and the number of reads returned per run (now around 500 000) will continue to increase, making this technique more attractive and the need for suitable tools greater.

Recently a novel DNA virus from the tsetse fly Glossina pallidipes (Diptera: Glossinidae) [14,15] (salivary gland hypertrophy virus, GpSGHV) has been sequenced by 454 sequencing; at this time paired-end ditag reads were not offered. The source material was not clonal, having been purified directly from infected salivary glands. The purified DNA was divided into two and submitted to 454 as two independent samples run together on a single reaction slide. After sequencing and assembling 454 were informed that the samples were from the same source and they were requested to assemble the combined fragment databases. The virus has a size of about 190 000 bases and the sequencing from the combined fragment database generated 321 953 reads averaging 108.8 bases each for a total of 35 023 571 bases read, giving about 185-fold coverage. These reads were assembled into 402 contigs by the Newbler software at 454. Prior to the pyrophosphate sequencing a number of clones had been sequenced by Sanger sequencing and subsequently targeted PCR amplification and sequencing were used for gap filling. However, it was realised that there was considerable additional information in the fragment database supplied by 454 if it could be suitably analysed.

Most individuals or small groups tempted by the low cost of pyrophosphate sequencing will probably use Windows based PCs. However, after the initial assembling, finishing the sequence and extracting additional information from the large database of individual reads presents a challenge when using Windows based programs. This is because the large datasets are too big to manipulate easily with standard Windows programs. There is, therefore, a clear need for assembly tools for this type of data that will permit detailed interrogation of the full dataset, but in the meantime we offer some simple tools that can be used to help to manually correct the problems encountered with the existing assembly programs and to investigate repeats and polymorphisms. We also discuss a number of issues related to the pyrophosphate sequencing revealed by examining our dataset.

Methods

The simple tools described here provided the ability to search a database of more than 300 000 fragments for matches, collect the fragments, sort and count them, count the bases adjacent to the match and to perform a simple assembly and associated functions. To produce these tools easily and to make them open source they are based on ISO standard C code compiled under GNU Linux and script files running standard GNU Linux utilities. To make them readily accessible to users of Windows systems unfamiliar with GNU Linux instructions are provided for running the tools using the Knoppix "live" CD [16]. This allows Debian GNU Linux to be run on any PC with a bootable CD or DVD drive without making any changes to the Windows installation. The Knoppix GNU Linux does not use the PC hard disk drive, leaving the original Windows system completely unchanged. The tools have been tested under Knoppix 5.1, Red Hat 9 and Fedora 7 on several 32-bit PCs.

The current native disk format from Windows 2000 (NTFS) cannot be reliably written by some versions of Linux so we recommend a USB flash memory drive in order to have storage available that is accessible both to the Linux system and to Windows. An external hard disk with a FAT32 partition or a FAT32 partition on the internal hard disk can also be used. Instructions for setting up the necessary file structure and running the tools that should allow anyone reasonably familiar with a Windows system to run the Knoppix system and use these tools is provided in Additional file 1. The supplementary material consists of three files, the instructions [see Additonal file 1], the sample files, source code and compiled objects for the tools [see Additonal file 2], and the supplementary tables [see Additional file 3]. The latest versions are available from our web page [17].

The tools

Performance figures are given for the tools run under Knoppix DVD 5.1 on an Intel Pentium 4, 1.5 GHz, with 256 kB cache and 512 MB RAM with the dataset already loaded into memory. Initial loading of the database takes about 30 seconds, and performance is dramatically reduced if the dataset needs to be reloaded into memory. As the tools all run in memory, disk or flash memory access speed is not an issue, affecting only the first loading of the database so long as the memory size is large enough to accommodate the dataset and any programs being run. For this reason, a minimum of 384 MB RAM is recommended. Unless otherwise stated the performance figures are based on the search string ACTGCTAAGTAATTTGTGAA, which produces 102 matching fragments with our full dataset of 331 953 fragments. The time taken for the tools is linear with respect to database size and to the number of matches produced.

The tools consist of the following.

fna_to_fnb output_file

fna_to_fnb converts the data file (fna) supplied by 454 into a format usable by these tools. The fna file is in fasta format and is converted by removing the line break between each header and sequence so that they are on one line to make searching easier and quicker. The header lines are converted to lower case, to distinguish them from the following sequence, which is converted to uppercase. See the supplementary material pyrophosphate_readme.pdf for further details.

build primer_file output_file

This routine takes a primer string (minimum 60 bases, maximum 250 bases) from the file primer_file and searches the database for exact matches to a 60 base string, forward and complement. The consensus base for position 61 is determined from the collected fragments, it is then added to the primer string, and the 60 base match frame moved forward one place. After 60 bases the result is written to output_file as a single line, with a progress message written to standard error. If the number of matching fragments falls below 8, or the proportion that the fragments containing the consensus base to all matches of the 60 base primer string is less than 70% (as at the end of a repeated sequence where the sequence becomes unique again) the program stops, writing the remaining data out and displaying the problem. The problem can then be further investigated using other tools. If the input primer string is more than 60 bases, the bases up to the last 60 are written to standard out before starting the search. The tool stops automatically after 1920 bases if not stopped before this by a mismatch, in order to avoid being caught in an endless loop.

This tool searches only for exact matches to the 60 base string in the primer. As the dataset is so large, at each position there will typically be 20 – 40 matching fragments. When the dataset is smaller it will fail more quickly, and as the average number of matches approaches the limit of 8 it will not work. Exact repeats of more than 60 bases will result in a loop being formed that has to be resolved with other tools. 1920 bases are assembled in about 23 minutes.

re [match_length] input_file output_file

re takes an existing draft sequence and checks that the draft is supported by the fragment dataset at each base position. There are two mandatory arguments; the input and output file names. The tool ignores all characters in the input file except uppercase A, C, G and T, so it will accept a standard FASTA format file. The output file is in fixed width text format, and can be loaded into other (Linux or Windows) programs, such as Excel, for further processing. The first parameter is optional; if present it specifies the match length of the search, if not present a default length of 65 bases is used.

re takes match_length bases of the input_file and counts the number of reads containing match_string with the immediately preceding base and the number containing match_string plus the immediately succeeding base. The output is in the form of a table (Supplementary Table 1; see Additional file 3). The third column is the match_string, the second column is the immediately preceding base and the first column is the position of this base in the sequence being re-built. The fourth column is the immediately succeeding base. Column 5 is the number of matches to match_string plus the preceding base, and column 6 is the percentage that column 5 represents of all matches to match_string with any preceding base (but not N). Columns 6 and 7 likewise present the data for match_string with the succeeding base. In the succeeding line the match_string is stepped on one base and the searches repeated. This output can easily be further processed in a spreadsheet to display the number of matches, with matches below a threshold highlighted (Fig. 1) indicating that further investigation is required. An abrupt drop to zero matches identifies a base error, an insert, deletion or substitution and a sudden increase in matches above the general level identifies an exact repeat region. re takes about 6 minutes 30 seconds to check a 4 100 base sequence with the default 65 base match.

Figure 1

Graphical representation of the output from re. The sequence of red points on the x-axis near 182400 represents a mismatch region where a single base was inserted by the assembler. The areas with many more reads than the average show identical repeats. The match length was set to 27 for this Fig.

fr match_string output_file

fr takes any length string from standard in and outputs all fragments with an exact match from the dataset to output_file. If output_file is omitted the output is sent to standard out. Reverse matches are also output, complemented to bring them into the same sense as the original match_string. The original match_string is included at the beginning of the output (Supplementary Table 2; see Additional file 3). fr takes about 1.3 seconds.

fnb_to_fas input_file >output_file

fnb_to_fas takes a file produced by fr and converts it into standard FASTA format for further processing in an assembler or other program, sending the output to standard out unless redirected. If no parameter is supplied it takes input from standard in (keyboard) until terminated by -d and outputs to standard out. fnb_to_fas takes about 30 seconds to convert 100 000 fragments to FASTA format.

fr2 match_string output_file

fr2 collects fragments as fr, but then truncates the beginning of each fragment to the first character of match_string and performs a simple alphabetic ascending sort to simplify visual comparison of the fragments. The forward and reverse matches are output separately, identified with > and < respectively and separated by one line (Supplementary Table 3; see Additional file 3). A related function, frx, produces the same output but with the match_string replaced by to make it easier to view longer sequences. fr2 takes about 1.3 seconds.

count match_string

count searches the dataset for exact matches to match_string and outputs two numbers, the forward and reverse match number. count takes about 1 second.

c4 match_string

c4 collects exact matches to match_string in the dataset and counts the number of occurrences of each of the four bases in the position immediately following match_string. The output is written to standard out and consists of twelve lines with match_string followed by each of the bases in the order A, C, G, T and the count of forward and reverse matches. c4 takes about 3 seconds.

c4p match_string

c4p collects the same information as c4 but with each of the four bases in the position immediately preceding match_string. The output is written to standard out and consists of twelve lines with match_string preceded by each of the bases in the order A, C, G, T and the count of forward and reverse matches. c4p takes about 3 seconds.

cmpl input_file

cmpl outputs the reverse complement of input_file to standard out. Only the uppercase characters A, C, G and T are complemented; all other characters are passed through unchanged, except for a limited set of wild-card characters (see Regular Expressions below) and > which is converted to <. If no input_file is specified cmpl waits for input from the keyboard and continues to take input from the keyboard until it is ended with -d.

A complete multi-line sequence can be reverse complemented by using cmpl followed by the Linux command tac, e.g.

cmpl input_file | tac >output_file

Note that any FASTA header will be mangled and will come at the end of its sequence. cmpl is intended to operate on files with a line length of 60 characters cmpl takes about 0.27 seconds to complement 4100 bases in lines of 60 characters; the time increases linearly with number of lines and as the square of the line length.

Where any of these tools expects input from standard in (the keyboard) or sends output to standard out (the monitor) the input or output can be redirected to or from a file with < and > respectively, e.g.

fr

cmpl input_file >output_file

Regular expressions

A very limited sub-set of regular expressions has been implemented to help with searching for partial matches and polymorphisms.

Wildcard.

The wildcard character. (period) can be used to substitute for any base. Two or more may appear in succession, e.g. AACT.TCC and AACT...TCC are both valid.

Multiplier *

The * character returns zero or more of the immediately preceding character, e.g. A* matches zero, one or more A, so AA* matches one or more A.

Alternatives [mn]

[mn] matches one of the characters contained in the square brackets, e.g. [AT] matches either A or T.

One or more of these may appear in any match_string or search_string, but they may not be combined in the same position; i.e. ACT [AT]CCG*A..A is valid but ACT [AT]*C is not as [mn] and * cannot be combined.

Results and Discussion

Additional information available from the 454 dataset

Although the 454 sequencing method produces characteristic errors as noted below, it also contains additional information that is not available from conventional sequencing. When the material to be sequenced is clonal, the result of the 454 method is essentially equivalent to Sanger sequencing. However when the original material is not clonal, the 454 sequencing will effectively sample the sequence many times at each position (approximately 185 times in our sample), permitting any polymorphisms to be investigated.

Single-base polymorphisms

If an alternate sequence exists in the sample being sequenced at the level of 10% or so, Sanger sequencing will always return the majority sequence as the minority sequence signal will be masked by the majority signal. The 454 system however effectively samples and clones a large number of individual sequences and returns a read for each. During assembly by the company the majority sequence will be selected, but the minority base can be discovered from the original reads so long as there is some indication of where to look, together with an estimate of the frequency.

For example in our data, we noted some possible variability in the following sequence at the underlined locations:

GGTATTTAGGTGTGTGTGAACAGGATTGTAACAAAATAACACCCCAT

Using the tool fr2 with the match_string:

GGTAT.TTAGGTGTG.TGTGAACAGGATTG.TAACAAAA.TAACAC.CCCAT

produces the output shown in Supplementary Table 3 [see Additional file 3]. At the foot of the table are four instances, out of 60 matches, of the alternate form with all five of the polymorphic bases consistently changed. Whilst the error rate in individual reads from the 454 sequencing is quite high, to get five different bases all changed consistently strongly indicates the presence of a minority sequence with a frequency of about 6%.

Assembling highly repeated regions

A number of tools are available to help with the assembly of repeat regions [18-22] but these are not specifically designed for the short reads from pyrophosphate sequencing. Pyrophosphate sequence data presents a particular challenge to assembling repeats because of the short read lengths and higher error rates. However, it also provides two advantages over Sanger sequencing:

1. The fragments generated by the shearing are incorporated at random onto the beads, avoiding the difficulties produced by sequences being under-represented because they are difficult to clone.

2. The large number of reads means that each sequence is read many times.

Repeats longer than about 100 bases become increasingly difficult to resolve, and in particular, exact repeats can only be assembled if there are reads longer than the repeat. However if the repeats are not exact the 454 data will provide information that can lead to at least partial assembly, and will indicate the number of copies even for exact repeats.

The commonest problem with assembling repeat regions is the collapsing of repeats. If the repeats are 100% identical this is difficult to avoid, but when the repeats differ by just one or two bases there is a risk that the assembler will treat the changed base as an error and misalign the sequences. As the tool build accepts only exact matches it can be used to force an assembly through a region of high repeat, so long as the individual base changes occur at least once in 60 bases. Where stretches of sequence occur that are identical for more than 60 bases build will stop at the split point, indicating the two or more possible variants. Each possibility can then be followed independently until they converge once more.

The power of the 454 data to identify misreads and to indicate the number of copies of a given repeat sequence can be illustrated from the GpSGHV data. The virus sequence contains a region of repeats between nucleotide position 178 000 and position 183 000. The repeats consist of alternating 27 base elements designated