phrap ("phragment assembly program", or "phil's revised assembly program") -- a program for assembling shotgun DNA sequence data. Key features: allows use of entire read (not just highest quality part); uses a combination of user-supplied and internally computed data quality information to improve accuracy of assembly in the presence of repeats; constructs contig sequence as a mosaic of the highest quality parts of reads (rather than a consensus). N.B. phrap does not provide editing or viewing capabilities; these are available with consed or phrapview. cross_match -- a general-purpose utility (based on SWAT, an efficient implementation of the Smith-Waterman algorithm) for comparing any two sets of (long or short) DNA sequences. For example, it can be used to compare a set of reads to a set of vector sequences and produce vector-masked versions of the reads; a set of cDNA sequences to a set of cosmids; contig sequences found by two alternative assembly procedures (e.g. phrap and xbap) to each other; or phrap contigs to the final edited cosmid sequence. It is slower but more sensitive (because it allows gaps) than BLASTN. I) INPUT FILES are of two types: sequence files (required), and quality files (optional). Ia) Sequence files. With both phrap and cross_match, either a single sequence file or multiple sequence files can be specified. If a single file is specified, all sequences in it are compared to each other; if more than one file is specified, sequences in the first file are compared to sequences in the second (and subsequent) files. Typically phrap is used with a single input sequence file containing the reads, unless one wants to assemble one set of reads "against" another sequence or set of sequences (see below). Phrap is designed to be able to use all of each read sequence in the assembly, not just the trimmed (highest quality) part, so the full sequence of each read should be provided. Cross_match is generally used with two input sequence files, the first containing the "query" sequences and the second containing the "subject" sequences. Sequence files must be in FASTA format: each sequence should have a single header line with leading character '>' followed immediately (i.e. with no intervening spaces) by the read name, and then (optionally), on the same line, a space or spaces followed by descriptive information about the read; the header line is followed by a separate line or lines containing the sequence. The read names, descriptive information, and sequence can be of arbitrary length (up to a combined total size of about 1 Gb), since space for them is allocated dynamically. Read names: It is assumed by phrap that the part of the read name up to the first '.', if any, is the name of the subclone from which the read is derived; this is used, for example, in checking for chimeric reads (for which purpose it is important to know when two reads come from the same subclone, since another read from the same subclone should not be regarded as independent confirmation of a join). IT IS IMPORTANT THAT DIFFERENT READS FROM THE SAME CLONE BE GIVEN NAMES THAT ADHERE TO THIS CONVENTION -- OTHERWISE THEY WILL BE INTERPRETED (INCORRECTLY) TO CONFIRM EACH OTHER, WHICH MAY ADVERSELY AFFECT ASSEMBLY. Orientation of the read (in cases where data is acquired from each end of a subclone insert) is assumed to be indicated by the first letter following the '.', namely 's' or 'x' for forward reads, and 'r' or 'y' for reverse reads. At present the read orientation is used only for consistency checking following assembly, and not in the assembly itself. In addition, it is assumed that dye terminator reads are indicated by an 'x' (for forward reads) or 'y' (for reverse reads) following the '.', and that 's' and 'r' refer to dye primer reads. This is used in the quality adjustment procedure (see below), where confirmation of a non-dye-terminator read by a terminator read counts as significantly as confirmation by an opposite-strand read. If the read name ends in ".seq" the .seq is removed. Descriptive information on the header line is not used by the program at present, although it is read in, and is written back out to the .ace and .singlets files. Sequence: Non-alphabetic characters (including '*') in the sequence lines are automatically stripped out when the file is read in, except for '-' which is converted to 'N'; '>' must not appear anywhere within the sequence since it is assumed to start the header of a new sequence (even if it is not at the beginning of a line). Lower case letters are converted to upper case on readin, so it is no longer possible (as in previous phrap versions) to use case as a crude indication of quality. Vector screening: Following creation of a FASTA file with the raw read sequences, it is important to remove or screen out sequencing (subclone) vector sequence. (It is unnecessary to remove the cloning (e.g. cosmid) vector, unless the inserts from multiple overlapping cosmids are being assembled simultaneously; however it is generally useful to remove at least the central part of the cosmid vector, since this then provides a natural point at which phrap can break the circular cosmid sequence into a linear one). Unremoved sequencing vector may cause reads to be identified as "possible chimeras", or otherwise interfere with proper assembly. Some vector removal programs may be unreliable at identifying vector sequences that are rearranged or found in the lower quality part of the read, and I recommend that you screen for sequencing vector using cross_match, whether or not you have already screened them using another program. First, create a FASTA file, say "vector", with all of the vector sequences you want to screen for (I use one that contains all vector sequences used in any of the ongoing sequencing projects, in case the vector is misidentified in the current project). If your read file is named "C05D11.reads" (for example), then run the command cross_match C05D11.reads vector -minmatch 12 -minscore 20 -screen > screen.out (This uses somewhat more sensitive parameter settings than the cross_match defaults). The '-screen' option causes a file named "C05D11.reads.screen" to be created, containing "vector masked" versions of the original sequences: i.e. any region that matches any part of a vector sequence is replaced by X's. This ".screen" file should be used in all subsequent analyses. The output from cross_match (which has been redirected to the file screen.out in the above example) lists the matches that were found (see below). N.B. If you have created a .qual file for your reads (say C05D11.reads.qual), be sure to rename or copy it to C05D11.reads.screen.qual so that it will be recognized by phrap when C05D11.reads.screen is used as the input FASTA file. Ib) Input of quality information. For each input sequence file in FASTA format, you may optionally include a corresponding file of data quality information. This is strongly recommended since it greatly improves the accuracy of assembly and of the consensus sequence. The name of this file should consist of the name of the FASTA file, with ".qual" appended. The format of the .qual file is similar to that of the corresponding FASTA file. For each read there should appear a header line identical (except possibly for the "description" field) to that in the FASTA file. This is followed by one or more lines giving the qualities for each base. Quality values should be integers between 0 and 99 (inclusive), and should be separated by spaces. No other (non-digit, non-white space) characters should appear. The total number of quality values for the read must match the number of bases that appear in the FASTA file. Also, the orders of reads in the FASTA and corresponding .qual file must be identical. Bases 'N' or 'X' are automatically assigned quality 0, whether or not a .qual file is provided. If provided, quality information is used by phrap both in the assembly itself (where discrepancies between high quality bases are used to discriminate repeats from true overlaps) and in determining the contig consensus sequence (which is formed as a mosaic of the highest quality parts of the reads). Because phrap generates its own "quality" measures (based on read-read confirmation) it performs reasonably well even if no input quality is provided, but when available the input data quality information can substantially increase the accuracy of the consensus, particularly in regions where there are strand-specific effects on quality (e.g. compressions) since in such cases the input quality information may constitute the only basis for choosing one strand over the other. If your sequencing data is from ABI 373 or 377 automated sequencers I strongly recommend that you generate the read sequences using the basecalling program phred, which will produce an appropriate quality file for input to phrap. Phred generally gives more accurate base calls overall, for a longer distance, than the ABI software does. On the basis of the trace characteristics, phred computes a probability p of an error in the base call at each position, and converts this to a quality value q using the transformation q = -10 log_10(p). Thus a quality of 30 corresponds to an error probability of 1 / 1000, a quality of of 40 to an error probability of 1 / 10000, etc. The quality value of 99 should be reserved for base calls that have been visually inspected and verified as "highly accurate" (during editing), and 98 for bases that have been edited but are not highly accurate (these are converted to quality 0 in phrap). If two reads appear to match but have discrepancies involving bases that have quality 99 in both reads, phrap does not allow them to be merged during assembly. At present this is the only way to break false joins made by phrap. If you wish to generate your own quality values you should be prepared to experiment. It is particularly important that possible insertion errors (which are fairly frequent late in the trace with ABI basecalls) receive quality 0, because in choosing the "consensus" phrap tends to favor the reads with the highest total quality in a given region. Also, if at all possible you should use quality values that are related to error probabilities in the manner described above, since to some extent phrap is tuned to expect these (and its output qualities will then have a similar interpretation). If all input quality values are relatively small (less than 15), phrap assumes that they do not correspond to error probabilities and attempts to rescale them so that the largest quality value is approximately 30. This allows different input scales to be used. Phrap adjusted quality values/error probabilities: Phrap computes adjusted quality values for each read on the basis of read-read confirmation information, as follows. If a read is confirmed by an opposite-strand or different-chemistry (dye terminator vs. dye primer) read at a given position, that position is given a quality which is the sum of the two input read qualities; when more than one opposite-strand read confirms a given position, only the single highest quality from all opposite-strand matching reads is used. If qualities are related to error probabilities as described above, this procedure can be interpreted as computing an error probability for each base which is the product of the error probabilities for the two reads; since error profiles for opposite-strand or different chemistry reads are essentially independent, this is reasonable, although it is somewhat conservative in that it does not take into account same-strand matches (which clearly should count for something, although they certainly are not independent). Contig positions are assigned quality values equal to the highest adjusted quality of any read at that position, and then adjusted downward to take into account any discrepancies with other reads. As a result, when phred input qualities are used, the output phrap qualities associated to each contig position have a natural interpretation as (conservative) error probabilities. Such error probabilities provide an extremely useful guide to where editing or additional data collection is needed. The phrap adjusted qualities are used in computing "LLR scores" for each pairwise match between two reads. These scores take into account the qualities of the base calls in the reads, and are (approximate) log-likelihood ratios for comparing the hypothesis that the reads truly overlap to the hypothesis that they are from 95% similar repeats. The point is that discrepancies between overlapping reads are due to base-calling errors and thus tend to occur in low-quality bases, whereas reads from different repeats can have high-quality discrepancies that are due to sequence differences between the repeats; and the probability of the observed data under each hypothesis can be quantified using the interpretation of the phrap qualities in terms of error probabilities. A pairwise match tends to have a positive LLR score if the two reads overlap, whereas it tends to have a negative LLR score if the two reads are from different repeats (unless the repeats are nearly identical). II) RUNNING phrap. The command line should be of the form phrap seq_file_name1 [seq_file_name2 ...] [-option value] [-option value] ... The output is extensive and should be redirected to a file. If more than one seq_file_name is provided, then reads in the first file are assembled only against sequences in the 2d (and later) files; no sequences are compared to other sequences in the same file. In this case, many features of phrap (e.g. checking for anomalies, vector, etc.) which rely on comprehensive read-read comparisons are turned off, SWAT scores are used in place of LLR_scores, and implied merges are not performed -- i.e. each read in the first file is merged with at most one sequence in the 2d file. This usage of phrap is convenient for aligning a set of reads against a known "finished" sequence in order to look for polymorphisms, for example. Example: phrap C05D11.reads.screen -minmatch 12 -ace > phrap.out The following options are available for both cross_match and phrap: option name value default value -penalty mismatch (substitution) penalty for SWAT comparisons -2 -gap_init gap initiating penalty for SWAT comparisons penalty - 2 -gap_ext gap extension penalty for SWAT comparisons penalty - 1 -minmatch minimum length of matching word to nucleate SWAT comparison; 14 if minmatch = 0, a full (non-banded) comparison is done -minscore minimum SWAT score 30 -bandwidth half band width for banded SWAT searches 14 (full width is 2 n + 1) -matrix score matrix for SWAT comparisons (supersedes -penalty) [default has +1 for matches, -penalty for mismatches] -raw use raw rather than complexity-adjusted Smith-Waterman scores. (default is to adjust scores to reflect complexity of matching sequence). (This is a flag; no value should be provided). -tags tags selected output lines, to facilitate parsing (This is a flag; no value should be provided). -indexwordsize size of indexing (hashing) words in words.c. The value of 10 this parameter has a (small) effect on run time and memory usage. The following options apply only to phrap: -maxgap maximum gap (unmatched region) allowed in merging contigs 30 -trim_start no. of bases to be removed at beginning of each read 0 -vector_bound no. of bases at beginning of each read, matches within which are assumed to be vector 60 -trim_penalty penalty used for identifying degenerate sequence at beginning/end of read -2 -trim_score minimum score for identifying degenerate sequence at beginning/end of read 20 -confirm_length minimum size of confirming segment (starts at 3d distinct 8 nuc following discrepancy) -confirm_trim amount by which confirming segments are trimmed at edges 1 -confirm_penalty penalty used in aligning against "confirming" reads -5 -confirm_score minimum alignment score for a read to be allowed to "confirm" part of another read 30 -qual_show LLR cutoff for flagging "low_quality" regions in 20 contig sequence and "high quality" discrepancies between read and contig. Bases in the .contigs, and .ace file are lower case iff their LLR-converted quality is below this value. -node_seg minimum segment size (for purposes of traversing weighted directed graph) 8 -node_space spacing between nodes (in weighted directed graph) 4 -revise_greedy splits initial greedy assembly into pieces at "weak joins", and then tries to reattach them to give higher overall score. (This is a flag; no value should be provided). -ace create ".ace" file with the assembly information padded reads, etc. (This is a flag; no value should be provided). -view create ".view" file suitable for input to phrapview (This is a flag; no value should be provided). The following options apply only to cross_match: -masklevel controls which matches (with scores at least minscore) are 80 displayed. A match is not displayed if the percentage of aligned query bases in the match that are contained within the aligned bases of higher-scoring matches is masklevel or higher. For example: if masklevel = 0, only the single highest scoring match for each query is displayed 100 any match that is not completely contained within a higher scoring match is displayed 101 all matches are displayed. -alignments display the alignment for each match (Warning -- may produce huge amount of output!) (This is a flag; no value should be provided). -discrep_lists give list of discrepancies for each match (This is a flag; no value should be provided). -discrep_tables give table of discrepancies (by quality) for each match (default is to display a single table, that combines data for all matches) (This is a flag; no value should be provided). -screen create a ".screen" file in which regions in the first input sequence file that match some region in the second (or later) file are replaced by X's. (This is a flag; no value should be provided). Normally you should just use the default values of the parameters; I haven't systematically explored optimal values for these, but the current ones appear to work reasonably well. To increase (potential) sensitivity of the pairwise comparisons, at the cost of a longer run time, try decreasing minmatch, minscore, and/or penalty, and increasing bandwidth. Phrap output files (with the exception of the standard output, all output files are given names consisting of the input read file, with an appropriate suffix (given below) appended): 1) .singlets file (named C05D11.reads.screen.singlets in the above example): A FASTA file containing the singlet sequences (i.e. the singleton reads with no match to any other read). 2) .contigs file (named C05D11.reads.screen.contigs in the above example): A FASTA file containing the sequences (without pads) of all contigs (upper case <=> quality is >= qual_show). These include singleton contigs consisting of single reads with a match to some other contig, but that couldn't be merged consistently with it. Thus reads that would go into xbap's "failed file" instead are considered separate contigs by phrap and appear with the other contigs. 3) .contigs.qual file: The corresponding quality file (see above for description of how phrap computes qualities for contig positions). 4) .log file: A summary of how the assembly proceeded, and how the contig sequences were constructed as mosaics of the individual reads. This is probably only of interest to me, for trouble-shooting purposes. 5) .ace file: (Produced only when the -ace option is used). A ".ace" formatted file containing (i) the sequences of the padded contigs and reads; (ii) the locations of the reads with respect to the contig sequence; (iii) information about how the contig sequence is pieced together as a mosaic of the high quality parts of the reads; (iv) an indication of where the reads can be "clipped" to remove low quality sequence at either end; (v) qualities of the bases in the contig sequence. Additional information will be added in the future. Sequence information is indicated by a header line consisting of the word DNA followed by the name of the sequence, followed by lines giving the sequence. Sequences are padded (with a '*' character) so that reads align directly against contigs, i.e. the alignments can be reconstructed exactly using only the information in this file. The end of the sequence is followed by a blank line. For reads in reverse orientation with respect to the contig, the sequence of the complement is given instead of the original sequence, and the suffix ".comp" is appended to the name of the read. The case (upper or lower) of bases may differ from that in C05D11.reads.screen, reflecting adjustment of quality values by phrap (upper case <=> LLR_quality is at least qual_show). Examples: DNA x72g12.s1.comp ccaatgtgacttggataatggcgtaaggcttgactttgacgtcccagtt* aa*ccctg*ctt*aaact*accttc*c*tta*tcctaaac*attaac*tt cttagctcaa*t*cccattcca*gcgaatcttca*aactgtcaaactana *tcgt*gaaatc*tcctccTGCTAGATtc*tcaaATCTGTCAATTGAAGA TCGTAGTTACTTGTTCCTCTACACAGATTTGCTTTTTGAATCTCCTGCCA TGATCGATGGAGTCTTGACATCTGCTGATGATGTTGCCAAACATTTCAC* GAAAGAT*CTCATTGATCACTCTATTCAAGT*GGGAGTC*TCCGGTCTTT ATGATCGATTTGTGAATCTTCGTATCAAAGTTGGAGCTGACAAGTATCCA TTGCTTGCGAAATGGGCTCAAATTTTCACTCAGGGAGTCGTCTTCGATCC TTCAAGAATTCATCAATGTGCTCAAAAGTTGGCTGGAGAAGCTCGTGATC GGAAGAGAGATGGATGTACTGTGGCATCAACTGCAGTAGCTTCAATGGTT TATGGAAAGAGTATGTTATTt DNA Contig35 GATAAGataATAAtGGAAAATaGAAaccGGAaAaATaATAAaATaaTTTc aGATcGcTGaAGAaGAaGaGAAGaGAATAGcAGccCaATGTGAGAAGCTC GGcAAAAAAGGACTCGAAGaaGcGGGAAAGAGTCTgGAAGcTGCCATTCT TGAGAACACTGcCAATcATcCaaGTGCTGAACTTTTGGATCAATTGATCG TCAAAGATcTTGAAGCTTTTGATCGTTTCCCAGTTCAATCCCTGACTTCA AACTCAcCTTCTC*TTACTCCTCAACAATCAACTTTCTTAGCTCAATT*C CCATTCCACGCGAATCTTCACAACTGTCCAACTAAATTCGtGGAAATCTT Clipping information for a read is indicated by a line consisting of the word "Sequence" followed by the name of the read, followed by two lines that each indicate the beginning and ending bases that remain after low quality (quality 0) sequence is removed from each end of the sequence. The line starting with the word "Clipping" gives the base coordinates with respect to the unpadded sequence, while the one starting with "Clipping*" refers to the padded sequence. If (and only if) the input fasta file included a "description" field (on the header line), that appears as a separate "Description" line following the clipping information. Example: Sequence x72g12.s1.comp Clipping 153 548 Clipping* 170 570 Description ALF run #5 The assembly information is indicated by a line consisting of the word "Sequence" followed by the name of the contig sequence, followed by lines giving the positions of the reads in this contig. For each read there are two lines: one starting with the word "Assembled_from" and giving the name of the read together with its starting and ending positions with respect to the unpadded contig sequence; and one starting with the word "Assembled_from*" which instead has positions with respect to the padded contig sequence. (In both cases the starting and ending positions are those of the first and last base of the read, not those of the "clipped" part). Note that it is possible for a read to overhang the end of the contig sequence, so that start positions may be negative and end positions may be greater than the number of the last base in the contig. Information about how the contig sequence is pieced together is given in the form of "Base_segment" lines: each such line gives the starting and ending positions of the segment within the contig sequence, the name of the read from which the sequence is taken, and the starting and ending positions of the corresponding segment within that read. Example: Sequence Contig35 Assembled_from z14f12.s1 1 596 Assembled_from* z14f12.s1 1 602 Assembled_from z13f5.s1 58 566 Assembled_from* z13f5.s1 58 572 Assembled_from z25a11.s1 115 738 Assembled_from* z25a11.s1 115 744 Assembled_from z14d4.s1.comp 145 741 Assembled_from* z14d4.s1.comp 145 747 Base_segment* 350 412 z26g6.s1 70 132 Base_segment 343 347 z13h11.s1 119 123 Base_segment* 345 349 z13h11.s1 121 125 Base_segment 339 342 z26g6.s1 60 63 Base_segment* 341 344 z26g6.s1 61 64 Base_segment 332 338 z13f5.s1 275 281 Base_segment* 334 340 z13f5.s1 277 283 Base_segment 331 331 z26g6.s1 52 52 Base_segment* 333 333 z26g6.s1 53 53 Base_segment 222 330 z25a11.s1 108 216 Base_segment* 222 332 z25a11.s1 108 218 Base_segment 166 221 z13f5.s1 109 164 Base_segment* 166 221 z13f5.s1 109 164 Base_segment 142 165 z25a11.s1 28 51 Base_segment* 142 165 z25a11.s1 28 51 Base_segment 77 141 z14f12.s1 77 141 Base_segment* 77 141 z14f12.s1 77 141 Base_segment 61 76 z13f5.s1 4 19 Base_segment* 61 76 z13f5.s1 4 19 Base_segment 1 60 z14f12.s1 1 60 Base_segment* 1 60 z14f12.s1 1 60 Quality information is indicated by a line with the word BaseQuality followed by the name of the contig or read; this is followed by a line or lines giving the numerical quality values. An arbitrary number of values may appear on each line. As with other data, the end of the BaseQuality information for a given sequence is indicated by a blank line. Example: BaseQuality Contig3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 31 21 24 23 23 22 22 22 22 30 31 20 20 20 20 20 20 30 30 31 30 33 33 34 34 34 34 31 33 31 34 34 30 30 25 25 25 26 30 30 33 37 30 33 36 35 23 23 24 24 24 24 30 30 30 30 30 30 30 30 30 33 33 33 33 30 30 33 33 33 30 33 33 29 29 29 29 33 34 33 33 21 21 21 21 23 23 23 23 23 23 24 24 26 26 26 26 26 26 26 40 33 30 30 30 30 33 33 33 33 33 30 33 33 33 33 33 33 33 31 31 31 31 34 37 40 40 45 37 52 52 52 52 52 52 55 55 59 64 65 68 58 57 49 49 49 49 49 49 49 56 60 53 60 53 49 49 49 49 49 34 45 45 45 45 30 25 25 25 28 30 45 45 49 49 49 49 49 70 60 60 59 59 53 56 53 53 53 55 53 53 60 45 49 49 42 42 42 45 45 36 33 34 49 46 46 46 54 54 59 53 49 49 42 41 43 40 44 40 49 49 49 54 58 53 52 57 51 51 41 39 15 15 35 35 [N.B. MUCH OF THE FOLLOWING HAS NOT BEEN RECENTLY UPDATED AND MAY BE OBSOLETE]. 6) standard output (redirected to phrap.out in the above example): contains a list of the contigs and the reads they contain, with information about how good the alignment of the read against the contig sequence is. Also identifies matching regions within and between contigs, suspect reads (probable deletions or chimeras), indicates possible problems with assembly or contig sequence (e.g. discrepancies between reads and contig sequence), and gives results of the forward/reverse read consistency checks, There is quite a lot of information in this output. It has the ultimate goal of giving a complete picture of the final assembly, including possible sites of assembly error; some of this information would probably be better represented in graphical form. The output begins by giving summary information about the parameter settings and input data file(s) (sequence and quality). It then lists (i) regions at beginnings or ends of reads which consist largely (>50 %) of a single nucleotide, usually due to degenerate signal. These are internally converted to N's, to prevent them from causing spurious matches which could interfere with assembly. (Note that the .pads file output by phrap will contain the N's rather than the original sequence). (ii) Exact and near-exact duplicate reads. Exact duplicates, which probably represent duplicate entry of the same data, are excluded from assembly. The near-exact duplicates (which are only listed if they appear to come from different clones, according to the nomenclature convention assumed by phrap) should be examined to make sure they really are from different clones. If not, you are not using the correct read nomenclature. (iii) probable unremoved sequencing vector, detected as regions at beginnings of reads which match only the beginnings of other reads. This method is not foolproof and should not be considered a substitute for screening out sequencing vector using cross_match prior to running phrap. The unremoved vector may simply be inaccurate vector sequence, that did not match the correct sequence in the original cross_match screen but does match other reads in which similar errors occur. In any case, since the procedure in phrap is not foolproof, you should be sure to add the (potentially erroneous) sequence that is displayed in the phrap output to your vector file, and redo the cross_match vector screen, to increase the likelihood that you catch all of the vector. If you find examples where sequence that is not vector is erroneously being labelled as such by phrap, please let me know. (iv) reads which have off-diagonal matches to themselves (due e.g. to short tandem repeats), in the same orientation. This information is ignored at present, but ultimately will be used to improve assembly of tandem repeat regions. Due to the use of complexity-corrected Smith-Waterman scores some commonly occurring microsatellites may not be detected. (v) probable chimeric reads, identified as reads whose confirmed parts break into two pieces, separated by at most 30 bp and each having a prematurely terminating alignment with some other read. Chimeric reads are excluded from the assembly and instead appear as singlet contigs. The condition that prematurely terminating alignments exist is necessary to avoid falsely identifying "single links" as chimeras. Since repeats also may cause prematurely terminating alignments, it is possible (but unlikely) for a single link non-chimeric read containing two or more repeats to be falsely identified as a chimera. Chimeras in which the first segment consists of thirty or so bases near the beginning of the read are probably due to unmasked sequencing vector. In fact, it is useful to test for unremoved vector in as sensitive a fashion as possible by running phrap with -minscore and -confirm_score both set to 25, and then examining the chimera list (once sequencing vector has been removed the assembly itself should be done with the default parameter values). (vi) probable deletion reads, identified as reads which match some other read but have a gap with respect to it of >10 bp (more specifically, read 1 is determined to have a deletion with respect to read 2 if there are two pairs of matching segments between the reads such that the segments in read 1 are nearly adjacent (within 20 bases) and the separation between the segments in read 2 is at least 10bp larger than the separation in read 1), and for which there is no other clone that shows the same deletion (the latter condition is needed to rule out the possibility of repeated sequences for which some copies have deletions). For each deletion clone and each such matching read, the apparent size of the deletion, location of the deletion and (in parentheses) name of the matching (undeleted) read and extent of the segment that was deleted from in it are given. As with the probable chimeras, the probable deletions are not used in assembly and instead appear as singlet contigs. (vii) Summary of "rejected" pairwise alignments; i.e. matches which either prematurely terminate (do not extend as far as they should, given the apparent sequence quality), or include mismatches between high quality bases. In such cases the reads are probably from different repeats, and so the alignments are not used in the assembly (N.B. at present incompletely extended matches are considered during the assembly, since a more sensitive criterion used there will reject them if warranted). (viii) (N.B. This part of output needs revision). Summary of information about accepted pairwise alignments, including no. of confirmed reads (i.e. reads matching some other read), their average lengths (total, confirmed, "strongly confirmed" (i.e. matching a reverse sense read) and trimmed); a crude estimate of the cosmid size (assuming all non-rejected matches are real -- if there are near-perfect repeats, the size may be underestimated) and depth of coverage; a histogram of depths (letting the depth of a read be the maximum depth that occurs in it); a summary of the discrepancies between matching reads (by strand sense, nucleotide and quality); and a histogram of spacings between adjacent indel pairs. By sense of the aligned pair (reads in same direction, reads in reverse directions). (ix) isolated singletons (reads having no match with score at least minscore to any other read). (x) contigs (including singlet contigs for which the read did have a match with some other read, but could not be consistently assembled with it). The contigs are sorted by (increasing) number of reads, so that the singlet contigs appear first. The following information appears for each contig: the number of reads; the total length of the contig sequence (not including pads); the length of the trimmed sequence (i.e. after regions consisting entirely of lower case letters, N, and X have been removed from either end of the contig); and a list of the reads in the contig and information about their alignments. For each read the following information is given: a 'C' if the read is in reverse orientation; the starting and ending positions for the read (i.e. its full length, not just the part that aligns with a given penalty) with respect to the padded sequence; the read name; the score of the alignment of the read against the contig sequence; (in parentheses) the highest score of a match of the read against some other read in a different contig or elsewhere in the same contig (so reads for which this number is non-zero are those which overlap a repeat, or an incorrectly or incompletely assembled region); the per cent mismatch, insertion, and deletion rates for the alignment of the read against the contig sequence; the number of bases at the beginning of the read which were not included in the alignment; (in parentheses) the number of bases at the beginning of the read which did not align against any other read (i.e. were not confirmed); and similarly for the bases at the end of the read. If the unparenthesized number is substantially greater than the parenthesized number that accompanies it, there could be a problem with the alignment. Note however that different penalties are used for aligning a read against another read as opposed to aligning it against the contig sequence. Following the read list is information to be used in assessing the quality of the assembly. This includes: A description of regions in the contig sequence whose quality has been adjusted on the basis of alignments of reads to the contig (e.g. regions that had double-stranded confirmation in the initial pairwise comparison may not have when the contig is assembled, due to resolution of repeats; and regions that were not confirmed in the initial pairwise comparisons may be confirmed in the assembled contig, due to the use of a lower gap penalty which gives more complete alignments). Histogram of the qualities of the contig bases; the extents of the leading and trailing quality 0 regions of the contig (such regions are likely to be highly inaccurate); and a list of the bases that have quality < qual_show. Slack and HQ Mismatch histograms (not yet documented here). A table showing "gaps" on each strand (i.e. regions for which there is no aligned part of any read, so that additional data, or possibly editing to permit alignment of existing data, is required); the closest read that could potentially be extended or edited to cover the gap; whether or not the inaccurate, unaligned part of that read already covers the gap (in which case editing might be sufficient to close it); and the length of an accurate extended read (from the same start in the same clone) that would be required to cover the gap. The table appears following the list of reads in the contig and their positions. An example (from a contig in C05D11): Gap Size Closest read (Start) Covers Read length required now? to cover Top strand: left - 244 244+ 1809 - 1880 72 x72c5.s1 (1382) No 499 2967 - 2993 27 ae82b07.s1 (2562) Yes 432 5392 - 5536 145 z17c5.s1 (4965) No 572 6036 - 6320 285 z13f6.s1 (5537) No 784 8137 - right 0+ z26c9.s1 (7667) No 469+ Bottom strand: left - 0 0+ z17b4.s1 ( 563) No 563+ 1138 - 1139 2 z14f3.s1 (1143) Yes 6 6940 - 7230 291 z17c10.s1 (7666) No 727 7717 - right 420+ The gaps that say "left" or "right" are the ones at the left and right ends of the contig (the gap size in this case is the part of the existing contig sequence at that end that is not double stranded, with the plus indicating that there is an additional gap to the left or the right of the contig; these gaps are always designated not covered ("No")), and the relevant read in this case is the one that should be extended to close gaps between contigs, while the other gaps are internal gaps and the relevant reads are the ones required to get double-stranding. Note one problem in this example: the bottom strand gap of size 2 lies in the inaccurate 5' part of the z14f3.s1 read, so that would probably not be an appropriate read to edit or get more data from. Need appropriate cutoff values (i.e. size of 5' part of read to ignore; this obviously depends on whether the 5' end has already been trimmed prior to inputting it to phrap). Also, note that since the full read lengths are used by phrap, the ends of the contigs may include some very inaccurate sequence; so the size of the end gaps (i.e. the total amount of inaccurate sequence there) could be larger than indicated; this isn't particularly important though since in any case one will always want to extend those reads to try to close the gaps between contigs. Phrap uses a lower gap penalty (-2) in aligning reads to the contig sequence (to allow more complete double-stranding than the -9 default used for aligning reads to each other). Gaps with "No" in the Covers? column will always require more data (an extended read or walk step, depending on the required length). Those with "Yes" may be coverable by editing of the appropriate read; but it may be simpler to automatically get longer reads for them also, particularly if the gap size is large, since substantial editing would probably be required in that case to use the existing data. Walking primers could be output directly from OSP analysis of the contig sequence (high quality part, i.e. all upper case letters) in cases where the gap looks too large to close with an extended read (presumably "finish" already does this?) Following the gap table, there is a table giving, for each quality value, the total number of bases of that quality in the reads for the contig; the number which are aligned (i.e. included in the SWAT alignments of the reads against the contig); the cumulative no. of aligned bases (for this and higher qualities); the number of bases not included in the alignment (but potentially alignable -- i.e. not lying before the beginning or after the end of the second sequence); the number of discrepancies of each type (substitution, deletion, insertion); the total # of discrepancies (for this quality level) and their percentage (of the aligned bases of the given quality), and the cumulative # of discrepancies and their percentage (of the cumulative no. of aligned bases). The regions at the ends of the sequence that consist entirely of quality 0 are reassigned quality values of -1, since these usually correspond to the lowest quality data at the extreme ends of reads. (Not done in the contig quality histogram.) An 'N' in either a read or the contig sequence is always counted as a substitution error. Depth 0 regions. Following this there is more specific information about possible problems with the assembly: Unaligned segments: parts of reads of length at least 10 that did not align to the contig sequence (the part of the read having quality -1 is not considered). High quality discrepancies: sites in the contig sequence for which at least one read with a quality at least qual_show disagrees with the contig sequence. Most errors in the contig sequence will be found either among these sites or in the regions where the contig quality is < qual_show. Low quality suspects: sites in the contig sequence of quality < qual_show (but greater than 0) for which there is a discrepant read of equal or greater quality. Such sites are particularly error-prone. The list gives the site position, base in the contig sequence and its quality, and quality of the discrepant read. [Formerly also printed out: for each of the two strands, the # of reads on that strand having a substitution, insertion or deletion (with respect to the contig sequence) at that position, and depth (# of reads whose alignments extend across that position). The site position is preceded by a '<' (resp. '=') if there is a read with a discrepant higher (resp. equivalent) quality base at that position.] There then is given a list of the regions that are not covered by "unique-reads" (i.e. reads with no non-rejected matches); and of regions that match (as detected by one or more non-rejected pairwise alignment between reads) other contigs or other regions in the same contig. These are likely to be either true near-perfect repeats or reflections of an error (of omission) in the assembly. To help distinguish these possibilities, such regions are designated spanned (if some read in the current contig extends all the way across it and includes sequence to either side of it; a list of all such reads is given) or "UNSPANNED" (if no such read exists). Misassembly errors of commission generally arise from true near-perfect unspanned repeats, and should result in an "UNSPANNED" flag for at least one of the two matching regions. However not all of them do -- while the fact that repeats exist should generally be detected by phrap, the full extent of the repeat may in some cases not be detected, and the apparent repeat may be "spanned" even though the full (true) repeat is not. Moreover in some cases of misassembly only one of the matching contig regions may be flagged as unspanned. Misassembly errors of omission should always result in "UNSPANNED" matching regions at the end of one or both contigs. Regions of a non-singlet contig matching an anomalous (chimeric or deleted) read are indicated only in the output for the singlet contig that contains the anomalous read. (ix) Following the contig information, there is a summary of the forward/reverse read consistency checks. For each pair of reads with the same clone name are given the contig no. and orientation within the contig (0 = top strand, 1 = bottom strand) for each read; an "L" (for "link") if the reads occur in different contigs; a '*' if there is an inconsistency (e.g. forward and reverse reads which are not pointing towards each other, or two forward reads which are not on the same strand); distance between read starts for forward/reverse pairs (should correspond to clone insert size, and be positive); and distance between starts for forward/forward or reverse/reverse pairs. N.B. IN ALL PLACES WHERE A POSITION IS GIVEN, THE POSITION IS WITH RESPECT TO THE UNPADDED SEQUENCE, NOT THE PADDED SEQUENCE. III) RUNNING cross_match. The command line should be of the form cross_match file_1 [file_2 ...] [-option value] [-option value] ... The sequence files should be in FASTA format. If there are at least two files, all sequences in the first file are compared to those in the second (and any subsequent) files; if there is only one file, all sequences in this file are compared to each other (N.B. at present when there is a single input file, entries are not compared to their own complements, but are compared to themselves in the same orientation and to both orientations of all other reasons). Any matches are printed out. Options are as described above. Output: (N.B. The following description is incomplete -- does not reflect all of the available options). The cross_match output lists matches (whose scores exceeding minscore) between any sequence in the first input file (the "query" sequences) and a sequence in the second (or later) input file (the "subject" sequences). By default (i.e. if allmatches is not set) a match is not shown if its domain is completely contained within the domains of higher scoring matches, where the "domain" of a match is the region in the query sequence that is defined by the alignment start and stop. In the output, matches are ordered by query, and for each query by position of the start of the alignment. For each such match, an initial output line gives summary information: Example: 440 2.38 1.39 0.79 hh44a1.s1 33 536 ( 0) C 00311 ( 3084) 8277 7771 * Here 440 = smith-waterman score of the match (complexity-adjusted, by default). 2.38 = %substitutions in matching region 1.39 = %deletions (in 1st seq rel to 2d) in matching region 0.79 = %insertions (in 1st seq rel to 2d) in matching region hh44a1.s1 = id of 1st sequence 33 = starting position of match in 1st sequence 536 = ending position of match in 1st sequence (0) = no. of bases in 1st sequence past the ending position of match (so 0 means that the match extended all the way to the end of the 1st sequence) C 00311 : match is with the Complement of sequence 00311 ( 3084) : there are 3084 bases in (complement of) 2d sequence prior to beginning of the match 8277 = starting position of match in 2d sequence (using top-strand numbering) 7771 = ending position of match in 2d sequence * indicates that there is a higher-scoring match whose domain partly includes the domain of this match. N.B. The following only appear if the appropriate option is selected. Following each match there appears a listing of the discrepancies between the aligned portions of the two sequences, giving for each discrepancy its type, position, nucleotide, and quality (in first sequence), and position in second sequence. This list is followed by a table giving, for each quality value, the total number of bases of that quality in the first sequence; the number which are included in the SWAT alignment; the cumulative no. of aligned bases (for this and higher qualities); the number of bases not included in the alignment (but potentially alignable -- i.e. not lying before the beginning or after the end of the second sequence); the number of discrepancies of each type (substitution, deletion, insertion); the total # of discrepancies (for this quality level) and their percentage (of the aligned bases of the given quality), and the cumulative # of discrepancies and their percentage (of the cumulative no. of aligned bases). This information is useful for computing accuracy as a function of quality, for automatically generated contig sequences. In the discrepancy listing and table, the regions at the ends of the sequence that consist entirely of quality 0 are reassigned quality values of -1, since these usually correspond to the lowest quality data at the extreme ends of reads. (Not done in the phrap output...) An 'N' in either sequence is always counted as a substitution error. [N.B. MUCH OF THE FOLLOWING IS NOW OBSOLETE]. Outline of phrap assembly: 0) Read in sequence & quality data, trim off any near-homopolymer runs at ends of reads (since this is likely to be garbage), construct read complements. 1) Find pairs of reads with matching words. Eliminate exact duplicate reads. Do swat comparisons of pairs of reads which have matching words, compute (complexity-adjusted) swat score. 2) Find probable vector matches (mark so they aren't used in assembly). 3) Find near duplicate reads. 4) Find reads with self-matches (probable tandem repeats). 5) Find matching read pairs that are "node-rejected" i.e. do not have "solid" segments. 6) Use pairwise matches to identify confirmed parts of reads; use these to compute revised quality values. 7) Compute LLR scores for each match (based on qualities of discrepant and matching bases). (Iterate above two steps). 8) Find best reads (highest LLR-scores among several overlapping). 9) Identify probable chimeric and deletion reads (withheld from assembly). 10) Construct contig layouts, using (non-rejected) pairwise matches in decreasing score order (greedy algorithm). Consistency of layout is checked at pairwise comparison level. 11) Construct contig sequence as mosaic of highest quality parts of the reads. 12) Align reads to contig; tabulate inconsistencies (read / contig discrepancies) & possible sites of misassembly. Adjust LLR-scores of contig sequence. Description of algorithms: 0) The procedure used in both phrap and cross_match to find sequence matches is the following. First, (in phrap only) any region at the beginning or end of a read that consists almost entirely of a single letter is converted to 'N's; such regions are highly likely to be of poor data quality which if not masked can lead to spurious matches. Reads are then converted to uppercase (in order to allow case-insensitive word matches). All matching words of length at least minmatch between any pair of sequences are then found, by (i) constructing a list of pointers to each position (in each sequence) that begins a word of at least minmatch letters not containing 'N' or 'X'; (ii) sorting the list (using a modified version of quicksort, with string comparison as the comparison function + a few tricks to improve speed by keeping track of the parts of the words that are known already to match); (iii) scanning the sorted list to find pairs of matching words. For each such pair, a band of a specified width (in the imaginary dot matrix for the two sequences), centered on the diagonal defined by the matching words, is defined. Overlapping bands (for the same pair of sequences) are merged. Following construction and merging of bands, a "recursive" SWAT search of each band is used to find matching segments with score greater than or equal to minscore: "recursive" here means that if such a match is found, the corresponding aligned segments in each sequence are (conceptually) X'd out and the process repeated on the remaining portions (if any) of each sequence. This procedure allows (at the cost of some redundant calculation) detection of multiple matching pieces in different locations, and will usually find most copies of repeats (since they generally occur in separate bands). By default, SWAT scores are complexity-adjusted (so that matches for which the collection of matching nucleotides has biassed composition have their scores significantly penalized.) 1) In phrap, sequence pairs with score >= minscore are considered for possible merging. A critical issue here is the appropriate score matrix for SWAT. [Given the generally high accuracy of reads and the desirability of minimizing false matches due to imperfect repeats, a relatively high penalty seems desirable. Currently I use +1 for a match, -9 for a mismatch involving A,C,G, or T, 0 for a match or mismatch involving N, -1 for a match or mismatch involving X, -11 for the gap initiating penalty (the first residue in a gap), and -10 for the gap extension penalty (each subsequent residue).] Setting the indel penalties slightly higher than the mismatch penalties gives better alignments by favoring mismatches in compression regions. Since SWAT uses profiles, one has the option of distinguishing different quality levels by use of different symbols (e.g. upper case for more accurate calls, lower case for less accurate calls) and setting the penalties appropriately; this would involve using the quality levels to adjust the symbols used in the reads, which should be done AFTER the word matching routine. (At present it does not seem particularly useful to use differing positive scores for different nucleotides to reflect their different frequencies). 1a) Determine "confirmed" part of each read (i.e. the part which appears in a SWAT alignment against some other read; a read with the same name -- up to an internal '.' if any -- is not considered confirming). Probable chimeras are detected as reads for which the confirmed part can be separated into two non-overlapping pieces, separated by at most MAX_CHIMERA_GAP bp (currently 30) such that the part confirmed by a given read lies in one or the other piece but not both; AND such that each piece has a prematurely terminating alignment with some other read. Chimeric reads may arise from i) chimeric clones; ii) gel mistracking across lanes; iii) unremoved sequencing vector; (iv) deletion clones. Reads having two non-overlapping confirmed pieces (but failing the "premature termination" condition) are often non-chimeras that are the only link between two non-overlapping contigs, and thus should be permitted in the assembly. Identify "strongly confirmed" regions in each read: having matches to a reverse sense read. Identify deletions. Identify "rejected" alignments -- those which don't extend as far as they should, or that have mismatches involving high quality bases. 2) Sort all matching pairs by decreasing LLR score, and assemble layout by progressively merging pairs with high score. Consistency of merge is required: the SWAT alignment implies relative offset of one contig with respect to another, and hence a potential implied overlap. The SWAT alignment(s) should extend over essentially the entire region of implied overlap, apart from an allowable gap (necessary to allow for the fact that the ends of the alignments are somewhat uncertain). Probable deletion clones (identified as reads which have two adjacent pieces matching two separated pieces of another read, and having no confirming read across the breakpoint), and chimeras are not used in any merges. [N.B. Following is obsolete: Potential chimeras are deferred to a second pass (as well as being flagged in output), so that true reads will assemble first.] 3) Construct contig "consensus" as a mosaic of individual reads. Strategy: ends of alignments, and midpoints of perfectly matching segments of sufficient length, define crosslinks - pursue crosslinks which increase accuracy and extend read. Formally this is done by constructing a weighted directed graph whose nodes consist of (selected) positions in reads; there are bidirectional edges with weight 0 between aligned bases in overlapping reads, and unidirectional edges from 5' to 3' positions within a single read, with weight equal to the total quality of the sequence between the two nodes. Standard C.S. algorithms (dating back to Tarjan) permit identification of a path with maximal weight, in time linear w.r.t. number of nodes. The quality values for the resulting sequence are inherited from the read segments of which it is composed.