c7caa3b6bb066265239e6695223a7a874d26c1ad
[jalview.git] / src / jalview / io / MegaFile.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8.0b1)
3  * Copyright (C) 2014 The Jalview Authors
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
10  *  
11  * Jalview is distributed in the hope that it will be useful, but 
12  * WITHOUT ANY WARRANTY; without even the implied warranty 
13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14  * PURPOSE.  See the GNU General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
17  * The Jalview Authors are detailed in the 'AUTHORS' file.
18  */
19 package jalview.io;
20
21 import jalview.datamodel.AlignmentI;
22 import jalview.datamodel.Sequence;
23 import jalview.datamodel.SequenceFeature;
24 import jalview.datamodel.SequenceI;
25
26 import java.io.IOException;
27 import java.util.ArrayList;
28 import java.util.HashMap;
29 import java.util.Iterator;
30 import java.util.LinkedHashMap;
31 import java.util.List;
32 import java.util.Map;
33 import java.util.Map.Entry;
34 import java.util.Set;
35
36 /**
37  * A parser for input or output of MEGA format files. <br>
38  * <br>
39  * Tamura K, Stecher G, Peterson D, Filipski A, and Kumar S (2013) MEGA6:
40  * Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology and
41  * Evolution 30: 2725-2729. <br>
42  * <br>
43  * 
44  * MEGA file format is supported as described in
45  * http://www.megasoftware.net/manual.pdf <br>
46  * Limitations:
47  * <ul>
48  * <li>nested comments (marked by [ ]) are accepted but not preserved</li>
49  * <li>to be completed</li>
50  * </ul>
51  * 
52  * @see http://www.megasoftware.net/
53  */
54 public class MegaFile extends AlignFile
55 {
56   private static final String WHITESPACE = "\\s+";
57
58   private static final int DEFAULT_LINE_LENGTH = 60;
59
60   private static final String INDENT = "    ";
61
62   private static final String N_SITES = "NSites";
63
64   private static final String N_SEQS = "NSeqs";
65
66   private static final String MISSING = "Missing";
67
68   private static final String IDENTICAL = "Identical";
69
70   private static final String INDEL = "Indel";
71
72   private static final String CODETABLE = "CodeTable";
73
74   private static final String PROTEIN = "Protein";
75
76   private static final String NUCLEOTIDE = "Nucleotide";
77
78   private static final String DATATYPE = "DataType";
79
80   private static final char COMMENT_START = '[';
81
82   private static final char COMMENT_END = ']';
83
84   private static final String HASHSIGN = "#";
85
86   private static final String SEMICOLON = ";";
87
88   private static final String BANG = "!";
89
90   private static final String EQUALS = "=";
91
92   private static final String MEGA_ID = HASHSIGN + "MEGA";
93
94   private static final String TITLE = "Title";
95
96   private static final String FORMAT = "Format";
97
98   private static final String DESCRIPTION = "Description";
99
100   private static final String GENE = "Gene";
101
102   private static final String DOMAIN = "Domain";
103
104   private static final String PROPERTY = "Property";
105
106   private static final String CODONSTART = "CodonStart";
107
108   /*
109    * names of properties to save to the alignment (may affect eventual output
110    * format)
111    */
112   static final String PROP_TITLE = "MEGA_TITLE";
113
114   static final String PROP_INTERLEAVED = "MEGA_INTERLEAVED";
115
116   static final String PROP_DESCRIPTION = "MEGA_DESCRIPTION";
117
118   static final String PROP_CODETABLE = "MEGA_CODETABLE";
119
120   static final String PROP_IDENTITY = "MEGA_IDENTITY";
121
122   static final String PROP_MISSING = "MEGA_MISSING";
123
124   static final String PROP_DATATYPE = "MEGA_DATATYPE";
125
126   // number of bases per line of file (value is inferred)
127   static final String PROP_LINELENGTH = "MEGA_LINELENGTH";
128
129   // TODO: need a controlled name for Gene as a feature if we want to be able to
130   // output the MEGA file with !Gene headers
131   // WTF do we do if the sequences get realigned?
132
133   // initial size for sequence data buffer
134   private static final int SEQBUFFERSIZE = 256;
135
136   private static final String SPACE = " ";
137
138   /*
139    * number of sequence positions output per line
140    */
141   private int positionsPerLine;
142
143   private String title;
144
145   // gap character may be explicitly declared, default is -
146   private char gapCharacter = '-';
147
148   // identity character if declared
149   private char identityCharacter = 0;
150
151   // this can be True, False or null (meaning not asserted in file)
152   private Boolean nucleotide;
153
154   // set once we have seen one block of interleaved data
155   private boolean firstDataBlockRead = false;
156
157   // this can be True, False or null (meaning we don't know yet)
158   private Boolean interleaved;
159
160   // write end of line positions as a comment
161   private boolean writePositionNumbers = true;
162
163   // id of sequence being processed
164   private String currentSequenceId;
165
166   /*
167    * Temporary store of {sequenceId, positionData} while parsing interleaved
168    * sequences; sequences are maintained in the order in which they are added
169    * i.e. read in the file
170    */
171   Map<String, StringBuilder> seqData;
172   
173   // number of residues read (so far) per sequence
174   Map<String, Integer> residuesRead;
175   
176   // current Gene if any we are parsing
177   private String currentGene;
178
179   // start residue (base 1) per sequence of current gene
180   Map<String, Integer> geneStart;
181
182   // current Domain if any we are parsing
183   private String currentDomain;
184
185   // start residue (base 1) per sequence of current domain
186   Map<String, Integer> domainStart;
187
188   // map of SequenceFeature's by sequence id
189   Map<String, List<SequenceFeature>> sequenceFeatures;
190
191   public MegaFile()
192   {
193   }
194
195   public MegaFile(String inFile, String type) throws IOException
196   {
197     super(inFile, type);
198   }
199
200   public MegaFile(FileParse source) throws IOException
201   {
202     super(source);
203   }
204
205   /**
206    * Parse the input stream.
207    */
208   @Override
209   public void parse() throws IOException
210   {
211     gapCharacter = '-';
212     sequenceFeatures = new HashMap<String, List<SequenceFeature>>();
213     geneStart = new HashMap<String, Integer>();
214     domainStart = new HashMap<String, Integer>();
215     residuesRead = new HashMap<String, Integer>();
216
217     /*
218      * Read and process MEGA and Title/Format/Description headers if present.
219      * Returns the first data line following the headers.
220      */
221     String dataLine = parseHeaderLines();
222
223     /*
224      * order-preserving map to hold sequences by id as they are built up during
225      * parsing
226      */
227     seqData = new LinkedHashMap<String, StringBuilder>();
228
229     /*
230      * The id of the sequence being read (for non-interleaved)
231      */
232     currentSequenceId = "";
233
234     while (dataLine != null)
235     {
236       dataLine = dataLine.trim();
237       if (dataLine.length() > 0)
238       {
239         if (dataLine.startsWith(BANG + GENE)
240                 || dataLine.startsWith(BANG + DOMAIN))
241         {
242           parseGeneOrDomain(dataLine);
243         }
244         else
245         {
246           currentSequenceId = parseDataLine(dataLine);
247         }
248       }
249       else if (!seqData.isEmpty())
250       {
251         /*
252          * Blank line after processing some data...
253          */
254         endOfDataBlock();
255       }
256       dataLine = nextNonCommentLine();
257     }
258
259     /*
260      * close off any features currently being parsed
261      */
262     createFeature(GENE, currentGene, geneStart);
263     createFeature(DOMAIN, currentDomain, domainStart);
264
265     // remember the (longest) line length read in, so we can output the same
266     setAlignmentProperty(PROP_LINELENGTH, String.valueOf(positionsPerLine));
267
268     deriveSequences();
269   }
270
271   /**
272    * Post-processing after reading one block of interleaved data
273    */
274   protected void endOfDataBlock()
275   {
276     this.firstDataBlockRead = true;
277     // TODO:
278     // (initialise and) populate arrays of sequence length so far (excluding
279     // gaps)
280     // On change or end of a denoted Gene or Domain, add sequence features for
281     // it
282   }
283
284   /**
285    * Parse a !Gene or !Domain command line. MEGA accepts
286    * <ul>
287    * <li>!Gene=name;</li>
288    * <li>!Gene=name Property=Coding/Noncoding CodonStart=1/2/3;</li>
289    * <li>!Gene=genename Domain=domainname Property= etc</li>
290    * <li>!Domain=domainname Gene=genename Property= etc</li>
291    * <li>!Domain=domainname Property= etc</li>
292    * <li>!domain=domainname property=domainend</li>
293    * </ul>
294    * Properly, a Gene should be composed of Domain segments, but MEGA accepts
295    * without. Note that keywords don't seem to be case sensitive.
296    * 
297    * @param dataLine
298    * @throws FileFormatException
299    */
300   protected void parseGeneOrDomain(String dataLine)
301           throws FileFormatException
302   {
303     String domain = null;
304     String gene = null;
305     String property = null;
306     String codonStart = null;
307     String errorMsg = "Unrecognized format: " + dataLine;
308
309     if (!dataLine.startsWith(BANG) || !dataLine.endsWith(SEMICOLON))
310     {
311       throw new FileFormatException(errorMsg);
312     }
313     String trimmed = dataLine.substring(1, dataLine.length() - 1).trim();
314     String[] tokens = trimmed.split(WHITESPACE);
315     for (String token : tokens)
316     {
317       String[] keyValue = token.split("=");
318       if (keyValue.length != 2)
319       {
320         throw new FileFormatException(errorMsg);
321       }
322       String key = keyValue[0];
323       if (GENE.equalsIgnoreCase(key))
324       {
325         gene = keyValue[1];
326       }
327       else if (DOMAIN.equalsIgnoreCase(key))
328       {
329         domain = keyValue[1];
330       }
331       else if (PROPERTY.equalsIgnoreCase(key))
332       {
333         property = keyValue[1];
334       }
335       else if (CODONSTART.equalsIgnoreCase(key))
336       {
337         codonStart = keyValue[1];
338       }
339       else
340       {
341         System.err.println("Unrecognised token: '" + key + "; in "
342                 + dataLine);
343       }
344     }
345
346     processGeneOrDomain(gene, domain, property, codonStart);
347   }
348
349   /**
350    * Process a statement containing one or both of Gene and Domain, and
351    * optionally Property or CodonStart commands.
352    * 
353    * @param gene
354    *          the Gene name if specified, else null
355    * @param domain
356    *          the Domain name if specified, else null
357    * @param property
358    *          the Property value if specified, else null
359    * @param codonStart
360    *          the CodonStart value if specified, else null
361    */
362   protected void processGeneOrDomain(String gene, String domain,
363           String property, String codonStart)
364   {
365     boolean domainEnd = "domainend".equalsIgnoreCase(property);
366
367     /*
368      * If we have been processing a Domain or Gene, and this does not continue
369      * it, then close it off (generate sequence features for it). Do Domain
370      * first as it is in the context of the enclosing gene if any.
371      */
372     if (this.currentDomain != null)
373     {
374       if (!this.currentDomain.equals(domain) || domainEnd)
375       {
376         String description = currentDomain
377                 + (currentGene == null ? "" : " (" + currentGene + ")");
378         createFeature(DOMAIN, description, domainStart);
379       }
380     }
381     if (this.currentGene != null && !this.currentGene.equals(gene))
382     {
383       createFeature(GENE, currentGene, geneStart);
384     }
385
386     /*
387      * and if we have declared a Gene or Domain which does not continue the
388      * current one, then record its start positions per sequence
389      */
390     if (gene != null && !gene.equals(currentGene))
391     {
392       startSequenceFeature(geneStart);
393     }
394     if (domain != null && !domain.equals(currentDomain))
395     {
396       startSequenceFeature(domainStart);
397     }
398
399     currentGene = gene;
400     currentDomain = domainEnd ? null : domain;
401   }
402
403   /**
404    * Start processing a new feature
405    * 
406    * @param startPositions
407    */
408   protected void startSequenceFeature(Map<String, Integer> startPositions)
409   {
410     /*
411      * If the feature declaration precedes all sequences, we will know in
412      * createFeature that it started with residue 1; otherwise note now where it
413      * starts in each sequence
414      */
415     if (!residuesRead.isEmpty())
416     {
417       for (Entry<String, Integer> entry : residuesRead.entrySet())
418       {
419         String seqId = entry.getKey();
420         Integer nextResidue = entry.getValue() + 1;
421         startPositions.put(seqId, nextResidue);
422       }
423     }
424   }
425
426   /**
427    * Add a SequenceFeature to each sequence, using the given start/end values
428    * per sequence
429    * 
430    * @param featureType
431    * @param featureValue
432    * @param featureStartResidues
433    */
434   protected void createFeature(String featureType, String featureValue,
435           Map<String, Integer> featureStartResidues)
436   {
437     if (featureValue == null)
438     {
439       return;
440     }
441
442     Iterator<String> seqids = this.seqData.keySet().iterator();
443     while (seqids.hasNext())
444     {
445       String seqid = seqids.next();
446       Integer startAt = featureStartResidues.get(seqid);
447       int sfstart = startAt == null ? 1 : startAt.intValue();
448       int sfend = residuesRead.get(seqid);
449       if (sfend >= sfstart)
450       {
451         /*
452          * don't add feature if entirely gapped in the sequence
453          */
454         // TODO: type="Gene" (but then all coloured the same) or
455         // type="GeneName"?
456         SequenceFeature sf = new SequenceFeature(featureValue, featureType,
457                 sfstart, sfend, 0f, null);
458         sequenceFeatures.get(seqid).add(sf);
459       }
460     }
461   }
462
463   /**
464    * Returns the next line that is not a comment, or null at end of file.
465    * Comments in MEGA are within [ ] brackets, and may be nested.
466    * 
467    * @return
468    * @throws IOException
469    */
470   protected String nextNonCommentLine() throws IOException
471   {
472     return nextNonCommentLine(0);
473   }
474
475   /**
476    * Returns the next non-comment line (or part line), or null at end of file.
477    * Comments in MEGA are within [ ] brackets, and may be nested. They may occur
478    * anywhere within a line (for example at the end with position numbers); this
479    * method returns the line with any comments removed.
480    * 
481    * @param depth
482    *          current depth of nesting of comments while parsing
483    * @return
484    * @throws IOException
485    */
486   protected String nextNonCommentLine(final int depth) throws IOException
487   {
488     String data = null;
489     data = nextLine();
490     if (data == null)
491     {
492       if (depth > 0)
493       {
494         System.err.println("Warning: unterminated comment in data file");
495       }
496       return data;
497     }
498
499     /*
500      * If we are in a (possibly nested) comment after parsing this line, keep
501      * reading recursively until the comment has unwound
502      */
503     int newDepth = commentDepth(data, depth);
504     if (newDepth > 0)
505     {
506       return nextNonCommentLine(newDepth);
507     }
508     else
509     {
510       /*
511        * not in a comment by end of this line; return what is left
512        */
513       String nonCommentPart = getNonCommentContent(data, depth);
514       return nonCommentPart;
515     }
516   }
517
518   /**
519    * Returns what is left of the input data after removing any comments, whether
520    * 'in progress' from preceding lines, or embedded in the current line
521    * 
522    * @param data
523    *          input data
524    * @param depth
525    *          nested depth of comments pending termination
526    * @return
527    * @throws FileFormatException
528    */
529   protected static String getNonCommentContent(String data, int depth)
530           throws FileFormatException
531   {
532     int len = data.length();
533     StringBuilder result = new StringBuilder(len);
534     for (int i = 0; i < len; i++)
535     {
536       char c = data.charAt(i);
537       switch (c)
538       {
539       case COMMENT_START:
540         depth++;
541         break;
542
543       case COMMENT_END:
544         if (depth > 0)
545         {
546           depth--;
547         }
548         else
549         {
550           result.append(c);
551         }
552         break;
553
554       default:
555         if (depth == 0)
556         {
557           result.append(c);
558         }
559       }
560     }
561     return result.toString();
562   }
563
564   /**
565    * Calculates new depth of comment after parsing an input line i.e. the excess
566    * of opening '[' over closing ']' characters. Any excess ']' are ignored (not
567    * treated as comment delimiters).
568    * 
569    * @param data
570    *          input line
571    * @param depth
572    *          current comment nested depth before parsing the line
573    * @return new depth after parsing the line
574    */
575   protected static int commentDepth(CharSequence data, int depth)
576   {
577     int newDepth = depth;
578     int len = data.length();
579     for (int i = 0; i < len; i++)
580     {
581       char c = data.charAt(i);
582       if (c == COMMENT_START)
583       {
584         newDepth++;
585       }
586       else if (c == COMMENT_END && newDepth > 0)
587       {
588         newDepth--;
589       }
590     }
591     return newDepth;
592   }
593
594   /**
595    * Convert the parsed sequence strings to objects and store them in the model.
596    */
597   protected void deriveSequences()
598   {
599     Set<Entry<String, StringBuilder>> datasets = seqData.entrySet();
600
601     for (Entry<String, StringBuilder> dataset : datasets)
602     {
603       String sequenceId = dataset.getKey();
604       StringBuilder characters = dataset.getValue();
605       SequenceI s = new Sequence(sequenceId, new String(characters));
606       this.seqs.addElement(s);
607
608       /*
609        * and add any derived sequence features to the sequence
610        */
611       for (SequenceFeature sf : sequenceFeatures.get(sequenceId))
612       {
613         s.addSequenceFeature(sf);
614       }
615     }
616   }
617
618   /**
619    * Process one line of sequence data. If it has no sequence identifier, append
620    * to the current id's sequence. Else parse out the sequence id and append the
621    * data (if any) to that id's sequence. Returns the sequence id (implicit or
622    * explicit) for this line.
623    * 
624    * @param dataLine
625    * @return
626    * @throws IOException
627    */
628   protected String parseDataLine(String dataLine)
629           throws IOException
630   {
631     String seqId = getSequenceId(dataLine);
632     if (seqId == null)
633     {
634       /*
635        * Just character data
636        */
637       parseNoninterleavedDataLine(dataLine);
638       return currentSequenceId;
639     }
640     else if ((HASHSIGN + seqId).trim().equals(dataLine.trim()))
641     {
642       /*
643        * Sequence id only - header line for noninterleaved data
644        */
645       return seqId;
646     }
647     else
648     {
649       /*
650        * Sequence id followed by data
651        */
652       parseInterleavedDataLine(dataLine, seqId);
653       return seqId;
654     }
655   }
656
657   /**
658    * Add a line of sequence data to the buffer for the given sequence id. Start
659    * a new one if we haven't seen it before.
660    * 
661    * @param dataLine
662    * @throws IOException
663    */
664   protected void parseNoninterleavedDataLine(String dataLine)
665           throws IOException
666   {
667     if (currentSequenceId == null)
668     {
669       /*
670        * Oops. Data but no sequence id context.
671        */
672       throw new IOException("No sequence id context at: " + dataLine);
673     }
674
675     assertInterleaved(false, dataLine);
676
677     dataLine = addSequenceData(currentSequenceId, dataLine);
678
679     setPositionsPerLine(Math.max(positionsPerLine, dataLine.length()));
680   }
681
682   /**
683    * Get the sequence data for this sequence id, starting a new one if
684    * necessary.
685    * 
686    * @param currentId
687    * @return
688    */
689   protected StringBuilder getSequenceDataBuffer(String currentId)
690   {
691     StringBuilder sb = seqData.get(currentId);
692     if (sb == null)
693     {
694       // first data met for this sequence id, start a new buffer
695       sb = new StringBuilder(SEQBUFFERSIZE);
696       seqData.put(currentId, sb);
697
698       // and a placeholder for any SequenceFeature found
699       sequenceFeatures.put(currentId, new ArrayList<SequenceFeature>());
700     }
701     return sb;
702   }
703
704   /**
705    * Parse one line of interleaved data e.g.
706    * 
707    * <pre>
708    * #TheSeqId CGATCGCATGCA
709    * </pre>
710    * 
711    * @param dataLine
712    * @param seqId
713    * @throws FileFormatException
714    */
715   protected void parseInterleavedDataLine(String dataLine, String seqId)
716           throws FileFormatException
717   {
718     /*
719      * New sequence found in second or later data block - error.
720      */
721     if (this.firstDataBlockRead && !seqData.containsKey(seqId))
722     {
723       throw new FileFormatException(
724               "Parse error: misplaced new sequence starting at " + dataLine);
725     }
726
727     String data = dataLine.substring(seqId.length() + 1).trim();
728
729     /*
730      * Do nothing if this line is _only_ a sequence id with no data following.
731      */
732     if (data != null && data.length() > 0)
733     {
734       data = addSequenceData(seqId, data);
735       setPositionsPerLine(Math.max(positionsPerLine, data.length()));
736       assertInterleaved(true, dataLine);
737     }
738   }
739
740   /**
741    * Remove spaces, and replace identity symbol, before appending the sequence
742    * data to the buffer for the sequence id. Returns the reformatted added data.
743    * Also updates a count of residues read for the sequence.
744    * 
745    * @param seqId
746    * @param data
747    * @return
748    */
749   protected String addSequenceData(String seqId, String data)
750   {
751     StringBuilder sb = getSequenceDataBuffer(seqId);
752     int len = sb.length();
753     String formatted = data.replace(SPACE, "");
754
755     /*
756      * If sequence contains '.' or other identity symbol; replace these with the
757      * same position from the first (reference) sequence
758      */
759     int nonGapped = 0;
760     StringBuilder referenceSequence = seqData.values().iterator().next();
761     StringBuilder sb1 = new StringBuilder(formatted.length());
762     for (int i = 0; i < formatted.length(); i++)
763     {
764       char nextChar = formatted.charAt(i);
765       if (nextChar != gapCharacter)
766       {
767         nonGapped++;
768       }
769       if (nextChar == identityCharacter
770               && len + i < referenceSequence.length())
771       {
772         sb1.append(referenceSequence.charAt(len + i));
773       }
774       else
775       {
776         sb1.append(nextChar);
777       }
778     }
779     formatted = sb1.toString();
780
781     data = formatted;
782     sb.append(data);
783
784     /*
785      * increment residue count for the sequence
786      */
787     if (nonGapped > 0)
788     {
789       Integer residueCount = residuesRead.get(seqId);
790       residuesRead.put(seqId, nonGapped
791               + (residueCount == null ? 0 : residueCount));
792     }
793
794     return data;
795   }
796
797   /**
798    * If the line begins with (e.g.) "#abcde " then returns "abcde" as the
799    * identifier. Else returns null.
800    * 
801    * @param dataLine
802    * @return
803    */
804   public static String getSequenceId(String dataLine)
805   {
806     // TODO refactor to a StringUtils type class
807     if (dataLine != null)
808     {
809       if (dataLine.startsWith(HASHSIGN))
810       {
811         int spacePos = dataLine.indexOf(" ");
812         return (spacePos == -1 ? dataLine.substring(1) : dataLine
813                 .substring(1, spacePos));
814       }
815     }
816     return null;
817   }
818
819   /**
820    * Read the #MEGA and Title/Format/Description header lines (if present).
821    * 
822    * Save as alignment properties in case useful.
823    * 
824    * @return the next non-blank line following the header lines.
825    * @throws IOException
826    */
827   protected String parseHeaderLines() throws IOException
828   {
829     String inputLine = null;
830     while ((inputLine = nextNonCommentLine()) != null)
831     {
832       inputLine = inputLine.trim();
833
834       /*
835        * skip blank lines
836        */
837       if (inputLine.length() == 0)
838       {
839         continue;
840       }
841
842       if (inputLine.toUpperCase().startsWith(MEGA_ID))
843       {
844         continue;
845       }
846
847       if (isTitle(inputLine))
848       {
849         this.title = getValue(inputLine);
850         setAlignmentProperty(PROP_TITLE, title);
851       }
852       else if (inputLine.startsWith(BANG + DESCRIPTION))
853       {
854         parseDescription(inputLine);
855       }
856
857       else if (inputLine.startsWith(BANG + FORMAT))
858       {
859         parseFormat(inputLine);
860       }
861       else if (!inputLine.toUpperCase().startsWith(MEGA_ID))
862       {
863
864         /*
865          * Return the first 'data line' i.e. one that is not blank, #MEGA or
866          * TITLE:
867          */
868         break;
869       }
870     }
871     return inputLine;
872   }
873
874   /**
875    * Parse a !Format statement. This may be multiline, and is ended by a
876    * semicolon.
877    * 
878    * @param inputLine
879    * @throws IOException
880    */
881   protected void parseFormat(String inputLine) throws IOException
882   {
883     while (inputLine != null)
884     {
885       parseFormatLine(inputLine);
886       if (inputLine.endsWith(SEMICOLON))
887       {
888         break;
889       }
890       inputLine = nextNonCommentLine();
891     }
892   }
893
894   /**
895    * Parse one line of a !Format statement. This may contain one or more
896    * keyword=value pairs.
897    * 
898    * @param inputLine
899    * @throws FileFormatException
900    */
901   protected void parseFormatLine(String inputLine)
902           throws FileFormatException
903   {
904     if (inputLine.startsWith(BANG + FORMAT))
905     {
906       inputLine = inputLine.substring((BANG + FORMAT).length());
907     }
908     if (inputLine.endsWith(SEMICOLON))
909     {
910       inputLine = inputLine.substring(0, inputLine.length() - 1);
911     }
912     if (inputLine.length() == 0)
913     {
914       return;
915     }
916     String[] tokens = inputLine.trim().split(WHITESPACE);
917     for (String token : tokens)
918     {
919       parseFormatKeyword(token);
920     }
921   }
922
923   /**
924    * Parse a Keyword=Value token. Possible keywords are
925    * <ul>
926    * <li>DataType= DNA, RNA, Nucleotide, Protein</li>
927    * <li>DataFormat= Interleaved, ?</li>
928    * <li>NSeqs= number of sequences (synonym NTaxa)</li>
929    * <li>NSites= number of bases / residues</li>
930    * <li>Property= Exon (or Coding), Intron (or Noncoding), End (of domain)</li>
931    * <li>Indel= gap character</li>
932    * <li>Identical= identity character (to first sequence) (synonym MatchChar)</li>
933    * <li>Missing= missing data character</li>
934    * <li>CodeTable= Standard, other (MEGA supports various)</li>
935    * </ul>
936    * 
937    * @param token
938    * @throws FileFormatException
939    *           if an unrecognised keyword or value is encountered
940    */
941   protected void parseFormatKeyword(String token)
942           throws FileFormatException
943   {
944     String msg = "Unrecognised Format command: " + token;
945     String[] bits = token.split(EQUALS);
946     if (bits.length != 2)
947     {
948       throw new FileFormatException(msg);
949     }
950     String keyword = bits[0];
951     String value = bits[1];
952
953     /*
954      * Jalview will work out whether nucleotide or not anyway
955      */
956     if (keyword.equalsIgnoreCase(DATATYPE))
957     {
958       if (value.equalsIgnoreCase("DNA") || value.equalsIgnoreCase("RNA")
959               || value.equalsIgnoreCase("Nucleotide"))
960       {
961         this.nucleotide = true;
962         // alignment computes whether or not it is nucleotide when created
963       }
964       else if (value.equalsIgnoreCase(PROTEIN))
965       {
966         this.nucleotide = false;
967       }
968       else
969       {
970         throw new FileFormatException(msg);
971       }
972       setAlignmentProperty(PROP_DATATYPE, value);
973     }
974
975     /*
976      * accept non-Standard code table but save in case we want to disable
977      * 'translate as cDNA'
978      */
979     else if (keyword.equalsIgnoreCase(CODETABLE))
980     {
981       setAlignmentProperty(PROP_CODETABLE, value);
982     }
983
984     /*
985      * save gap char to set later on alignment once created
986      */
987     else if (keyword.equalsIgnoreCase(INDEL))
988     {
989       this.gapCharacter = value.charAt(0);
990     }
991
992     else if (keyword.equalsIgnoreCase(IDENTICAL)
993             || keyword.equalsIgnoreCase("MatchChar"))
994     {
995       setAlignmentProperty(PROP_IDENTITY, value);
996       this.identityCharacter = value.charAt(0);
997       if (!".".equals(value))
998       {
999         System.err.println("Warning: " + token
1000                 + " not supported, Jalview uses '.' for identity");
1001       }
1002     }
1003
1004     else if (keyword.equalsIgnoreCase(MISSING))
1005     {
1006       setAlignmentProperty(PROP_MISSING, value);
1007       System.err.println("Warning: " + token + " not supported");
1008     }
1009
1010     else if (keyword.equalsIgnoreCase(PROPERTY))
1011     {
1012       // TODO: can Property appear in a Format command?
1013       // suspect this is a mistake in the manual
1014     }
1015
1016     else if (!keyword.equalsIgnoreCase(N_SEQS)
1017             && !keyword.equalsIgnoreCase("NTaxa")
1018             && !keyword.equalsIgnoreCase(N_SITES))
1019     {
1020       System.err.println("Warning: " + msg);
1021     }
1022   }
1023
1024   /**
1025    * Returns the trimmed data on the line following either whitespace or '=',
1026    * with any trailing semi-colon removed<br>
1027    * So
1028    * <ul>
1029    * <li>Hello World</li>
1030    * <li>!Hello: \tWorld;</li>
1031    * <li>!Hello=World</li>
1032    * <ul>
1033    * should all return "World"
1034    * 
1035    * @param inputLine
1036    * @return
1037    */
1038   protected static String getValue(String inputLine)
1039   {
1040     if (inputLine == null)
1041     {
1042       return null;
1043     }
1044     String value = null;
1045     String s = inputLine.replaceAll("\t", " ").trim();
1046
1047     /*
1048      * KEYWORD = VALUE should return VALUE
1049      */
1050     int equalsPos = s.indexOf("=");
1051     if (equalsPos >= 0)
1052     {
1053       value = s.substring(equalsPos + 1);
1054     }
1055     else
1056     {
1057       int spacePos = s.indexOf(' ');
1058       value = spacePos == -1 ? "" : s.substring(spacePos + 1);
1059     }
1060     value = value.trim();
1061     if (value.endsWith(SEMICOLON))
1062     {
1063       value = value.substring(0, value.length() - 1).trim();
1064     }
1065     return value;
1066   }
1067
1068   /**
1069    * Returns true if the input line starts with "TITLE" or "!TITLE" (not case
1070    * sensitive). The latter is the official format, some older data file
1071    * examples have it without the !.
1072    * 
1073    * @param inputLine
1074    * @return
1075    */
1076   protected static boolean isTitle(String inputLine)
1077   {
1078     if (inputLine == null)
1079     {
1080       return false;
1081     }
1082     String upper = inputLine.toUpperCase();
1083     return (upper.startsWith(TITLE.toUpperCase()) || upper.startsWith(BANG
1084             + TITLE.toUpperCase()));
1085   }
1086
1087   /**
1088    * Reads lines until terminated by semicolon, appending each to the
1089    * Description property value.
1090    * 
1091    * @throws IOException
1092    */
1093   protected void parseDescription(String firstDescriptionLine)
1094           throws IOException
1095   {
1096     StringBuilder desc = new StringBuilder(256);
1097     desc.append(getValue(firstDescriptionLine));
1098     if (!firstDescriptionLine.endsWith(SEMICOLON))
1099     {
1100       String line = nextNonCommentLine();
1101       while (line != null)
1102       {
1103         if (line.endsWith(SEMICOLON))
1104         {
1105           desc.append(line.substring(0, line.length() - 1));
1106           break;
1107         }
1108         else if (line.length() > 0)
1109         {
1110           desc.append(line).append(newline);
1111         }
1112         line = nextNonCommentLine();
1113       }
1114     }
1115     setAlignmentProperty(PROP_DESCRIPTION, desc.toString());
1116   }
1117
1118   /**
1119    * Returns the alignment sequences in Mega format.
1120    */
1121   @Override
1122   public String print()
1123   {
1124     return MEGA_ID + newline + print(getSeqsAsArray());
1125   }
1126
1127   /**
1128    * Write out the alignment sequences in Mega format - interleaved unless
1129    * explicitly noninterleaved.
1130    */
1131   protected String print(SequenceI[] s)
1132   {
1133     String result;
1134     if (this.interleaved != null && !this.interleaved)
1135     {
1136       result = printNonInterleaved(s);
1137     }
1138     else
1139     {
1140       result = printInterleaved(s);
1141     }
1142     return result;
1143   }
1144
1145   /**
1146    * Print to string in Interleaved format - blocks of next N characters of each
1147    * sequence in turn.
1148    * 
1149    * @param s
1150    */
1151   protected String printInterleaved(SequenceI[] s)
1152   {
1153     int maxIdLength = getMaxIdLength(s);
1154     int maxSequenceLength = getMaxSequenceLength(s);
1155     int numLines = maxSequenceLength / positionsPerLine + 3; // approx
1156
1157     int numDataBlocks = (maxSequenceLength - 1) / positionsPerLine + 1;
1158     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1159     int chunksPerLine = (positionsPerLine + spaceEvery - 1) / spaceEvery;
1160
1161     /*
1162      * Roughly size a buffer to hold the whole output
1163      */
1164     StringBuilder sb = new StringBuilder(numLines
1165             * (maxIdLength + positionsPerLine + chunksPerLine + 10));
1166
1167     /*
1168      * Output as: #Seqid CGT AGC ACT ... or blocks of 10 for peptide
1169      */
1170     int from = 0;
1171     for (int i = 0; i < numDataBlocks; i++)
1172     {
1173       sb.append(newline);
1174       boolean first = true;
1175       int advancedBy = 0;
1176       for (SequenceI seq : s)
1177       {
1178         int seqFrom = from;
1179         String seqId = String.format("#%-" + maxIdLength + "s",
1180                 seq.getName());
1181
1182         /*
1183          * output next line for this sequence
1184          */
1185         sb.append(seqId);
1186         int lastPos = seqFrom + positionsPerLine; // exclusive
1187         for (int j = 0; j < chunksPerLine; j++)
1188         {
1189           char[] subSequence = seq.getSequence(seqFrom,
1190                   Math.min(lastPos, seqFrom + spaceEvery));
1191           if (subSequence.length > 0)
1192           {
1193             sb.append(SPACE).append(subSequence);
1194           }
1195           seqFrom += subSequence.length;
1196           if (first)
1197           {
1198             // all sequences should be the same length in MEGA
1199             advancedBy += subSequence.length;
1200           }
1201         }
1202         // write last position as a comment
1203         if (writePositionNumbers)
1204         {
1205           sb.append(SPACE).append(COMMENT_START).append(from + advancedBy)
1206                   .append(COMMENT_END);
1207         }
1208         sb.append(newline);
1209         first = false;
1210       }
1211       from += advancedBy;
1212     }
1213
1214     return new String(sb);
1215   }
1216
1217   /**
1218    * Outputs to string the MEGA header and any other known and relevant
1219    * alignment properties
1220    * 
1221    * @param al
1222    */
1223   protected String printHeaders(AlignmentI al)
1224   {
1225     StringBuilder sb = new StringBuilder(128);
1226     sb.append(MEGA_ID).append(newline);
1227     String propertyValue = (String) al.getProperty(PROP_TITLE);
1228     if (propertyValue != null)
1229     {
1230       sb.append(BANG).append(TITLE).append(SPACE).append(propertyValue)
1231               .append(SEMICOLON).append(newline);
1232     }
1233     propertyValue = (String) al.getProperty(PROP_DESCRIPTION);
1234     if (propertyValue != null)
1235     {
1236       sb.append(BANG).append(DESCRIPTION).append(newline)
1237               .append(propertyValue).append(SEMICOLON)
1238               .append(newline);
1239     }
1240
1241     /*
1242      * !Format DataType CodeTable
1243      */
1244     sb.append(BANG).append(FORMAT).append(newline);
1245     String dataType = (String) al.getProperty(PROP_DATATYPE);
1246     if (dataType == null)
1247     {
1248       dataType = al.isNucleotide() ? NUCLEOTIDE : PROTEIN;
1249     }
1250     sb.append(INDENT).append(DATATYPE).append(EQUALS).append(dataType);
1251     String codeTable = (String) al.getProperty(PROP_CODETABLE);
1252     sb.append(SPACE).append(CODETABLE).append(EQUALS)
1253             .append(codeTable == null ? "Standard" : codeTable)
1254             .append(newline);
1255     
1256     /*
1257      * !Format NSeqs NSites (the length of sequences - they should all be the
1258      * same - including gaps)
1259      */
1260     sb.append(INDENT).append(N_SEQS).append(EQUALS).append(al.getHeight());
1261     sb.append(SPACE).append(N_SITES).append(EQUALS)
1262             .append(String.valueOf(al.getWidth()));
1263     sb.append(newline);
1264
1265     /*
1266      * !Format Indel Identical Missing
1267      */
1268     sb.append(INDENT);
1269     sb.append(INDEL).append(EQUALS).append(al.getGapCharacter());
1270     String identity = (String) al.getProperty(PROP_IDENTITY);
1271     if (identity != null)
1272     {
1273       sb.append(SPACE).append(IDENTICAL).append(EQUALS).append(identity);
1274     }
1275     String missing = (String) al.getProperty(PROP_MISSING);
1276     if (missing != null)
1277     {
1278       sb.append(SPACE).append(MISSING).append(EQUALS).append(missing);
1279     }
1280     sb.append(SEMICOLON).append(newline);
1281
1282     return sb.toString();
1283   }
1284
1285   /**
1286    * Get the longest sequence id (to allow aligned printout).
1287    * 
1288    * @param s
1289    * @return
1290    */
1291   protected static int getMaxIdLength(SequenceI[] s)
1292   {
1293     // TODO pull up for reuse
1294     int maxLength = 0;
1295     for (SequenceI seq : s)
1296     {
1297       int len = seq.getName().length();
1298       if (len > maxLength)
1299       {
1300         maxLength = len;
1301       }
1302     }
1303     return maxLength;
1304   }
1305
1306   /**
1307    * Get the longest sequence length
1308    * 
1309    * @param s
1310    * @return
1311    */
1312   protected static int getMaxSequenceLength(SequenceI[] s)
1313   {
1314     // TODO pull up for reuse
1315     int maxLength = 0;
1316     for (SequenceI seq : s)
1317     {
1318       int len = seq.getLength();
1319       if (len > maxLength)
1320       {
1321         maxLength = len;
1322       }
1323     }
1324     return maxLength;
1325   }
1326
1327   /**
1328    * Print to string in noninterleaved format - all of each sequence in turn, in
1329    * blocks of 50 characters.
1330    * 
1331    * @param s
1332    * @return
1333    */
1334   protected String printNonInterleaved(SequenceI[] s)
1335   {
1336     int maxSequenceLength = getMaxSequenceLength(s);
1337     // approx
1338     int numLines = maxSequenceLength / positionsPerLine + 2 + s.length;
1339
1340     /*
1341      * Roughly size a buffer to hold the whole output
1342      */
1343     StringBuilder sb = new StringBuilder(numLines * positionsPerLine);
1344
1345     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1346     int chunksPerLine = positionsPerLine / spaceEvery;
1347     for (SequenceI seq : s)
1348     {
1349       sb.append(newline);
1350       sb.append(HASHSIGN + seq.getName()).append(newline);
1351       int startPos = 0;
1352       while (startPos < seq.getLength())
1353       {
1354         boolean firstChunk = true;
1355         /*
1356          * print next line for this sequence
1357          */
1358         int lastPos = startPos + positionsPerLine; // exclusive
1359         for (int j = 0; j < chunksPerLine; j++)
1360         {
1361           char[] subSequence = seq.getSequence(startPos,
1362                   Math.min(lastPos, startPos + positionsPerLine));
1363           if (subSequence.length > 0)
1364           {
1365             if (!firstChunk)
1366             {
1367               sb.append(SPACE);
1368             }
1369             sb.append(subSequence);
1370             firstChunk = false;
1371           }
1372           startPos += subSequence.length;
1373         }
1374         sb.append(newline);
1375       }
1376     }
1377
1378     return new String(sb);
1379   }
1380
1381   /**
1382    * Flag this file as interleaved or not, based on data format. Throws an
1383    * exception if has previously been determined to be otherwise.
1384    * 
1385    * @param isIt
1386    * @param dataLine
1387    * @throws IOException
1388    */
1389   protected void assertInterleaved(boolean isIt, String dataLine)
1390           throws FileFormatException
1391   {
1392     if (this.interleaved != null && isIt != this.interleaved.booleanValue())
1393     {
1394       throw new FileFormatException(
1395               "Parse error: mix of interleaved and noninterleaved detected, at line: "
1396                       + dataLine);
1397     }
1398     this.interleaved = new Boolean(isIt);
1399     setAlignmentProperty(PROP_INTERLEAVED, interleaved.toString());
1400   }
1401
1402   public boolean isInterleaved()
1403   {
1404     return this.interleaved == null ? false : this.interleaved
1405             .booleanValue();
1406   }
1407
1408   /**
1409    * Adds saved parsed values either as alignment properties, or (in some cases)
1410    * as specific member fields of the alignment
1411    */
1412   @Override
1413   public void addProperties(AlignmentI al)
1414   {
1415     super.addProperties(al);
1416     al.setGapCharacter(gapCharacter);
1417     
1418     /*
1419      * warn if e.g. DataType=DNA but data is protein (or vice versa)
1420      */
1421     if (this.nucleotide != null && this.nucleotide != al.isNucleotide()) {
1422       System.err.println("Warning: " + this.title + " declared "
1423               + (nucleotide ? "" : " not ") + "nucleotide but it is"
1424               + (nucleotide ? " not" : ""));
1425     }
1426   }
1427
1428   /**
1429    * Print the given alignment in MEGA format. If the alignment was created by
1430    * parsing a MEGA file, it should have properties set (e.g. Title) which can
1431    * influence the output.
1432    */
1433   @Override
1434   public String print(AlignmentI al)
1435   {
1436     this.nucleotide = al.isNucleotide();
1437
1438     String lineLength = (String) al.getProperty(PROP_LINELENGTH);
1439     this.positionsPerLine = lineLength == null ? DEFAULT_LINE_LENGTH : Integer
1440             .parseInt(lineLength);
1441
1442     /*
1443      * round down to a multiple of 3 positions per line for nucleotide
1444      */
1445     if (nucleotide)
1446     {
1447       positionsPerLine = positionsPerLine - (positionsPerLine % 3);
1448     }
1449
1450     String interleave = (String) al.getProperty(PROP_INTERLEAVED);
1451     if (interleave != null)
1452     {
1453       this.interleaved = Boolean.valueOf(interleave);
1454     }
1455
1456     String headers = printHeaders(al);
1457     return headers + print(al.getSequencesArray());
1458   }
1459
1460   /**
1461    * Returns the number of sequence positions output per line
1462    * 
1463    * @return
1464    */
1465   public int getPositionsPerLine()
1466   {
1467     return positionsPerLine;
1468   }
1469
1470   /**
1471    * Sets the number of sequence positions output per line. Note these will be
1472    * formatted in blocks of 3 (nucleotide) or 10 (peptide).
1473    * 
1474    * @param p
1475    */
1476   public void setPositionsPerLine(int p)
1477   {
1478     this.positionsPerLine = p;
1479   }
1480 }