/***************************************************************************** # Copyright (C) 1994-1999 by Phil Green. # All rights reserved. # # This software is part of a beta-test version of the swat/cross_match/phrap # package. It should not be redistributed or # used for any commercial purpose, including commercially funded # sequencing, without written permission from the author and the # University of Washington. # # This software is provided ``AS IS'' and any express or implied # warranties, including, but not limited to, the implied warranties of # merchantability and fitness for a particular purpose, are disclaimed. # In no event shall the authors or the University of Washington be # liable for any direct, indirect, incidental, special, exemplary, or # consequential damages (including, but not limited to, procurement of # substitute goods or services; loss of use, data, or profits; or # business interruption) however caused and on any theory of liability, # whether in contract, strict liability, or tort (including negligence # or otherwise) arising in any way out of the use of this software, even # if advised of the possibility of such damage. # #*****************************************************************************/ DOCUMENTATION FOR PHRAP AND CROSS_MATCH (VERSION 0.990319) phrap ("phragment assembly program", or "phil's revised assembly program"; a homonym of "frappe" = French for "swat") -- a program for assembling shotgun DNA sequence data. Key features: allows use of entire read (not just trimmed high 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); provides extensive information about assembly (including quality values for contig sequence) to assist trouble-shooting; able to handle very large datasets. N.B. phrap does not provide editing or viewing capabilities; these are available with consed and phrapview. It is strongly recommended that phrap be used in conjunction with the base calls and base quality values produced by the basecaller, phred; and with the sequence editor/assembly viewer, consed. cross_match -- a general-purpose utility (based on a "banded" version of 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. CONTENTS I. INSTALLATION II. RUNNING PHRAP & CROSS_MATCH 1. Phrap 2. Cross_match III. INPUT DATA FILES 1. Read naming convention 2. Sequence files 3. Quality files 4. Vector screening IV. COMMAND LINE OPTIONS 1. Scoring of pairwise alignments 2. Banded search 3. Filtering of matches 4. Input data interpretation 5. Assembly 6. Consensus sequence construction 7. Output 8. Miscellaneous V. VIEWING PHRAP ASSEMBLIES WITH OTHER PROGRAMS 1. Consed 2. Phrapview 3. Gap VI. PHRAP OUTPUT VII. CROSS_MATCH OUTPUT VIII. EST ASSEMBLIES (TO BE ADDED) IX. PROBLEMS 1. Insufficient memory 2. Other phrap- or cross_match-generated error messages (TO BE ADDED) 3. Phrap- or cross_match-generated warning messages (TO BE ADDED) 4. "Crashes" reported by operating system 5. Long running time 6. Misassemblies, incomplete assemblies, incorrect consensus sequence (TO BE ADDED) 7. Reporting problems APPENDIX: ALGORITHMS I. INSTALLATION The source code for the swat/cross_match/phrap package will be sent to you in the form an email message containing a uuencoded .tar.Z file; you will need to have access to a Unix system for the initial unpacking, but once you've uudecoded it and unpacked the .tar file (steps i and ii below), you should be able to compile the programs on computers running other operating systems -- they should be portable to almost anything with a decent C compiler and adequate memory (64 Mb RAM or more is desirable). Here are the steps needed to unpack and install the programs: i. Save the email message as a file (for example, "temp.mail"). If possible, do this using the Unix mail command, rather than another mail program -- some mail programs (e.g. Pine) remove trailing spaces on each line of incoming messages, which will corrupt a uuencoded message. Do not attempt to modify the saved mail message in any way. That is unnecessary and may corrupt the message. ii. To unpack the saved email message, execute the following two commands on a Unix workstation, in the directory containing the file created in step i above: > uudecode temp.mail > zcat distrib.tar.Z | tar xvf - If either of these commands results in an error message, it is likely that the email message was corrupted by your mail program (see step i above). iii. To produce working versions of the programs, move (if necessary) all of the files produced by the above command to the computer on which you wish to run the programs (which must have a C compiler!), and execute the following command in the directory that contains the files: > make If your compiler does not recognize the -O2 optimization flag (which should be evident from warning messages it produces), you should change the line CFLAGS= -O2 in the file "makefile" to CFLAGS= -O Then remove all files ending in .o produced by the original make, and recompile. iv. If you have datasets with more than 64,000 reads, or that include sequences longer than 64,000 bp, you will need the .manyreads and/or .longreads versions of phrap and cross_match. These are created using the command > make manyreads v. If you are operating a non-commercial (academic or government) computer facility which provides access to several independent investigators, you are required by the licensing agreement to set the permissions on the executables and source code to allow execute but not read access, so that the programs may not be copied. II. RUNNING PHRAP & CROSS_MATCH 1. Phrap. i. PhredPhrap script. Generally if you are assembling reads basecalled by phred you should run phrap via the phredPhrap script rather than directly. PhredPhrap runs phred; runs another script, phd2fasta, to make the fasta format sequence file from the .phd files generated by phred; runs cross_match to mask vector sequence; runs phrap; and finally runs some other programs that provide useful information to consed (e.g. tagging of repeats). It will accept any of the phrap command line options described below (section IV). PhredPhrap and its documentation are included in the consed distribution. Although phredPhrap automates most of the steps in the assembly, you will still need to do several things: (a) set up the directory structure that is assumed by phredPhrap and make sure your chromatogram files are appropriately located in this structure; (b) adhere to the read naming convention assumed by phrap, or make other arrangements to get template/chemistry/read direction information into the .phd or fasta file (THIS IS IMPORTANT -- SEE section III.1 below); and (c) make sure that appropriate clone and subclone (sequencing) vector sequences are included in a fasta file used for screening with cross_match (see section III.4 below). ii. The following instructions pertain to "standalone" operation of phrap. Before the run, you may want to increase the amount of memory available to phrap: see the instructions for using the Unix "limit" command under IX.1, "Insufficient memory" below. The command line should be of the form > phrap seq_file1 [seq_file2 ...] [-option value] [-option value] ... where each seq_file is the name of a sequence file in FASTA format as described below in III.2. Available options are described in section IV below. The standard output is extensive and should be redirected to a file (phrap.out in the following example). Example: > phrap C05D11.reads.screen -minmatch 20 -new_ace > phrap.out Note that the quality file is not specified on the command line, but instead must be named appropriately (as described below in III.3) to be recognized by phrap; in this example its name would need to be C05D11.reads.screen.qual. If more than one seq_file 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. [N.B. THE ABILITY TO USE MULTIPLE SEQUENCE FILES IS NOT AVAILABLE AT PRESENT -- YOU SHOULD PROVIDE ONLY A SINGLE SEQUENCE FILE TO PHRAP]. 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. 2. Cross_match. The command line should be of the form > cross_match seq_file1 [seq_file2 ...] [-option value] [-option value] ... where each seq_file is the name of a sequence file in FASTA format as described below in III.2. Available options are described in section IV below. The standard output should generally be redirected to a file. 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 entries are not compared to their own complements, but are compared to themselves in the same orientation, and to all other entries in both orientations). Matches meeting relevant criteria (see section IV) are written to the standard output. III. INPUT FILES Input files to cross_match and phrap are of two types: sequence files (required), and quality files (optional, but strongly recommended for phrap). If you are doing phrap assemblies of reads generated by the basecalling program phred, then it is recommended that you generate these input files automatically using the phd2fasta or phredPhrap scripts distributed with phred and consed (see II.1.i above). Even in this case however you should read the following section on the read naming convention, since you will either need to make sure that the read names attached to your chromatogram files adhere to the naming convention, or you will need to modify the phredPhrap script (or create your own script) so that correct template names and read orientation information is appended to the .phd files after they are produced by phred. 1. Read naming convention (phrap only) In addition to the sequence and quality data, phrap needs to know three things for each read: (i) the name of the subclone (or other template) from which the read is derived; this is used, for example, in checking for chimeric subclones (for which purpose it is necessary to know when two reads are from the same subclone so that they are not regarded as independently confirming each other) and for certain other data anomalies; (ii) the orientation of the read (forward or reverse) within the subclone, in cases where data is acquired from each end of a subclone insert (at present this information is used only for consistency checking following assembly, and not in the assembly itself, but this will change in future versions); (iii) the chemistry used to generate the read (which influences phrap's decisions as to how to treat discrepancies between potentially overlapping reads, and as to how to adjust qualities of reads which confirm each other -- confirmation of a read by another read with different chemistry counts as significantly as confirmation by an opposite-strand read). The above information is generally conveyed via the read names using the "St. Louis" read naming convention described below. If your laboratory uses a different naming system, you have several options: (i) provide the relevant information in the description field of the .fasta file instead (see below), or (if you are using the phd2Fasta script) as tags in the .phd files for each read; (ii) see whether appropriate values of the command line options -subclone_delim and -n_delim (described below) can be used to bring your own naming convention into conformity with the St. Louis convention -- if your naming convention is such that the subclone name always occurs as the first part of the read name, and the end of the subclone name occurs at the first (or n-th for some fixed value of n) occurrence of a particular character, then you may be able to use these options without having to change your read names; (iii) create a script which translates your read names into St. Louis form (this does not mean that you cannot retain the original names as well -- for example, for purposes of running phrap you could create a directory of "links" to your chromatogram files, such that the link names adhere to the naming convention, and then run the phredPhrap script on this directory); (iv) alter the subroutines in the source module "names.c" in order to parse your read names correctly; or (v) ignore the issue and hope that it does not cause problems. This may have unpredictable adverse effects on assembly and is NOT recommended. The "St. Louis" read naming convention assumed by phrap is as follows: The portion of the read name up to the first '.', if any, is the name of the subclone from which the read is derived (i.e. the DNA template used in the sequencing reaction). If desired, the command line option -subclone_delim can be used to change the character which delimits the end of the subclone name to something other than '.', and/or the option -n_delim may be used to indicate that the subclone name contains a fixed number of occurrences of this character within it. See section IV below for a description of command line options. The orientation of the read within the subclone, and the chemistry, are indicated by the first letter following the '.', as follows: "s" forward direction read on single stranded (SS) template, dye primer chemistry "f" forward read on double stranded (DS) template, dye primer chemistry "r" DS reverse read, dye primer chemistry "x" SS forward read, standard dye terminator chemistry "z" DS forward read, standard dye terminator chemistry "y" DS reverse read, standard dye terminator chemistry "i" SS forward read, big dye terminator chemistry "b" DS forward read, big dye terminator chemistry "g" DS reverse read, big dye terminator chemistry "t" for T7 (cDNAs) "p" for SP6 (cDNAs) "e" for T3 (cDNAs) "d" for special "c" consensus pieces "a" assembly pieces The remainder of the name is ignored. If the read name ends in ".seq" the .seq is removed. Example: the name BAC112_a11b3.x3_pg indicates a forward read with (old) dye terminator chemistry from the subclone BAC112_a11b3. 2. Sequence files. With both phrap and cross_match, either a single sequence file or multiple sequence files can be specified (N.B. WITH CURRENT VERSION OF PHRAP ONLY A SINGLE SEQUENCE FILE SHOULD 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). Cross_match is typically used with two input sequence files, the first containing the "query" sequences and the second containing the "subject" sequences, but it may also be used to compare the sequences in a single input file to each other. 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 when available. If phred-generated qualities are not available however then the default quality parameter -default_qual should be set appropriately (section IV). Sequence files must be in FASTA format: (i) Each sequence entry has a single header line as follows: The first character is '>'. This is followed immediately (i.e. with no intervening spaces or other characters) by the sequence name, and then (optionally) by a space or spaces, followed by descriptive information about the sequence. (ii) The header line is followed by a separate line or lines containing the sequence. The sequence names, descriptive information, and sequences can be of arbitrary length (up to a combined total size for each type of about 2 billion characters), since space for them is allocated dynamically. If there are more than 64,000 entries in the first (query) file, or if any of the sequences are longer than about 64,000 characters, the .manyreads version of phrap or cross_match must be used. Descriptive information on the header line may optionally be used to indicate template (subclone), read orientation, read chemistry and dye. If present, this information overrides whatever would be implied by the read name. The template name is indicated by including the word "TEMPLATE:", followed by a space, followed by the name of the template. Orientation is indicated by the word "DIRECTION:", followed by the word "fwd" or "rev". Chemistry is indicated by the word "CHEM:", followed by the word "prim" (for dye primer), "term" for dye terminator, or "unknown". Dye is indicated by the word "DYE:" followed by "rhod", "big", "ET", "d-rhod", or "unknown". The chemistry and dye information is provided automatically by newer versions of phred (version 0.980904.a or later), which obtain it from the chromatogram file. (Use the script phd2fasta distributed with phred and consed to create the .fasta file from the .phd files created by phred). Example: >read#5 DIRECTION: rev CHEM: prim DYE: ET TEMPLATE: a23f1 any_other_information GAAAGATCTCATTGATCACTCTATTCAAGTGGGAGTCTCCGGTCTTT ATGATCGATTTGTGAATCTTCGTATCAAAGTTGGAGCTGACAAGTATCCA TTGCTTGCGAAATGGGCTCAAATTTTCACTCAGGGAGTCGTCTTCGATCC TTCAAGAATTCATCAATGTGCTCAAAAGTTGGCTGGAGAAGCTCGTGATC GGAAGAGAGATGGATGTACTGTGGCATCAACTGCAGTAGCTTCAATGGTT TATGGAAAGAGTATGTTATTt All descriptive information in the header line is written back out to the .ace and .singlets files. Non-alphabetic characters (including '*', and digits) 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 not possible (as it was in old phrap versions) to use case as a crude indication of quality. 3. Quality files. For each input sequence file in FASTA format, you may optionally include a corresponding file of data quality information. This is strongly recommended for phrap runs since it greatly improves the accuracy of assembly and of the consensus sequence; with cross_match, the qualities are at present used only for purposes of annotating and tabulating discrepancies and not in the alignment scoring. The name of this file must 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 each read must match the number of bases for that read in the FASTA file. Also, the orders of the reads in the FASTA and corresponding .qual file must be identical. N's are automatically assigned quality 0, overriding any value present in the .qual file; however X's retain whatever value is present in the .qual file. 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; however when it is available, it is important to provide input data quality information since this can substantially increase the accuracy of the consensus, particularly in regions where there are strand-specific effects on quality (e.g. compressions) -- 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 automated sequencing machines it is strongly recommended 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 99 is 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 low quality, 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 defined in terms of 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. 4. 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 before running phrap. (It is unnecessary to remove the cloning (e.g. cosmid or BAC) vector, unless the inserts from multiple overlapping clones are being assembled simultaneously; however it is generally useful to remove at least the central part of the cloning vector, since this then provides a natural point at which phrap can break the circular 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 as described here, whether or not you have already screened them using another program. To carry out the screening, 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 10 -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 is what should be provided as input to phrap. 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. The phredPhrap script automatically carries out the above steps (i.e. vector screening using cross_match, and renaming of the .qual file); however it is your responsibility to make sure the vector sequence file includes vector sequences appropriate for your cloning and sequencing vectors. IV. COMMAND LINE OPTIONS This section describes all command line options available with cross_match and phrap, grouped by category. Each option name is followed immediately by the "default value" assumed by the program, and a brief description. A '*' appearing as the default value indicates that the option is a flag rather than a parameter; in this case no value should be included on the command line after the option name, and by default it is turned off. [None] means there is no default, i.e. the parameter is not relevant when the option is omitted. Unless otherwise indicated, an option can be used with either cross_match or phrap. 1. Scoring of pairwise alignments The following options control the scoring of pairwise sequence alignments. By default, matching residues receive a reward of +1, mismatches get a penalty of -2, gap opening residues a penalty of -4 (= -2 - 2), and gap extension residues a penalty of -3 (= -2 - 1). The highest scoring word-nucleated local alignment between a pair of sequences is found, and its score is then complexity-adjusted (i.e. scores of matches between sequence regions of biassed nucleotide composition are adjusted downwards). These parameters are chosen such that regions of two sequences that are about 70% or more identical will tend to have a positive alignment score. For more stringent comparisons, -penalty should be made more negative (e.g. a value of -9 will tend to find alignments that are 90% identical), and/or -minscore should be increased; for more sensitive comparisons, a score matrix should be used instead. option name & default value -penalty -2 Mismatch (substitution) penalty for SWAT comparisons. -gap_init penalty-2 Gap initiation penalty for SWAT comparisons. -gap_ext penalty-1 Gap extension penalty for SWAT comparisons. -ins_gap_ext gap_ext Insertion gap extension penalty for SWAT comparisons (insertion in subject relative to query). -del_gap_ext gap_ext Deletion gap extension penalty for SWAT comparisons (deletion in subject relative to query). -matrix [None] Score matrix for SWAT comparisons (if present, supersedes -penalty) Matrix format: (TO BE ADDED) -raw * Use raw rather than complexity-adjusted Smith-Waterman scores. 2. Banded search The following options control the region of the Smith-Waterman matrix that is searched. In phrap and cross_match (as opposed to swat, which performs full Smith-Waterman comparisons), word-nucleated banded Smith-Waterman comparisons are performed: the set of input sequences is scanned to find pairs of perfectly matching subsequences ("word matches") meeting certain criteria, and for each such word match a band in the Smith-Waterman matrix that is centered on the diagonal defined by the match is searched to find an optimal-scoring alignment. If there are multiple matching words for a given pair of sequences, the union of the corresponding bands is searched. The search is "recursive" in the sense that if an optimal-scoring alignment is found whose score exceeds -minscore the remainder of the band is searched again; as a result multiple alignments may be found within a single band. The word matches that are considered are "maximal" in the sense that they cannot be extended in either direction without introducing a discrepancy in one of the two sequences. If the length of the word match is at least -maxmatch, it is accepted; otherwise if the complexity-adjusted length is less than -minmatch, or if the matching sequence occurs more than -max_group_size times on the forward strands of the query file, it is rejected. The latter criterion effectively downweights frequently occurring words. It can be turned off by setting -maxmatch equal to -minmatch. -minmatch 14 Minimum length of matching word to nucleate SWAT comparison. if minmatch = 0, a full (non-banded) comparison is done [N.B. NOT PERMITTED IN CURRENT VERSION]. Increasing -minmatch can dramatically decrease the time required for the pairwise sequence comparisons; in phrap, it also tends to have the effect of increasing assembly stringency. However it may cause some significant matches to be missed, and it may increase the risk of incorrect joins in phrap in certain situations (by causing implied overlaps between reads with high-quality discrepancies to be missed). -maxmatch 30 Maximum length of matching word. For cross_match, the default value is equal to minmatch, instead of 30. -max_group_size 20 Group size (query file, forward strand words) -word_raw * Use raw rather than complexity-adjusted word length, in testing against minmatch (N.B. maxmatch always refer to raw lengths). (The default is to adjust word length to reflect complexity of matching sequence). -bandwidth 14 1/2 band width for banded SWAT searches (full width is 2 times bandwidth + 1). Decreasing bandwidth also decreases running time at the expense of sensitivity. Phrap assemblies of clones containing long tandem repeats of a short repeat unit (< 30 bp) may be more accurately assembled by decreasing -bandwidth; -bandwidth should be set such that 2 bandwidth + 1 is less than the length of a repeat unit. -bandwidth 0 can be used to find gap-free alignments. 3. Filtering of matches The following options control which matches are displayed (cross_match) or used in assembly (phrap). -minscore 30 Minimum alignment score. -vector_bound 80 Number of potential vector bases at beginning of each read. Matches that lie entirely within this region are assumed to represent vector matches and are ignored. For cross_match, the default value is 0 instead of 80. -masklevel 80 (cross_match only). A match is reported only if at least (100 - masklevel)% of the bases in its "domain" (the part of the query that is aligned) are not contained within the domain of any higher-scoring match. Special cases: -masklevel 0 report only the single highest scoring match for each query -masklevel 100 report any match whose domain is not completely contained within a higher scoring match -masklevel 101 report all matches 4. Input data interpretation The following options influence how input data (including read names) are interpreted. -default_qual 15 Quality value to be used for each base, when no input .qual file is provided. Note that a quality value of 15 corresponds to an error rate of approximately 1 in 30 bases, i.e. relatively accurate sequence. If you are using sequence that is substantially less accurate than this and do not have phred-generated quality values you should be sure to decrease the value of this parameter. -subclone_delim . (phrap only). Subclone name delimiter: Character used to indicate end of that part of the read name that corresponds to the subclone name -n_delim 1 (phrap only). Indicates which occurrence of the subclone delimiter character denotes the end of the subclone name (so for example -subclone_delim _ -n_delim 2 means that the end of the subclone name occurs at the second occurrence of the character '_'). Must be the same for all reads! -group_delim _ (phrap only). Group name delimiter: Character used to indicate end of that part of the read name that corresponds to the group name (relevant only if option -preassemble is used); this character must occur before the subclone delimiter (else it has no effect, and the read is not assigned to a group). -trim_start 0 (phrap only). No. of bases to be removed at beginning of each read. 5. Assembly The following options are used to control completeness and stringency of assembly; note that the options -minmatch, -bandwidth, -penalty, and -minscore discussed above also affect assembly stringency. -forcelevel 0 (phrap only). Relaxes stringency to varying degree during final contig merge pass. Allowed values are integers from 0 (most stringent) to 10 (least stringent), inclusive. -bypasslevel 1 (phrap only). Controls treatment of inconsistent reads in merge. Currently allowed values are 0 (no bypasses allowed; most stringent) and 1 (a single conflicting read may be bypassed). -maxgap 30 (phrap only). Maximum permitted size of an unmatched region in merging contigs, during first (most stringent) merging pass. -repeat_stringency .95 (phrap only). Controls stringency of match required for joins. Must be less than 1 (highest stringency), and greater than 0 (lowest stringency). -revise_greedy * (phrap only). Splits initial greedy assembly into pieces at "weak joins", and then tries to reattach them to give higher overall score. Use of this option should correct some types of missassembly. -shatter_greedy * (phrap only). Breaks assembly at weak joins (as with -revise_greedy) but does not try to reattach pieces. -preassemble * (phrap only). Preassemble reads within groups, prior to merging with other groups. This is useful for example when the input data set consists of reads from two distinct but overlapping clones, and it is desired to assemble the reads from each clone separately before merging in order to reduce the risk of incorrect joins due to repeats. The preassemble merging pass is relatively stringent and not guaranteed to merge all of the reads from a group. Groups are indicated by the first part of the read name, up to the character specified by -group_delim. -force_high * (phrap only). Causes edited high-quality discrepancies to be ignored during final contig merge pass. This option may be useful when it is suspected that incorrect edits are causing a misassembly. 6. Consensus sequence construction The following options affect the weighted directed graph that is used to find the consensus sequence. Increasing their values will reduce the size of the graph, which should reduce memory requirements for the phrap run (substantially in some cases) but may decrease the accuracy of the consensus sequence that is found. -node_seg 8 (phrap only). Minimum segment size (for purposes of traversing weighted directed graph). -node_space 4 (phrap only). Spacing between nodes (in weighted directed graph) . 7. Output The following options control the creation or formatting of certain output files, and of the standard output. -tags * Tag selected lines in the standard output, to facilitate parsing. (N.B. this does NOT refer to the phrap-generated tags that are included in the .ace file and viewable in consed.) -screen * (cross_match): Create a ".screen" file. This is a FASTA format sequence file, in which any sequence region in a sequence in the first input sequence file that matches some region in a sequence in the second (or later) file is replaced by X's. This option is primarily used to create "vector-masked" copies of the reads prior to phrap assembly. (phrap): when the -old_ace or -new_ace option is specified (see below), this option causes parts of the read sequences that consist of phrap-inferred sequencing vector and chimeric segments to be replaced by X's in the .ace file. (The "phrap-inferred sequencing vector" consists of the beginning part of the read that either does not match any other read, or matches only the beginning parts of other reads; "phrap- inferred chimeric segments" consist of segments of the read that match other reads but do not match the consensus at that location.) -old_ace * (phrap only). Create ".ace" file in old style format. -new_ace * (phrap only). Create ".ace" file in a new style format (STRONGLY recommended over old_ace !!) -ace * (phrap only). Same as -new_ace. -view * (phrap only). Create ".view" file suitable for input to phrapview. -qual_show 20 (phrap only). Cutoff for flagging "low_quality" regions in contig sequence and "high quality" discrepancies between read and contig. Bases in the .contigs and .ace file are lower case if and only if their LLR-converted quality values are below this value. -print_extraneous_matches * (phrap only). Print information about non-local matches between contigs. -exp [None] (gcphrap only). Name of a directory in which output experiment files are to be placed. -alignments * (cross_match only). Display the alignment for each match. -discrep_lists * (cross_match only). Give list of discrepancies for each match. -discrep_tables * (cross_match only). Give table of discrepancies (by quality) for each match (default is to display a single table, that combines results for all matches). 8. Miscellaneous -retain_duplicates * (phrap only). Retain exact duplicate reads, rather than eliminating them. -max_subclone_size 5000 (phrap only). Maximum subclone size -- for forward-reverse read pair consistency checks. -trim_penalty -2 (phrap only). Penalty used for identifying degenerate sequence at beginning & end of read. -trim_score 20 (phrap only). Minimum score for identifying degenerate sequence at beginning & end of read. -trim_qual 13 (phrap only). Quality value used in to define the "high-quality" part of a read, (the part which should overlap; this is used to adjust qualities at ends of reads. -confirm_length 8 (phrap only). Minimum size of confirming segment (segment starts at 3d distinct nuc following discrepancy). -confirm_trim 1 (phrap only). Amount by which confirming segments are trimmed at edges. -confirm_penalty -5 (phrap only). Penalty used in aligning against "confirming" reads. -confirm_score 30 (phrap only). Minimum alignment score for a read to be allowed to "confirm" part of another read. -indexwordsize 10 Size of indexing (hashing) words, used in finding word matches between sequences. The value of this parameter has a generally minor effect on run time and memory usage. V. VIEWING PHRAP ASSEMBLIES WITH OTHER PROGRAMS 1. Consed Consed is the recommended sequence editor for use with phrap. Its development has been closely tied to the development of phrap and it has been designed to take advantage of the information regarding the assembly that is generated by phrap (e.g. phrap's consensus sequence and quality values, and tags indicating potential assembly problems and various data anomalies such as chimeras and compressions). To use consed with phrap, ace files must be generated during the phrap run using the flag -new_ace or -old_ace on the command line. This is done automatically if you use the script phredPhrap to carry out all of the data processing steps. For additional information, consult consed documentation. 2. Phrapview Phrapview is a graphical tool that provides a "global" view of the phrap assembly, complementary to the "local" view provided by the sequence editor CONSED; it will become obsolete when a similar globabl view is added to CONSED. a. Installation Phrapview is written in perl/Tk; to run it you will need to have installed on your system a recent version of perl (perl5 or later), together with the perl/Tk perl module written by Nick Ing-Simmons. The perl/Tk scripts can be obtained from any CPAN site in modules/by-authors/Nick_Ing-Simmons/, for example: ftp://ftp.perl.org/pub/CPAN/modules/by-authors/Nick_Ing-Simmons/ The latest version of perl/Tk is the most recent Tk800.xxx.tar.gz file (where xxx is the most recent version number). The latest version of perl can also be obtained from any CPAN site in /src, for example ftp://ftp.perl.org/pub/CPAN/src/perl5.005_02.tar.gz (Perl 5.005_02 was the most recent version when this was written.) If you have problems with running phrapview, please make sure that you have both a recent version of perl installed and a recent version of perl/Tk. If you are unsure whether your version is recent enough, you or your system administrator should download the latest sources and install them. b. Running phrapview Phrapview is invoked with the command phrapview filename where filename is the name of a ".view" file that was created by running phrap with the -view option. N.B. The format of the .view file is still under development. To ensure that phrapview will read correctly the .view file produced by phrap, make sure that both programs (phrap and phrapview) were part of the same release. The following information can be displayed: (i) Basic information concerning the numbers of reads, singletons, contigs, chimeras, etc. Each contig, and each chimeric read, is represented by a horizontal black line proportional in length to the contig or read size. Note that although they are shown separately in this display, chimeric reads in fact generally are incorporated by phrap into one of the contigs, in a region that corresponds to one of the two chimeric pieces; the location of this is indicated by the black "chimera match" line connecting the chimera to a contig (see below). The following information requires pressing an appropriate "button" to display: (ii) Depth of coverage. Graphs (in yellow) above and below the contig line indicate the depth of coverage (i.e. number of reads aligned against the contig) on each strand at each position; the horizontal orange lines correspond to a depth of 10. In addition to the ordinary depth of coverage (which counts all reads in the contig, and is mainly useful to indicate regions where an abnormal deficit or excess occurs), the "reduced depth" can be displayed. This counts only reads that are non-chimeric, have a positive LLR score against the contig, and have no positive LLR scoring match to a read elsewhere in the assembly -- i.e. the "confidently placed" reads. (See phrap.doc for a description of LLR scores). Misjoins usually occur in places where the reduced depth is 0 on both strands. Regions of the contig where the depth is 0 on one or both strands are indicated in red. (iii) "Matches" and "links". "Matches" are pairwise read matches between reads in different parts of the assembly, and are indicated by a single (curved or straight) line connecting the midpoints of the two matching regions. Contig matches (between two non-chimeric reads in the assembled contigs) and chimera matches (between a chimeric read and itself or another read) can be displayed separately. "Links" connect the startpoints of two reads that derive from the same subclone: forward-reverse links correspond to reads from opposite ends of a subclone insert, while same-strand links correspond to walking or duplicated reads in the same direction on the insert. Links will only be indicated properly if the read naming convention assumed by phrap is used (see section III.1 above). Each match or link is considered to be in one of three classes, indicated by color: "problems" (red), which indicate serious discrepancies or a significant possibility of incorrect, incomplete, or non-unique assembly; "ok" (black), which tend to confirm the assembly or are otherwise consistent with it; and "grayzone" (blue or green), which may in some cases indicate problems but are probably ok. Links are considered "ok" if the two read starts meet the expected spacing and orientation constraints, otherwise they are marked as "problems" (in many cases however these reflect subclone tracking errors rather than assembly errors). Matches are considered "ok", and are not shown, if they are between reads assembled in the same location, or have LLR scores below the cutoff. Matches are considered "problems" if they have LLR scores above the cutoff, are non-local, and involve regions where the reduced depth is 0. "Problem" matches between or within large contigs often indicate the presence of near-perfect repeats, which may or may not have been misassembled; or (in some cases) a join that was missed. "Problem" matches between a small (e.g. singleton) contig and a large one typically represent cases where a read could not be assembled into its proper location by phrap, due generally to some anomaly involving it (e.g. an internal region of low quality data). Non-local matches with LLR scores above the cutoff which do not lie in regions of reduced depth 0 are considered to be "grayzone". Usually these involve isolated pairs of reads each of which extends part way into a repeat, such that the overlapping region is too small to contain high-quality discrepancies. The fact that the reduced depth is non-zero means that there are other reads in the same region which do not have any positive LLR scores to the other region, so it is unlikely that there is a significant misassembly in such cases, although it is possible that one of the two matching reads may have been incorrectly placed. Chimeras generally have one "ok" match (to the site where the read was assembled) and one or more "problem" matches to the region from which the other part of the chimera derives (sometimes these matches all have negative LLR scores and are thus not shown with the default LLR cutoff). Some chimeras actually represent subclones with internal deletions; these tend to be obvious since the two pieces map to nearby locations in some contig and are in the same orientation there. Double-headed arrows on matches indicate that the matching regions are on opposite strands; double-headed arrows on links indicate that the orientation of the two reads is inconsistent with the read names (i.e. forward-reverse pairs that point in the same direction, or same strand pairs that point in opposite directions; these are colored as "problems" if the reads are in the same contig). (iv) Low quality contig bases (bases having phrap qualities that are lower than a specified threshold), and high quality discrepancies (positions where there is an aligned read discrepant with the contig base and having a phrap quality that exceeds the threshold). A horizontal green line above the contig indicates the quality threshhold value ("qual cutoff"). Low quality contig bases are indicated by a vertical blue line drawn from a height corresponding to the contig quality value, up to the threshhold line; high quality discrepancies are indicated by a vertical blue line drawn from a height corresponding to the discrepant read's quality, down to the threshhold line. Thus in each case longer lines are more "serious" (reflect a larger deviation from the threshhold, in the wrong direction). Passing the mouse pointer (slowly!) over an item of types ii-iv above results in its being highlighted (converted to yellow); it remains highlighted until the pointer is moved over another item, or a contig line. Textual information about the highlighted item is shown in a box in the top right region of the display. Several parameters can be changed by entering new values in appropriate boxes at the bottom of the display. Four different magnifications can be changed: horizontal magnification; the "spacing magnification" which determines the vertical spacing between contigs; and the depth and quality magnifications, which affect the scales used for coverage depth and quality. Values of each of these are specified as percents of the original (startup) magnification. The role of LLR cutoff, and qual cutoff are described above. The unalign parameter refers to the size of the largest segment involved in the pairwise read match that is unaligned against the contig sequence. If large, such segments may indicate a misplaced read or (occasionally) an otherwise undetected chimera. Min fwd-rev and max fwd-rev refer to the minimum and maximum allowed separation (in bp) between (the starts of) forward-reverse read pairs; i.e. these represent the minimum and maximum allowed subclone insert size. Separations outside this range (or in the wrong orientation) are marked as problems. Min ss and max ss are the minimum and maximum spacing between same-strand reads (i.e. size of a walking step). The display is updated to reflect changes made in the above parameters only after a relevant button (Show Contig Matches, etc.) is pressed. Colors can be changed by altering the variable definitions that occur near the beginning of the phrapview source file, and re-running the program. 3. Gap A version of phrap ("gcphrap", for "gap-compatible phrap") suitable for use with the Staden gap4 package may be created as follows: Obtain relevant source files from the Staden web site (http://www.mrc-lmb.cam.ac.uk/pubseq/index.html). Put these in the same directory as the phrap source files, and enter the command > make gcphrap Consult the Staden web site for instructions on using gcphrap in conjunction with gap4, including input and output file formats. In addition to the usual phrap command line options, gcphrap accepts an option -exp, which is used to specify the name of a directory in which the output experiment files will be placed. VI. PHRAP OUTPUT Phrap output includes the "standard output", the "standard error", and a series of output files. The latter receive names that consist of the input read file name, with an appropriate suffix (e.g. ".ace") appended. The standard output is described below. The other files are (i) The standard error (which is normally written to the terminal, but may be redirected to a file). This contains information indicating the point in the run that phrap has reached, summary results for some of the steps, and various warning and error messages. (ii) The .contigs file. This is a FASTA file containing the contig sequences (without pads; bases in this file are upper case if and only if the 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. (iii) The .contigs.qual file. This has the phrap-generated qualities for the contig bases. (iv) The .singlets file. A FASTA file containing the singlet reads (i.e. the reads with no match to any other read). (v) The .log file. This includes various diagnostic information, and a summary of aspects of the assembly; it is probably only of interest to me, for trouble-shooting purposes. (vi) The .problems file. Again this is probably only of interest to me. (vii) The .ace file (produced when the options -new_ace or -old_ace is used). This file is required for viewing the assembly using consed. It's format is described in the consed documentation. (vii) The .view file (produced when the option -view is used). Required for viewing the assembly using phrapview. Standard output: The standard output has the goal of giving as complete a summary as possible of the final assembly, including data anomalies and possible sites of assembly error. The information it provides includes a list of the contigs and the reads they contain, information about the quality of the alignments of the reads against the contig sequence, matching regions within and between contigs, suspect reads (probable deletions or chimeras), sites of possible errors in the assembly or contig sequence (e.g. discrepancies between reads and contig sequence), and results of the forward/reverse read consistency checks. The output includes, in order: (i) Summary information about the run conditions (command line, phrap version number, and parameter settings) and the input data file(s) (sequence and quality, template/read counts and read orientations, chemistry). It is worth reviewing this information routinely for possible problems with the input data, such as parameter misspecification, read nomenclature that is causing faulty inferences by phrap regarding templates or chemistry, or a missing quality file. (ii) Low-quality regions at beginnings or ends of reads which consist largely (>50 %) of a single nucleotide, due (usually) 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 .ace and .singlets files output by phrap will contain the N's in place of the original sequence). (iii) Exact and near-exact duplicate reads. Exact duplicates, which probably represent duplicate entry of the same data, are excluded from assembly (this feature may be turned off using the option -retain_duplicates). The near-exact duplicates are only listed if they are from different subclones (insofar as this information is available to phrap!). (iv) Probable unremoved sequencing vector, detected as regions at beginnings of reads which match only the beginnings of other reads. Phrap's method for identifying these is not foolproof and should not be considered a substitute for screening out sequencing vector as described above (III.4). If vector has been removed in that manner, the "unremoved vector" detected here 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. You may want to consider adding 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. I would appreciate hearing of any examples of sequence that is not vector being erroneously labelled as such by phrap. (v) "Internal read matches", i.e. reads which have off-diagonal same orientation matches to themselves (due e.g. to short tandem repeats). 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 short microsatellites may not be detected here. Node-rejected pairs. (vi) Multi-segment reads, defined as reads whose confirmed parts break into two or more pieces. These include chimeras and non-chimeric "single links". They are omitted from the initial merging passes to reduce the risk of false joins in the assembly, but are incorporated in later merging passes; as a result, chimeric reads will generally be assembled into a contig, but should be tagged as chimeric. (vii) 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. The probable deletions are not used in assembly and instead appear as singleton contigs. Revised quality (= phrap quality), and LLR score histograms, for two passes. (viii) 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). (ix) (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 clone 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). Blocked reads # reads lacking a high-quality segment Bypassed reads (x) Isolated singlets. These are reads having no non-vector match (with score at least minscore) to any other read. Shown for each read are its length, and the number of high-quality non-X bases. For most reads, the latter number should be 0, indicating that the read is either entirely vector (its high-quality part has been completely X'd out) or is of very low quality; reads for which this value is positive may be contaminants. (xi) Contigs (including "singleton contigs" consisting of a single read -- these are cases in which the read did have a match with some other read(s), but could not be consistently assembled with it). The contigs are sorted by (increasing) number of reads, so that the singleton 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 including vector and low-quality data, not just the part that is alignable) with respect to the contig 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-singleton contig matching an anomalous (chimeric or deleted) read are indicated only in the output for the singleton contig that contains the anomalous read. (xii) 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. VII. CROSS_MATCH OUTPUT In addition to the standard output, the files produced include the standard error and .log files (as for phrap); and the .screen file (produced when the option -screen is used). Standard output The standard output lists matches between any sequence in the first input file (the "query" sequences) and a sequence in the second (or later) input file (the "subject" sequences); or, if a single input file is provided, the matches between any two sequences in this file. The matches that are reported are controlled by the command line options -minscore and -masklevel, as well as by the options that control scoring of the alignments and the band search (see section IV above). The reported matches are ordered by query, and for each query by the position of the start of the alignment within the query. For each reported 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 * Interpretation: 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. Following this line there appears (if the flag -discrep_lists is specified) a listing of the discrepancies between the aligned regions of the two sequences, giving for each discrepancy its type, position, nucleotide, and quality (in the query sequence), and its position in second sequence. This list is followed (if -discrep_tables is specified) 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. If the flag -alignments is specified, a complete alignment for each reported match is displayed. VIII. EST ASSEMBLIES (TO BE ADDED) IX. PROBLEMS 1. Insufficient memory. If the run stops prematurely, displaying the message FATAL ERROR: REQUESTED MEMORY UNAVAILABLE then you need to increase the amount of memory alloted to you by the operating system. This can usually be done without obtaining additional physical RAM for your computer -- it is usually enough just to have the operating system allocate additional virtual memory to your process. (However if you are substantially exceeding physical RAM the run may take an exceptionally long time to complete). On Unix systems this can be done using the "limit" command, as follows: i. Run "limit" with no arguments to list the system resources currently available to you. For example (this output is for a DEC alpha computer -- the output on other computers may look somewhat different): > limit cputime unlimited filesize unlimited datasize 131072 kbytes stacksize 2048 kbytes coredumpsize unlimited memoryuse 699392 kbytes descriptors 4096 files addressspace 1048576 kbytes The relevant parameter here is datasize, currently set to about 131 megabytes. ii. Now run limit with the option -h to see what the maximum possible values for these parameters are: > limit -h cputime unlimited filesize unlimited datasize 1048576 kbytes stacksize 32768 kbytes coredumpsize unlimited memoryuse 699392 kbytes descriptors 4096 files addressspace 1048576 kbytes This shows that datasize can be increased to 1048576 kbytes, or about 1 Gb. It may be possible to increase this even further by altering the operating system configuration (this may require rebuilding the kernel and is best left to an expert). iii. Run limit once again, providing the desired value for datasize, e.g.: > limit datasize 500000 which increases datasize to 500 megabytes. Any value smaller than the value returned in step ii can be used. iv. Now rerun phrap. Note that the new datasize value will only apply to the particular shell (window) in which you execute the limit command, so you may need to repeat step iii next time you log in or if you switch to a different window. 2. Other phrap- or cross_match-generated error messages (TO BE ADDED) 3. Phrap- or cross_match-generated warning messages These can generally be ignored, unless there is something obviously wrong or incomplete about the assembly (TO BE ADDED) 4. "Crashes" reported by operating system True "crashes" (e.g. segmentation fault, floating point exception) generally indicate a bug in the program. I would greatly appreciate having any such problem reported to me, as described below. 5. Long running time. If the run appears to be "stuck" or takes an exceptionally long time to complete, try using a larger value for the parameter -minmatch (above, section IV.2) 6. Misassemblies, incomplete assemblies, incorrect selection of read for consensus sequence. Marked discrepancy not split ... (TO BE ADDED) 7. Reporting problems If you have a problem that is not addressed by any of the above, first be sure you have read sections I-IV of the documentation and that you are using a current version of the programs. Then repeat the run in which the problem occurred, capturing the standard error to a file. For example, the following command, run on a Unix computer in a C shell, captures the standard output to a file phrap.out and the standard error to a file stderr.file: (phrap reads.screen > phrap.out) >& stderr.file Then send a copy of the standard error (only! -- not the standard output), i.e. the file stderr.file in the above example, to me at the following email address: email@example.com (If the standard error is more than 1 or 2 pages you have probably inadvertently captured the standard output.) Please do NOT send the standard output, datasets, or any other large file without making prior arrangements to do so with me! Large files cannot be sent to the above address, which has a limited disk quota on a Univ. of Washington computer. APPENDIX: ALGORITHMS [N.B. MUCH OF THE FOLLOWING DESCRIPTION IS SOMEWHAT OUT-OF-DATE]. Outline of phrap assembly: 0) Read in sequence & quality data, trim off any near-homopolymer runs at ends of reads, 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 and mark so they aren't used in assembly. 3) Find near duplicate reads. 4) Find reads with self-matches. 5) Find matching read pairs that are "node-rejected" i.e. do not have "solid" matching 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 alignment for each matching pair of reads that have more than one significant alignment in a given region (highest LLR-scores among several overlapping). 9) Identify probable chimeric and deletion reads (the latter are withheld from assembly). 10) Construct contig layouts, using consistent pairwise matches in decreasing score order (greedy algorithm). Consistency of layout is checked at pairwise comparison level. 11) Construct contig sequence as a mosaic of the 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. 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). 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.