Muscle 3.7 Create Multiple Alignments from Sequences or Profiles Robert C. Edgar Edgar, Robert C. (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput, Nucleic Acids Research 32(5), 1792-97. Edgar, R.C. (2004) MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5(1): 113. Phylogeny / Alignment muscle muscle perl "muscle_expanse" 0 scheduler_input scheduler.conf perl !$more_memory perl "ChargeFactor=1.0\\n" . "nodes=1\\n" . "mem=4G\\n" . "node_exclusive=0\\n" . "threads_per_process=1\\n" scheduler_input2 scheduler.conf perl $more_memory perl "ChargeFactor=1.0\\n" . "nodes=1\\n" . "mem=16G\\n" . "node_exclusive=0\\n" . "threads_per_process=1\\n" runtime 1 scheduler.conf Maximum Hours to Run (click here for help setting this correctly) 1.0 Estimate the maximum time your job will need to run (up to 166 hrs). Your job will be killed if it doesn't finish within the time you specify, however jobs with shorter maximum run times are often scheduled sooner than longer jobs. Maximum Hours to Run must be between 0.1 - 166 perl $runtime < 0.1 || $runtime > 166 perl "runhours=$value\\n" The job will run on 2 processors as configured; it will consume 2 X $runtime cpu hours perl !$more_memory The job will run on 8 processors as configured; it will consume 8 X $runtime cpu hours perl $more_memory infile Sequences File (must be in fasta format) (-infile) perl " -in infile.fasta" 2 infile.fasta data_type Data type (-seqtype) perl "-seqtype $value" protein dna rna auto auto 3 more_memory I need more memory 3 make_scorefile Write a Matrix Score File (-scorefile) perl ($value)? "-scorefile matrixscorefile.txt" : "" 0 5 Write a score file. This contains one line for each column in the alignment. Each line contains the letters in the column followed by the average BLOSUM62 score over pairs of letters in the column. matrixscorefile matrixscorefile.txt perl $make_scorefile clusteronly_para Cluster Only Parameters (-cluster) clusteronly Use clustering only (upgma) perl ($value)? "-cluster upgma":"" 0 1 The fast clustering algorithm in Muscle can be used to prepare data for other applications. Typical usage is: muscle -cluster -tree1 tree.phy -maxiters 1 The sequences will be clustered, and a tree written to an tree.phy. Options -weight1, -distance1, -cluster1 and -root1 can be applied if desired. UPGMA clustering is used. clusteronlymaxiters Number of iterations for Clustering (-maxiters) perl $clusteronly perl " -maxiters $value" 1 6 At least one iteration is required perl $clusteronlymaxiters < 1 runparams Alignment Parameters maxiterations Number of iterations (-maxiters) perl !$clusteronly perl " -maxiters $value" 16 8 At least one iteration must be specified perl $maxiterations < 1 If you specify 1, 2 or 3, then this is exactly the number of iterations that will be performed. If the value is greater than 3, then Muscle will continue up to the maximum you specify or until convergence is reached, which ever happens sooner. The default is 16. If you have a large number of sequences, refinement may be rather slow. If you have a large number of sequences (a few thousand), or they are very long, then the default settings of may be too slow for practical use. A good compromise between speed and accuracy is to run just the first two iterations of the algorithm. On average, this gives accuracy comparable to T-Coffee and speeds much faster than CLUSTALW. This is done by setting maxiters = 2. The default settings may be too slow for practical use if you have a few thousand sequences or more, or if they are very long. In this situation, use the option maxiters 2, running just the first two iterations. This offers a good compromise between speed and accuracy (accuracy with these settings is comparable to T-Coffee, with speeds much faster than ClustalW). max_hours Maximum time (hours) to run the program (-maxhours) perl !$clusteronly perl (defined $value) ? "-maxhours $value":"" Sorry, the value for maxhours cannot be greater than 71.5 perl $max_hours > 71.5 Sorry, the value for maxhours must be 15 minutes less than the value for runhours perl $max_hours > ($runtime - 0.25) 9 If you have a large alignment, muscle may take a long time to complete. It is sometimes convenient to say "I want the best alignment I can get in 24 hours" rather than specifying a set of options that will take an unknown length of time. This is done by using –maxhours, which specifies a floating-point number of hours. If this time is exceeded, muscle will write out current alignment and stop. Note that the actual time may exceed the specified limit by a few minutes while muscle finishes up on a step. It is also possible for no alignment to be produced if the time limit is too small max_mem Set the Maximum Memory to be allocated for this job (maxmb) perl "-maxmb $value" 30000000 10 Specifies the amount of memory (in mb) that MUSCLE uses. Allowed values are 500-2000 mb. If the limit is exceeded, MUSCLE quits, saving the best alignment so far produced (if any). This feature has been hacked on top of code that wasn't really designed for it. So it doesn't always work perfectly, but is better than nothing. The ideal solution would be to implement linear space dynamic programming code (e.g., the Myers-Miller algorithm) for situations where memory is tight. protprofilefunction Specify Protein Profile Scoring Function perl $data_type eq "protein" perl "$value" -le -sp -sv 11 Three different protein profile scoring functions are supported: the log-expectation score (–le option) and 2 sum of pairs scores: either the PAM200 (–sp) or the VTML240 matrix (–sv). The log-expectation score is the default as it gives better results on our tests, but is typically somewhere between two or three times slower than the sum-of-pairs score. For nucleotides, –spn is currently the only option (which is of course the default for nucleotide data, so you don't need to specify this option) hydrowindow Window size for determining whether a region is hydrophobic. Protein only. (-hydro) perl !$clusteronly && $data_type eq "protein" perl "-hydro $value" 5 12 The hydrophobic window size must be must be greater than or equal to 1 perl $hydrowindow < 1 hydrofactormultiplier Multiplier for gap open/close penalties in hydrophobic regions, protein only (-hydrofactor) perl !$clusteronly && $data_type eq "protein" perl " -hydrofactor $value" 1.2 13 The hydrophobic multiplier must be greater than 0 perl $hydrofactormultiplier <= 0 logfile Produce a logfile containing the command line? (-log) perl ($value) ? " -log logfile.txt":"" 1 15 You can specify a log file by using this option. A message will be written to the log file when muscle starts and stops. Error and warning messages will also be written to the log. verboselogfile Produce a more verbose logfile? (-verbose) perl ($value) ? " -verbose":"" 0 16 If –verbose is specified, then more information will be written, including the command line used to invoke muscle, the resulting internal parameter settings, and also progress messages.IThe author highly recommends using –verbose and –log;. This enables a user to verify whether a particular alignment was completed and to review any errors or warnings that occurred. logfileout perl $logfile logfile.txt verboselogfileout perl $verboselogfile vlogfile.txt User_matrix User Supplied Matrices add_matrix User Supplied Matrix (-matrix) perl !$clusteronly perl ($value)? "-matrix" : "" 0 18 You can specify your own substitution matrix by using the -matrix option. This reads a protein substitution matrix in NCBI or WU-BLAST format. The alphabet is assumed to be amino acid, and sum-of-pairs scoring is used. The -gapopen, -gapextend and -center parameters should be specified; normally you will specify a zero value for the center. Note that gap penalties MUST be negative. The environment variable MUSCLE_MXPATH can be used to specify a path where the matrices are stored. For example: muscle -in seqs.fa -out seqs.afa -matrix blosum62 -gapopen -12.0 -gapextend -1.0 -center 0.0 You can hack a nucleotide matrix by pretending that AGCT are amino acids and making a 20x20 matrix out of the original 4x4 matrix. mymatrix Matrix File (NCBI or WU-BLAST format) perl $add_matrix perl "" 19 user_matrix gpopen Gap Open Penalty (must be negative) perl $add_matrix perl (defined $value) ? " -gapopen $value":"" The gap open penalty must be a negative number perl $gpopen >= 0 20 gpextend Gap Extend Penalty (must be negative) perl $add_matrix perl (defined $value) ? " -gapextend $value":"" The gap extend penalty must be a negative number perl $gpextend >= 0 21 ctr Center Value (must be zero or negative) perl $add_matrix perl (defined $value) " -center $value":"" 0 The center value is usually 0, MUST be less than or equal to zero perl $ctr > 0 22 Diagonal_Options Diagonal Options use_diags Use Diagonals perl !$clusteronly perl ($value) ? "" : "" 0 23 diagfunction Diagonal Functions (-diags) perl $use_diags perl "$value" -diags -diags1 -diags2 -diags 24 The –diags option optimizes for speed by finding common words (6-mers in a compressed amino acid alphabet) between the two sequences as seeds for diagonals. For large numbers of closely related sequences, this option works very well, but sometimes has lower average accuracy. If you want the fastest possible speed, then the following example shows the applicable options for proteins: Maxiters=1; Diagonal Iterations 1 and 2; Profile Matrix=VTML240; Distance measure =kbit20_3 If the number of sequences exceeds, say, several thousand, or if they are very long, the kbit20_3 option may cause a memory issue. If so, you can try using the default distance settings. Maxiters=1; Diagonal Iteration 1; Profile Matrix=VTML240 For nucleotides, use: Maxiters=1; Diagonal Iterations 1 and 2 Muscle uses k-mer extension to find diagonals. It is disabled by default, as it causes a slight reduction in average accuracy. You may select k-mer extension separately for iteratons 1 or 2 (or both) because it would be reasonable to enable diagonals in the first iteration but not the second. The main goal of the first iteration is to construct a multiple alignment quickly in order to improve the distance matrix, which is not very sensitive to alignment quality. The goal of the second iteration is to make the best possible progressive alignment maxdiagdist Maximum distance between two diagonals that allows them to merge into one diagonal. (-diagbreak) perl $use_diags perl " -diagbreak $value" 1 The maximum distance between two diagonals must be greater than or equal to 1 perl $maxdiagdist < 1 25 mindiaglength Minimum length of diagonal (-diaglength) perl $use_diags perl " -diaglength $value" 24 26 The minimum length of the diagonal must not be less than 1 perl $mindiaglength < 1 diagmarginpos Discard this many positions at ends of diagonal (-diagmargin) perl $use_diags perl " -diagmargin $value" 5 27 You must discard at least 1 position at the end of a diagonal perl $diagmarginpos < 1 tree_options Tree parameters use_guide_tree Provide a Guide Tree perl !$clusteronly perl ($value)? "-usetree_nowarn" : "" 28 0 Users can provide a guide tree to help with their sequence alignments. guidetree Guide Tree File; Newick/Phylip format (-usetree) perl $use_guide_tree perl "guidetree.dnd" 29 guidetree.dnd You can specify a guide tree for your alignment. It must be in Phylip/Newick format (.dnd file). WARNING: Do not use this option just because you believe that you have an accurate evolutionary tree for your sequences. The best guide tree for multiple sequence alignment is not the correct evolutionary tree. This can be understood by the following argument. Alignment accuracy decreases with lower sequence identity. It follows that a given set of profiles the two that can be aligned most accurately will tend to be the pair with the highest identity, i.e. at the shortest evolutionary distance. This is exactly the pair selected by the nearest-neighbor criterion that Muscle uses by default. When mutation rates are variable, the evolutionary neighbor may not be the nearest neighbor. This explains why a nearest-neighbor tree may be superior to the true evolutionary tree for for guiding a progressive alignment. make_tree1file Save tree produced in first iteration in Newick (Phylip-compatible) format perl ($value)? "-tree1 treeiteration1.dnd" : "" 30 maketree1filename Save tree produced in second iteration to given file in Newick (Phylip-compatible) format perl $make_tree1file treeiteration1.dnd make_tree2file Save tree produced in second iteration in Newick (Phylip-compatible) format perl !$clusteronly perl ($value)? "-tree2 treeiteration2.dnd" : "" 31 maketree2filename Save tree produced in second iteration to given file in Newick (Phylip-compatible) format perl $make_tree2file treeiteration2.dnd options1 First Stage Options use_weights_1 Use Weights for Iterations 1 and 2 perl ($value)? "" : "" 30 0 weight1 Sequence weighting scheme for Iteration 1: (-weight1) perl " -weight1 $value" clustalw henikoff henikoffpb gsc clustalw threeway 32 weight1 is used in iterations 1 and 2; none=all sequences have equal weight; Henokoff = Henikoff and Henikoff, Modified Henikoff = Henikoff as used in Psi-BLAST[ ClustalW=CLUSTALW method, Gotoh = Gotoh three-way method. distance1p Distance measure for iteration 1; Protein only (-distance1) perl $data_type eq "protein" perl "-distance1 $value" kmer6_6 kmer6_6 kmer20_3 kmer20_4 kbit20_3 kmer4_6 33 cluster1 Clustering for Iterations 1 and 2: (-cluster1) perl " -cluster1 $value " upgmb upgma upgmb neighborjoining 35 UPGMB is a method based on a combination of both UPGMA and Neighbor Joining (called here UPGMB). sueff1 SUEFF(must be between 0 and 1) perl $cluster1 eq "upgmb" perl " -sueff $value" 0.1 The SUEFF Value Must be > 0 and 1 or less perl $sueff1 <= 0 The SUEFF Value Must be > 0 and 1 or less perl $sueff1 > 1 36 root1 Method used to root tree in iterations 1 and 2 (-root1) perl "-root1 $value" pseudo pseudo midlongestspan minavgleafdist 37 options2 Later Stage Options maxtrees Maximum trees in Iteration 2 (-maxtrees) perl !$clusteronly perl "-maxtrees $value" 1 38 The value must be greater than or equal to 1 perl $maxtrees < 1 This option controls the maximum number of new trees to create in iteration 2. A point of diminishing returns is typically reached after the first tree. If a larger value is given, the process will repeat until convergence or until this number of trees has been created, whichever comes first. use_weights_2 Use Weights for Iterations 3 and beyond perl ($value)? "" : "" 30 0 weight2 Sequence weighting for tree-dependent refinement. perl !$clusteronly perl " -weight2 $value" clustalw henikoff henikoffpb gsc clustalw threeway 39 weight2 is used in iterations 1 and 2; none=all sequences have equal weight; Henokoff = Henikoff and Henikoff, Modified Henikoff = Henikoff as used in Psi-BLAST[ ClustalW=CLUSTALW method, Gotoh = Gotoh three-way method. distance2 Distance Measure for Iterations 2 and above (-distance2) perl !$clusteronly && $data_type eq "protein" perl " -distance2 $value" kmer20_3 kmer20_4 kbit20_3 pctidkimura pctidlog pctidkimura 40 cluster2 Clustering method for later iterations: (-cluster2) perl !$clusteronly perl " -cluster2 $value" upgmb upgma upgmb neighborjoining 41 sueff2 SUEFF(must be between 0 and 1) perl $cluster2 eq "upgmb" perl " -sueff $value" 0.1 The SUEFF Value Must be > 0 and 1 or less perl $sueff2 <= 0 The SUEFF Value Must be > 0 and 1 or less perl $sueff2 > 1 42 root2 Method used to root tree in iterations 3 and beyond (-root2) perl !$clusteronly perl "-root2 $value" pseudo pseudo midlongestspan minavgleafdist 43 objectivescore Objective score used by tree dependent refinement (-objscore) perl !$clusteronly perl "-objscore $value" sp sp spf spm ps dp xp 44 anchor Use Anchor Optimization (-noanchors) This option has no effect if -maxiters 1 or -maxiters 2 is specified. perl ($value)? "":"-noanchors" 0 45 Tree-dependent refinement (iterations 3, 4 ... ) can be speeded up by dividing the alignment vertically into blocks. Block boundaries are found by identifying high-scoring columns (e.g., a perfectly conserved column of Cs or Ws would be a candidate). Each vertical block is then refined independently before reassembling the complete alignment, which is faster because of the L2 factor in dynamic programming (e.g., suppose the alignment is split into two vertical blocks, then 2 x (0.5)e2 = 0.5, so the dynamic programming time is roughly halved). The –noanchors option is used to disable this feature. This option has no effect if –maxiters 1 or –maxiters 2 is specified. On benchmark tests, enabling anchors has little or no effect on accuracy, but if you want to be very conservative and are striving for the best possible accuracy then –noanchors is a reasonable choice. minanchorcolumnspacing Minimum spacing between anchor columns perl $anchor 32 perl " -anchorspacing $value" 46 The minimum spacing between anchor columns must not be less than 1 perl $minanchorcolumnspacing <1 minanchorcolumnscore Minimum score a column must have to be an anchor perl $anchor 1 perl " -anchorspacing $value" 47 The minimum score a column to be an anchor cannot be less than 1 perl $minanchorcolumnscore <1 minsmoothanchorscore Minimum smoothed score a column must have to be an anchor perl $anchor 1 48 The minimum smoothed score to be an anchor cannot be less than 1 perl $minsmoothanchorscore <1 anchorsmoothingwindow Window Used for Anchor Smoothing (-smoothwindow; default=7) perl $anchor 7 perl " -smoothwindow $value" 49 The window smoothing value must be at least 1 perl $anchorsmoothingwindow < 1 output_options Output Options (select as many as you like) out_default Default submit string perl !$multi_output perl "-out output.fasta" 50 multi_output Specify output formats 0 If you specify output format AND set -maxhours, your output might appear in the stdout.txt file perl defined $max_hours 51 output_fasta Fasta output perl $multi_output perl ($value) ? " -fastaout output.fasta":"" 0 51 output_clustal Clustal output (.aln) perl $multi_output perl ($value) ? " -clwout output.aln":"" 0 52 output_clustals Strict Clustal output (.aln; has the Clustal header) perl $multi_output perl ($value) ? " -clwstrictout outputs.aln":"" 0 53 output_html HTML output (.html) perl $multi_output perl ($value) ? " -htmlout output.html":"" 0 54 output_gcg GCG output (.msf) perl $multi_output perl ($value) ? "-msfout output.msf":"" 0 55 output_phyi Interleaved Phylip output perl $multi_output perl ($value) ? " -phyiout outputi.phy":"" 0 56 output_phys Sequential Phylip output (.phy) perl $multi_output perl ($value) ? " -physout outputs.phy":"" 0 57 fasta_output *.fasta clustal_output perl $output_clustal output.aln clustal_output_strict perl $output_clustals outputs.aln html_output perl $output_html *.html gcg_output perl $output_gcg *.msf phylip_output_interleaved perl $output_phyi outputi.phy phylip_ouput_sequential perl $output_phyi outputs.phy all_txt *.txt jobinfo *.TXT