* Each DNA sequence consists of a sequence of sites. Each site has a * state, which is a set of bases. The four bases are adenine, * cytosine, guanine, and thymine. For textual I/O, each state is represented by * a single character as follows: *
*
| Char. | * | Meaning | * | Set |
| A | * | Adenine | * | (A) |
| C | * | Cytosine | * | (C) |
| G | * | Guanine | * | (G) |
| T | * | Thymine | * | (T) |
| Y | * | pYrimidine | * | (C or T) |
| R | * | puRine | * | (A or G) |
| W | * | "Weak" | * | (A or T) |
| S | * | "Strong" | * | (C or G) |
| K | * | "Keto" | * | (G or T) |
| M | * | "aMino" | * | (A or C) |
| B | * | not A | * | (C or G or T) |
| D | * | not C | * | (A or G or T) |
| H | * | not G | * | (A or C or T) |
| V | * | not T | * | (A or C or G) |
| X | * | unknown | * | (A or C or G or T) |
| - | * | deletion | * | () |
* The DNA sequence file format is that used by Joseph Felsenstein's Phylogeny * Inference Package (PHYLIP). While the file is a plain text file, it often has * the extension ".phy" to indicate that it is in PHYLIP format. For * further information, see: *
* Here is an example of an input file: *
*
* 5 42 * Turkey AAGCTNGGGC ATTTCAGGGT * Salmo gair AAGCCTTGGC AGTGCAGGGT * H. Sapiens ACCGGTTGGC CGTTCAGGGT * Chimp AAACCCTTGC CGTTACGCTT * Gorilla AAACCCTTGC CGGTACGCTT * * GAGCCCGGGC AATACAGGGT AT * GAGCCGTGGC CGGGCACGGT AT * ACAGGTTGGC CGTTCAGGGT AA * AAACCGAGGC CGGGACACTC AT * AAACCATTGC CGGTACGCTT AA* |
*
* The first line contains the number of species S and the number of * sites N in each sequence. S must be >= 2. N must be * >= 1. *
* The next S lines contain the initial data for each species. The first * ten characters contain the sequence name. This must be exactly ten * characters, padded with blanks if necessary. Then comes one character for * each site in the sequence. Uppercase and lowercase are considered the same. * Characters other than those for the states listed above are ignored. Often, a * blank is inserted every ten characters for readability, but this is not * necessary. After these S lines come zero or more blank lines for * readability, which are ignored. If there is more sequence data, the next * S lines give the states for the next sites in the sequences. This * continues for the rest of the file. *
* This is known as the "interleaved" file format. There is also a "sequential" * file format, but the sequential file format is not supported. *
* Thus, the complete sequence for each species in the example is: *
*
| Species | * | Sequence |
| Turkey | * | AAGCTNGGGCATTTCAGGGTGAGCCCGGGCAATACAGGGTAT |
| Salmo gair | * | AAGCCTTGGCAGTGCAGGGTGAGCCGTGGCCGGGCACGGTAT |
| H. Sapiens | * | ACCGGTTGGCCGTTCAGGGTACAGGTTGGCCGTTCAGGGTAA |
| Chimp | * | AAACCCTTGCCGTTACGCTTAAACCGAGGCCGGGACACTCAT |
| Gorilla | * | AAACCCTTGCCGGTACGCTTAAACCATTGCCGGTACGCTTAA |
* In the input file, the following alternate characters can be used: X, N, and * ? all mean "unknown." O (capital letter O) and - (hyphen) both mean * "deletion." The character . (period) means "the same as the corresponding * site in the first species." Here is another input file with the same * sequences as the one above: *
*
* 5 42 * Turkey AAGCTNGGGC ATTTCAGGGT * Salmo gair ..G.CTT... AG.G...... * H. Sapiens .CCGGTT... .G........ * Chimp ..A.CCTT.. .G..AC.CT. * Gorilla ..A.CCTT.. .GG.AC.CT. * * GAGCCCGGGC AATACAGGGT AT * .....GT... CGGG..C... .. * ACAGGTT... CG.T...... .A * A.A..GA... CGGGACACTC .. * A.A..ATT.. CGGTAC.CT. .A* |
*
* Here are some more example DNA sequence files: *
* Note: The DNA sequences in the new list are copies of (not * references to) the DNA sequences in the given list. * * @param list DNA sequence list to copy. * * @exception NullPointerException * (unchecked exception) Thrown if list is null. */ public DnaSequenceList (DnaSequenceList list) { int N = list.mySequence.length; this.mySequence = new DnaSequence [N]; for (int i = 0; i < N; ++ i) { this.mySequence[i] = new DnaSequence (list.mySequence[i]); } if (list.isInformative != null) { this.isInformative = (boolean[]) list.isInformative.clone(); } this.nInformative = list.nInformative; this.nChanges = list.nChanges; } // Exported operations. /** * Obtain this DNA sequence list's length. * * @return Length N (number of DNA sequences). */ public int length() { return mySequence.length; } /** * Get the DNA sequence at the given index in this DNA sequence list. * * @param i Index, 0 ≤ i ≤ N−1. * * @return DNA sequence. * * @exception ArrayIndexOutOfBoundsException * (unchecked exception) Thrown if i is out of bounds. */ public DnaSequence seq (int i) { return mySequence[i]; } /** * Read a DNA sequence list from the given input file. The input file must * be in interleaved PHYLIP format. *
* The DNA sequences' sites and names are read from the input file. The DNA
* sequences' scores are set to 0.
*
* @param file File.
*
* @return DNA sequence list.
*
* @exception NullPointerException
* (unchecked exception) Thrown if file is null.
* @exception IOException
* Thrown if an I/O error occurred. Thrown if the input file's contents
* were invalid.
*/
public static DnaSequenceList read
(File file)
throws IOException
{
Scanner filescanner = new Scanner (file);
Scanner linescanner;
int S, N;
DnaSequenceList list;
int[] sitecount;
String line;
try
{
// Read number of species and number of sites from first line.
if (! filescanner.hasNextLine())
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Empty file");
}
linescanner = new Scanner (filescanner.nextLine());
if (! linescanner.hasNextInt())
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Number of species invalid or missing");
}
S = linescanner.nextInt();
if (S < 2)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Number of species must be >= 2");
}
if (! linescanner.hasNextInt())
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Number of sites invalid or missing");
}
N = linescanner.nextInt();
if (N < 1)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Number of sites must be >= 1");
}
// Set up DNA sequence list and site count array.
list = new DnaSequenceList();
list.mySequence = new DnaSequence [S];
sitecount = new int [S];
// Read sequence data from groups of S lines until EOF.
fileloop: for (;;)
{
speciesloop: for (int s = 0; s < S; ++ s)
{
// Get a line of sequence data for species s.
if (filescanner.hasNextLine())
{
}
else if (s != 0 || sitecount[s] == 0)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Missing a line of sequence data for species " +
(s+1));
}
else
{
break fileloop;
}
line = filescanner.nextLine();
// Ignore blank lines.
if (line.trim().equals (""))
{
-- s;
continue;
}
// The first time, extract sequence name and create
// DnaSequence object.
if (sitecount[s] == 0)
{
if (line.length() < 10)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Name must be 10 characters for species " +
(s+1));
}
list.mySequence[s] =
new DnaSequence
(N, 0, line.substring (0, 10) .trim());
line = line.substring (10);
}
// Parse characters in sequence data.
int len = line.length();
byte[] seq = list.mySequence[s].mySites;
byte[] seq0 = list.mySequence[0].mySites;
int count = sitecount[s];
for (int i = 0; i < len; ++ i)
{
switch (line.charAt(i))
{
case 'O': case 'o': case '-':
verifyCount (count, N, file, s);
seq[count] = (byte) 0; // ----
++ count;
break;
case 'A': case 'a':
verifyCount (count, N, file, s);
seq[count] = (byte) 1; // ---A
++ count;
break;
case 'C': case 'c':
verifyCount (count, N, file, s);
seq[count] = (byte) 2; // --C-
++ count;
break;
case 'M': case 'm':
verifyCount (count, N, file, s);
seq[count] = (byte) 3; // --CA
++ count;
break;
case 'G': case 'g':
verifyCount (count, N, file, s);
seq[count] = (byte) 4; // -G--
++ count;
break;
case 'R': case 'r':
verifyCount (count, N, file, s);
seq[count] = (byte) 5; // -G-A
++ count;
break;
case 'S': case 's':
verifyCount (count, N, file, s);
seq[count] = (byte) 6; // -GC-
++ count;
break;
case 'V': case 'v':
verifyCount (count, N, file, s);
seq[count] = (byte) 7; // -GCA
++ count;
break;
case 'T': case 't':
verifyCount (count, N, file, s);
seq[count] = (byte) 8; // T---
++ count;
break;
case 'W': case 'w':
verifyCount (count, N, file, s);
seq[count] = (byte) 9; // T--A
++ count;
break;
case 'Y': case 'y':
verifyCount (count, N, file, s);
seq[count] = (byte) 10; // T-C-
++ count;
break;
case 'H': case 'h':
verifyCount (count, N, file, s);
seq[count] = (byte) 11; // T-CA
++ count;
break;
case 'K': case 'k':
verifyCount (count, N, file, s);
seq[count] = (byte) 12; // TG--
++ count;
break;
case 'D': case 'd':
verifyCount (count, N, file, s);
seq[count] = (byte) 13; // TG-A
++ count;
break;
case 'B': case 'b':
verifyCount (count, N, file, s);
seq[count] = (byte) 14; // TGC-
++ count;
break;
case 'X': case 'x': case 'N': case 'n': case '?':
verifyCount (count, N, file, s);
seq[count] = (byte) 15; // TGCA
++ count;
break;
case '.':
verifyCount (count, N, file, s);
if (s == 0)
{
throw new IOException
("DnaSequenceList.read(\"" + file +
"\"): " +
"'.' not allowed in species 1");
}
if (count >= sitecount[0])
{
throw new IOException
("DnaSequenceList.read(\"" + file +
"\"): " +
"'.' in species " + (s+1) +
" has no corresponding site in species 1");
}
seq[count] = seq0[count];
++ count;
break;
}
}
sitecount[s] = count;
}
}
// Verify correct site count for all species.
for (int s = 0; s < S; ++ s)
{
if (sitecount[s] < N)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Too few sites for species " + (s+1));
}
else if (sitecount[s] > N)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Too many sites for species " + (s+1));
}
}
// Return DNA sequence list.
return list;
}
finally
{
filescanner.close();
}
}
private static void verifyCount
(int count,
int N,
File file,
int s)
throws IOException
{
if (count >= N)
{
throw new IOException
("DnaSequenceList.read(\"" + file + "\"): " +
"Too many sites for species " + (s+1));
}
}
/**
* Write this DNA sequence list to the given output file. The output file is
* in interleaved PHYLIP format. There are 70 sites on each output line.
* Periods are not used. Informative sites are not marked in bold.
*
* @param file File.
*
* @exception NullPointerException
* (unchecked exception) Thrown if file is null.
* @exception IOException
* Thrown if an I/O error occurred.
*/
public void write
(File file)
throws IOException
{
write (file, 70, false, false);
}
/**
* Write this DNA sequence list to the given output file. The output file is
* in interleaved PHYLIP format.
*
* @param file File.
* @param sites Number of sites per output line.
* @param periods True to use periods, false not to use periods.
* @param bold True to mark informative sites in bold, false not to.
*
* @exception NullPointerException
* (unchecked exception) Thrown if file is null.
* @exception IllegalArgumentException
* (unchecked exception) Thrown if sites <= 10.
* @exception IOException
* Thrown if an I/O error occurred.
*/
public void write
(File file,
int sites,
boolean periods,
boolean bold)
throws IOException
{
PrintStream ps =
new PrintStream
(new BufferedOutputStream
(new FileOutputStream (file)));
try
{
write (ps, sites, periods, bold);
}
finally
{
ps.close();
}
}
/**
* Write this DNA sequence list to the given print stream in interleaved
* PHYLIP format.
*
* @param ps Print stream.
* @param sites Number of sites per output line.
* @param periods True to use periods, false not to.
* @param bold True to mark informative sites in bold, false not to.
*
* @exception NullPointerException
* (unchecked exception) Thrown if ps is null.
* @exception IllegalArgumentException
* (unchecked exception) Thrown if sites <= 10.
* @exception IOException
* Thrown if an I/O error occurred.
*/
public void write
(PrintStream ps,
int sites,
boolean periods,
boolean bold)
throws IOException
{
if (sites <= 10)
{
throw new IllegalArgumentException
("DnaSequenceList.write(): sites = " + sites + " illegal");
}
// Determine informative sites if necessary.
if (bold) computeInformativeSites();
// Print number of species and number of sites.
int S = mySequence.length;
int N = mySequence[0].myLength;
ps.print (S);
ps.print (' ');
ps.print (N);
ps.println();
// Print groups of sites for each species. On the first line, print
// sequence name, padded or truncated to 10 characters.
int lb = 0;
int ub = Math.min (sites-10, N);
byte[] seq0 = mySequence[0].mySites;
while (lb < N)
{
for (int s = 0; s < S; ++ s)
{
byte[] seq = mySequence[s].mySites;
if (lb == 0) ps.print (padName (mySequence[s].myName));
for (int i = lb; i < ub; ++ i)
{
if ((lb == 0 || i > lb) && i % 10 == 0)
{
ps.print (' ');
}
if (periods && s > 0 && seq[i] == seq0[i])
{
printSite (ps, i, '.', bold);
}
else
{
printSite
(ps, i, DnaSequence.state2char[seq[i]], bold);
}
}
ps.println();
}
ps.println();
lb = ub;
ub = Math.min (ub+sites, N);
}
// Check for I/O errors.
if (ps.checkError())
{
throw new IOException ("DnaSequenceList.write(): I/O error");
}
}
private static String padName
(String name)
{
if (name == null) return "
* Each site in the DNA sequences is either "uninformative" or
* "informative," defined as follows:
*
* Since the uninformative sites do not affect the outcome of a maximum
* parsimony phylogenetic tree search, the uninformative sites can be
* omitted from the tree scoring process to save time. The informative sites
* do affect the outcome and must be included in the tree scoring process.
*
* The exciseUninformativeSites() removes the uninformative sites
* from the DNA sequences in this list. The DNA sequences' scores and names
* are unchanged.
*
* @return Number of state changes the (excised) uninformative sites
* contribute to the parsimony score.
*/
public int exciseUninformativeSites()
{
int S = mySequence.length;
int N = mySequence[0].length();
// Determine which sites are informative.
computeInformativeSites();
// Excise uninformative sites from sequences.
for (int s = 0; s < S; ++ s)
{
byte[] oldSites = mySequence[s].mySites;
mySequence[s] =
new DnaSequence
(nInformative,
mySequence[s].myScore,
mySequence[s].myName);
byte[] excSites = mySequence[s].mySites;
int j = 0;
for (int i = 0; i < N; ++ i)
{
if (isInformative[i])
{
excSites[j++] = oldSites[i];
}
}
}
// Mark all sites as informative.
isInformative = new boolean [nInformative];
Arrays.fill (isInformative, true);
// Return number of state changes.
return nChanges;
}
/**
* Returns the number of informative sites in this DNA sequence list.
*
* @return Number of informative sites.
*/
public int informativeSiteCount()
{
computeInformativeSites();
return nInformative;
}
/**
* Determine the number of absent states after adding each sequence in this
* DNA sequence list to a tree. The return value A is an
* N-element array, where N is the length of this DNA sequence
* list. As sequences from this list are added to a tree in order from
* i = 0 to N−1, A[i] is the number of
* character states that do not yet appear in the tree. Thus, the number of
* state changes in the tree must increase by at least A[i]
* when the sequences after sequence i are added to the tree. This
* can be used to prune a branch-and-bound search.
*
* @return Array A.
*/
public int[] countAbsentStates()
{
int N = mySequence.length;
int L = mySequence[0].length();
int[] A = new int [N];
// Compute the union of all the DNA sequences.
byte[] sites = new byte [L];
for (int i = 0; i < N; ++ i)
{
byte[] mysites_i = mySequence[i].mySites;
for (int j = 0; j < L; ++ j)
{
sites[j] |= mysites_i[j];
}
}
// Subtract each sequence from the union, count and record states.
for (int i = 0; i < N; ++ i)
{
byte[] mysites_i = mySequence[i].mySites;
int count = 0;
for (int j = 0; j < L; ++ j)
{
sites[j] &= ~ mysites_i[j];
count += DnaSequence.state2bitCount [sites[j]];
}
A[i] = count;
}
return A;
}
/**
* Create a DNA sequence tree from this DNA sequence list and the given tree
* signature. The tree signature is an array of indexes of length N,
* where N is the length of this list. To construct the tree, for all
* i from 0 to N−1, the DNA sequence at index i
* in this list is added to the tree at index signature[i] using
* the DnaSequenceTree.add() method. For all i,
* signature[i] must be in the range 0 ..
* 2(i − 1), except signature[0] is 0.
*
* Note: The returned tree has references to (not copies of) the DNA
* sequences in this list.
*
* @param signature Tree signature (array of tree indexes).
*
* @return Tree.
*/
public DnaSequenceTree toTree
(int[] signature)
{
int N = mySequence.length;
DnaSequenceTree tree = new DnaSequenceTree (2*N - 1);
for (int i = 0; i < N; ++ i)
{
tree.add (signature[i], mySequence[i]);
}
return tree;
}
/**
* Returns an iterator for the DNA sequences in this list.
*
* @return Iterator.
*/
public Iterator
*
*