PineSAP--sequence alignment and SNP identification pipeline.

UNLABELLED: The Pine Alignment and SNP Identification Pipeline (PineSAP) provides a high-throughput solution to single nucleotide polymorphism (SNP) prediction using multiple sequence alignments from re-sequencing data. This pipeline integrates a hybrid of customized scripting, existing utilities and machine learning in order to increase the speed and accuracy of SNP calls. The implementation of this pipeline results in significantly improved multiple sequence alignments and SNP identifications when compared with existing solutions. The use of machine learning in the SNP identifications extends the pipeline's application to any eukaryotic species where full genome sequence information is unavailable. AVAILABILITY: All code used for this pipeline is freely available at the Dendrome project website (http://dendrome.ucdavis.edu/adept2/resequencing.html)

1 INTRODUCTION

Single nucleotide polymorphism (SNP) detection involves looking across multiple sequence alignments and identifying base discrepancies. The higher the sequence coverage and quality score at a given site, the more confident the SNP prediction. In high-throughput studies, we rely on existing tools to make confident identifications as visual confirmation is not an option. Poor initial alignments can greatly increase both the false positive and false negative rate of SNP predictions. The automated alignment and SNP detection from re-sequencing data where a reference genome sequence is not available is an on-going challenge. Pine (Pinus) sequence data presents an even greater obstacle as it is highly polymorphic (average of one SNP per 50 bases) (Neale, 2007).

Existing programs such as Phrap (Lee and Vega, 2004) are heavily used for aligning genomic re-sequencing data and providing direct input to SNP prediction programs. Phrap, by definition, is a DNA sequence assembler and does not perform well when paired with the task of aligning highly polymorphic re-sequencing samples. It will often place individuals with different haplotypes into separate contigs. If the stringency of Phrap is reduced (in an attempt to force creation a single contig), misalignments of indels leads to poor overall sequence alignments. Quality benchmarks evaluated across several DNA and RNA aligners found ProbconsRNA (Do et al., 2005) to be highly accurate (Carroll et al., 2007; Wilm et al., 2006), however, it proves to be inhibitive in terms of speed for high-throughput studies.

SNP identification solutions that can accommodate fluorescence-based re-sequencing reads, often require genomic reference sequence. The utilities PolyPhred (Nickerson et al., 1997) and Polybayes (Marth et al., 1999) rely primarily on quality scores and sequence coverage. In tests utilizing Polybayes and Polyphred in loblolly pine (Pinus taeda L.), we received at best, 78% prediction accuracy with the majority of the discrepancy resulting from false positives. Recently, machine learning techniques have been applied to the problem of computational SNP discovery (Matukumalli et al., 2006; Unneberg et al., 2005; Zhang et al., 2005). Two of these applications rely on either an existing reference genome sequence (Zhang et al., 2005) or process EST sequences directly rather than raw tracefiles (Unneberg et al., 2005). The SNP-Phage application (Matukumalli et al., 2006) can call SNPs from fluorescent reads without a reference sequence, however, attention to the alignments of highly polymorphic organisms prior to SNP calling is not available.

PineSAP was developed as a high-throughput solution to analyze re-sequencing data in the form of chromatogram files for forward and reverse reads of multiple individuals. The pipeline presented here runs on the Unix/Linux platform and was written in Perl. PineSAP implements a combination of Phred, Phrap and ProbconsRNA to efficiently and accurately call bases and align re-sequencing reads. Following alignment, SNPs and indels are identified through the Polyphred and Polybayes packages. Sequence-based information is extracted and processed through a supervised machine learning algorithm for the purpose of accepting or rejecting the SNP predictions.

2 ALIGNMENT

The alignment section of the pipeline implements a de-coupled and modified version of phredPhrap on the original chromatogram files. Phred (Ewing et al., 1998) is responsible for the base calls and the assignment of quality scores. Phred is integrated into this pipeline with parameters to trim low quality ends in order to prevent alignment issues from bases in these regions. Conservative Phrap parameters are applied to prevent any misalignments in the resulting contigs and to ensure that all reads are retained.

Contig consensus sequences are exported from the ace format file generated from phredPhrap and aligned with ProbconsRNA. For each contig, an aligned FASTA file is created with each read in the contig aligned to the consensus sequence. These files along with the aligned file of all contig consensus sequences are used to generate a single multi-sequence FASTA file. In this file, each read is aligned to an overall consensus for the amplicon based on the alignment in the ProbconsRNA output and each read's alignment to the consensus sequence of its member contig. The final multi-sequence aligned FASTA file is converted back to an ace formatted file suitable for input to Polybayes and Polyphred.

3 SNP CALLING

The purpose of the classifier is to evaluate the accuracy of the SNP calls resulting from Polyphred and Polybayes. Sequence-based statistics were derived through a customized feature extraction program and fed as a vector for each polymorphism to the J48 classification tree available in the WEKA classifier package. The final set of features identified fully represents the local and global sequence variation, alignment depth and quality, local and global base quality and sequence alignment quality. From the set of 300 used for training, all of the true positive and false negative vectors are represented as real SNPs and the true negative and false positive vectors as non-SNPs. The resulting SNP calls were extracted from their respective output files flanking sequence, quality scores and a normalized confidence score.

4 VALIDATION

Manual validation of the alignments was completed with the same 300 loblolly pine amplicons. Each amplicon had between 8 and 36 reads. With Phrap alone and default parameters, 23% of the amplicons were placed in a single contig, 27% into two and 20% into three. Amplicons with two or more contigs could be forced into a single contig 82% of the time, however there are problems with the alignment in 34% of cases. The alignment method implemented in PineSAP improves the success rate to 98%.

When evaluated in terms of speed, the alignment method described above was run with a straight ProbconsRNA implementation. We determined it would take ∼25 times longer to process 36 sequences per amplicon.

The classification tree generated from the training sequences was tested against a unique set of 120 sequences with 563 manually validated SNPs. All SNP calls were identified as based on visual inspection of Polyphred and Polybayes predictions in Consed (Gordon et al., 1998). The classification tree resulted in a significant overall improvement with a calculated accuracy of 93.6% (Table 1).

Table 1.

Results of SNP prediction on the test sequence data

Evaluation	J48	Polyphred	Polybayes
Accuracy	93.6	76.25	78.02
Sensitivity	88.21	83.22	86.54
Specificity	98.73	N/A	N/A