|
Home > |
|
BLASTCLUSTBLAST score-based single-linkage clustering1. Clustering procedure.BLASTCLUST automatically and systematically clusters protein or DNA sequences based on pairwise matches found using the BLAST algorithm in case of proteins or Mega BLAST algorithm for DNA. In the latter case a single Mega BLAST search is performed for all the sequences combined against a database created from the same sequences. BLASTCLUST finds pairs of sequences that have statistically significant matches and clusters them using single-linkage clustering. BLASTCLUST uses the default values for the BLAST and Mega BLAST parameters. For protein sequences these are: matrix BLOSUM62; gap opening cost 11; gap extension cost 1; no low-complexity filtering. For DNA sequences: match reward 1, mismatch penalty -3, non-affine gapping costs (see README.mbl document for explanation), wordsize 28. In both cases e-value threshold is set to 1e-6. ,p>For each pair of sequences the top-scoring alignment is evaluated according to the following criteria:
x1 x2 HSP length on seqX: Hx = x2-x1+1
| | gaps in seqX: Gx
seqX ---======================----- seqX length: Lx
\\|||||||||||||||||// BLAST score: S
seqY ----====================------ number of identical residues: N
| | seqY length: Ly
y1 y2 gaps in seqY: Gy
HSP length on seqY: Hy = y2-y1+1
coverage of seqX: Cx = Hx/Lx
coverage of seqY: Cy = Hy/Ly
coverage: max(Cx,Cy) or min(Cx,Cy),
depending on the value of -b option
alignment length Al = Hx+Gx = Hy+Gy
score density: S/min(Hx,Hy) or N/Al*100%
If the coverage is above a certain threshold AND the score density is above a certain threshold, these two sequences are considered to be neighbored. Thus determined neighbor relationships is considered symmetric and provides the base for clustering by a single-linkage method (which puts a sequence to a cluster if the sequence is a neighbor to at least one sequence in the cluster). 2. Input formats.The primary input format for BLASTCLUST is a FASTA-format sequence file. Each sequence should have a unique identifier (as defined by formatdb). BLASTCLUST formats this sequence set into a BLASTable database (in the directory pointed to by the environment variable TMPDIR or in the current directory), then removes the database. Instead of a FASTA file, a database prepared by formatdb with -o option set to TRUE can be supplied as an input. Another type of input is a sequence hit-list previously saved by BLASTCLUST (in this case BLASTCLUST will use pre-computed HSP data instead of making de novo comparisons). You can restrict clustering to a subset of your data by supplying an ID list file (IDs separated by spaces, tabs, newlines, commas or semicolons). This is supposed to be used for re-clustering subsets of sequences using the previously computed hit-list file. 3. Output format.BLASTCLUST prints out clusters of sequence IDs, sorted from largest to smallest cluster (alphabetically by ID of the first sequence if of the same size), separating clusters by a newline character. Sequence identifiers within a cluster are space-separated and sorted from longest to shortest sequence (alphabetically by IDs if of the same length). 4. Crash recovery.If the program crashed because of system error you can restart it using crash recovery mode. This works only if you were saving hit-list during the clustering. Start the job with the same command line as before, specifying the hit-list saving to the same file but also set the "continue unfinished clustering" option to TRUE. The process will restart from the last saved point and will append the hit-list file. 5. Environment.BLASTCLUST is supposed to work in a normal NCBI environment, in particular: BLOSUM62 matrix is available via .ncbirc or BLASTMAT environment variable. 6. Program optionsInput: -i <file> sequence file in the FASTA format (default = stdin) -d <file> sequence database name -r <file> name of a hit-list file saved by BLASTCLUST These three options are mutually exclusive.
-l <file> a file with a list of IDs to restrict the clustering,
applicable only when reclustering from a saved hit-list.
Thresholds:
-S <threshold> similarity threshold
if <3 then the threshold is set as a BLAST score density
(0.0 to 3.0; default = 1.75)
if >=3 then the threshold is set as a percent of identical
residues (3 to 100)
-L <threshold> minimum length coverage (0.0 to 1.0; default = 0.9)
-b <T|F> require coverage as specified by -L and -S on both (T) or
only one (F) sequence of a pair (default = TRUE)
Output:
-o <file> file to save cluster list (default = stdout)
-s <file> file to save hit-list (this file may be not portable
across platforms)
-p <T|F> protein (T) or nucleotide (F) sequences in the input
(default = TRUE)
Misc:
-C <T|F> continue unfinished clustering (crash recovery mode).
(default = FALSE)
-a <number> Number of CPU's to use in a multi-thread mode
(default = 1).
-v <logfile> Progress report destination (printed every 1000
sequences).
Set to F to suppress report messages (default = stderr).
-e <T|F> Enable sequence id parsing in database formatting.
Set to F if multiple sequences have identical ids
(default = TRUE).
-W Word size to use for initial matches (default = 0, translates
to 3 for proteins and 32 for nucleotides).
-c <config file> Configuration file with advanced options,
containing any of the following options with their values,
separated by whitespace:
-r, -q, -G, -E - match, mismatch, gap open and gap extension
scores respectively,
-e - e-value cut off,
-y, -X - the dropoff values for the ungapped and gapped extension
respectively,
-A - window size for two-hit version,
-I - hitlist size,
-Y, -z - effective search space and database length respectively,
to be used for e-value and bit score calculations,
-F - filter string,
-s - raw score cut off for nucleotide search,
-S - strand option.
7. Credits and complaints:Ilya Dondoshansky (dondosha@ncbi.nlm.nih.gov) APPENDIX A.Format of the hit-list file. The hit-list file consists of the following parts: A.1. Header.
4-byte integer IDtype 1 if numeric IDs; 0 if string IDs
4-byte integer ListSz size of the ID list; if IDs are numeric this
is the number of SeqID records, otherwise this
is the length of the ID list (in bytes)
A.2. Sequence ID list.
If IDtype is 1 (numeric IDs) then the list is ListSz records of
4-byte integer SeqID sequence ID (numeric)
If IDtype is 0 (string IDs) then the list is a list of records of
var-length char SeqID sequence ID (string)
space (' ') separator
(total length is ListSz bytes; the number of sequences is equal to the number of spaces). A.3. Sequence length list. This is a list of
4-byte integer SeqLen sequence length
A.4. Hit list. The list consists of the following records going to the end of file: 4-byte integer N1 ordinal number of the 1st sequence 4-byte integer N2 ordinal number of the 2nd sequence 4-byte integer HSPL1 HSP length on the 1st sequence 4-byte integer HSPL2 HSP length on the 2nd sequence 8-byte float Score BLAST score 8-byte float PercId Percent of identical residues |
