58961f824c301ced6846577485b97f618e509d38
[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     /*
366      * If we have been processing a Gene or Domain, and this does not continue
367      * it, then close it off (generate sequence features for it)
368      */
369     if (this.currentGene != null && !this.currentGene.equals(gene))
370     {
371       createFeature(GENE, currentGene, geneStart);
372     }
373     if (this.currentDomain != null)
374     {
375       if (!this.currentDomain.equals(domain)
376               || "domainend".equalsIgnoreCase(property))
377       {
378         createFeature(DOMAIN, currentDomain, domainStart);
379       }
380     }
381
382     /*
383      * and if we have declared a Gene or Domain which does not continue the
384      * current one, then record its start positions per sequence
385      */
386     if (gene != null && !gene.equals(currentGene))
387     {
388       startSequenceFeature(geneStart);
389     }
390     if (domain != null && !domain.equals(currentDomain))
391     {
392       startSequenceFeature(domainStart);
393     }
394
395     currentGene = gene;
396     currentDomain = domain;
397   }
398
399   /**
400    * Start processing a new feature
401    * 
402    * @param startPositions
403    */
404   protected void startSequenceFeature(Map<String, Integer> startPositions)
405   {
406     /*
407      * If the feature declaration precedes all sequences, we will know in
408      * createFeature that it started with residue 1; otherwise note now where it
409      * starts in each sequence
410      */
411     if (!residuesRead.isEmpty())
412     {
413       for (Entry<String, Integer> entry : residuesRead.entrySet())
414       {
415         String seqId = entry.getKey();
416         Integer nextResidue = entry.getValue() + 1;
417         startPositions.put(seqId, nextResidue);
418       }
419     }
420   }
421
422   /**
423    * Add a SequenceFeature to each sequence, using the given start/end values
424    * per sequence
425    * 
426    * @param featureType
427    * @param featureValue
428    * @param featureStartResidues
429    */
430   protected void createFeature(String featureType, String featureValue,
431           Map<String, Integer> featureStartResidues)
432   {
433     if (featureValue == null)
434     {
435       return;
436     }
437
438     Iterator<String> seqids = this.seqData.keySet().iterator();
439     while (seqids.hasNext())
440     {
441       String seqid = seqids.next();
442       Integer startAt = featureStartResidues.get(seqid);
443       int sfstart = startAt == null ? 1 : startAt.intValue();
444       int sfend = residuesRead.get(seqid);
445       if (sfend >= sfstart)
446       {
447         /*
448          * don't add feature if entirely gapped in the sequence
449          */
450         // TODO: type="Gene" (but then all coloured the same) or
451         // type="GeneName"?
452         SequenceFeature sf = new SequenceFeature(featureValue, featureType,
453                 sfstart, sfend, 0f, null);
454         sequenceFeatures.get(seqid).add(sf);
455       }
456     }
457   }
458
459   /**
460    * Returns the next line that is not a comment, or null at end of file.
461    * Comments in MEGA are within [ ] brackets, and may be nested.
462    * 
463    * @return
464    * @throws IOException
465    */
466   protected String nextNonCommentLine() throws IOException
467   {
468     return nextNonCommentLine(0);
469   }
470
471   /**
472    * Returns the next non-comment line (or part line), or null at end of file.
473    * Comments in MEGA are within [ ] brackets, and may be nested. They may occur
474    * anywhere within a line (for example at the end with position numbers); this
475    * method returns the line with any comments removed.
476    * 
477    * @param depth
478    *          current depth of nesting of comments while parsing
479    * @return
480    * @throws IOException
481    */
482   protected String nextNonCommentLine(final int depth) throws IOException
483   {
484     String data = null;
485     data = nextLine();
486     if (data == null)
487     {
488       if (depth > 0)
489       {
490         System.err.println("Warning: unterminated comment in data file");
491       }
492       return data;
493     }
494
495     /*
496      * If we are in a (possibly nested) comment after parsing this line, keep
497      * reading recursively until the comment has unwound
498      */
499     int newDepth = commentDepth(data, depth);
500     if (newDepth > 0)
501     {
502       return nextNonCommentLine(newDepth);
503     }
504     else
505     {
506       /*
507        * not in a comment by end of this line; return what is left
508        */
509       String nonCommentPart = getNonCommentContent(data, depth);
510       return nonCommentPart;
511     }
512   }
513
514   /**
515    * Returns what is left of the input data after removing any comments, whether
516    * 'in progress' from preceding lines, or embedded in the current line
517    * 
518    * @param data
519    *          input data
520    * @param depth
521    *          nested depth of comments pending termination
522    * @return
523    * @throws FileFormatException
524    */
525   protected static String getNonCommentContent(String data, int depth)
526           throws FileFormatException
527   {
528     int len = data.length();
529     StringBuilder result = new StringBuilder(len);
530     for (int i = 0; i < len; i++)
531     {
532       char c = data.charAt(i);
533       switch (c)
534       {
535       case COMMENT_START:
536         depth++;
537         break;
538
539       case COMMENT_END:
540         if (depth > 0)
541         {
542           depth--;
543         }
544         else
545         {
546           result.append(c);
547         }
548         break;
549
550       default:
551         if (depth == 0)
552         {
553           result.append(c);
554         }
555       }
556     }
557     return result.toString();
558   }
559
560   /**
561    * Calculates new depth of comment after parsing an input line i.e. the excess
562    * of opening '[' over closing ']' characters. Any excess ']' are ignored (not
563    * treated as comment delimiters).
564    * 
565    * @param data
566    *          input line
567    * @param depth
568    *          current comment nested depth before parsing the line
569    * @return new depth after parsing the line
570    */
571   protected static int commentDepth(CharSequence data, int depth)
572   {
573     int newDepth = depth;
574     int len = data.length();
575     for (int i = 0; i < len; i++)
576     {
577       char c = data.charAt(i);
578       if (c == COMMENT_START)
579       {
580         newDepth++;
581       }
582       else if (c == COMMENT_END && newDepth > 0)
583       {
584         newDepth--;
585       }
586     }
587     return newDepth;
588   }
589
590   /**
591    * Convert the parsed sequence strings to objects and store them in the model.
592    */
593   protected void deriveSequences()
594   {
595     Set<Entry<String, StringBuilder>> datasets = seqData.entrySet();
596
597     for (Entry<String, StringBuilder> dataset : datasets)
598     {
599       String sequenceId = dataset.getKey();
600       StringBuilder characters = dataset.getValue();
601       SequenceI s = new Sequence(sequenceId, new String(characters));
602       this.seqs.addElement(s);
603
604       /*
605        * and add any derived sequence features to the sequence
606        */
607       for (SequenceFeature sf : sequenceFeatures.get(sequenceId))
608       {
609         s.addSequenceFeature(sf);
610       }
611     }
612   }
613
614   /**
615    * Process one line of sequence data. If it has no sequence identifier, append
616    * to the current id's sequence. Else parse out the sequence id and append the
617    * data (if any) to that id's sequence. Returns the sequence id (implicit or
618    * explicit) for this line.
619    * 
620    * @param dataLine
621    * @return
622    * @throws IOException
623    */
624   protected String parseDataLine(String dataLine)
625           throws IOException
626   {
627     String seqId = getSequenceId(dataLine);
628     if (seqId == null)
629     {
630       /*
631        * Just character data
632        */
633       parseNoninterleavedDataLine(dataLine);
634       return currentSequenceId;
635     }
636     else if ((HASHSIGN + seqId).trim().equals(dataLine.trim()))
637     {
638       /*
639        * Sequence id only - header line for noninterleaved data
640        */
641       return seqId;
642     }
643     else
644     {
645       /*
646        * Sequence id followed by data
647        */
648       parseInterleavedDataLine(dataLine, seqId);
649       return seqId;
650     }
651   }
652
653   /**
654    * Add a line of sequence data to the buffer for the given sequence id. Start
655    * a new one if we haven't seen it before.
656    * 
657    * @param dataLine
658    * @throws IOException
659    */
660   protected void parseNoninterleavedDataLine(String dataLine)
661           throws IOException
662   {
663     if (currentSequenceId == null)
664     {
665       /*
666        * Oops. Data but no sequence id context.
667        */
668       throw new IOException("No sequence id context at: " + dataLine);
669     }
670
671     assertInterleaved(false, dataLine);
672
673     dataLine = addSequenceData(currentSequenceId, dataLine);
674
675     setPositionsPerLine(Math.max(positionsPerLine, dataLine.length()));
676   }
677
678   /**
679    * Get the sequence data for this sequence id, starting a new one if
680    * necessary.
681    * 
682    * @param currentId
683    * @return
684    */
685   protected StringBuilder getSequenceDataBuffer(String currentId)
686   {
687     StringBuilder sb = seqData.get(currentId);
688     if (sb == null)
689     {
690       // first data met for this sequence id, start a new buffer
691       sb = new StringBuilder(SEQBUFFERSIZE);
692       seqData.put(currentId, sb);
693
694       // and a placeholder for any SequenceFeature found
695       sequenceFeatures.put(currentId, new ArrayList<SequenceFeature>());
696     }
697     return sb;
698   }
699
700   /**
701    * Parse one line of interleaved data e.g.
702    * 
703    * <pre>
704    * #TheSeqId CGATCGCATGCA
705    * </pre>
706    * 
707    * @param dataLine
708    * @param seqId
709    * @throws FileFormatException
710    */
711   protected void parseInterleavedDataLine(String dataLine, String seqId)
712           throws FileFormatException
713   {
714     /*
715      * New sequence found in second or later data block - error.
716      */
717     if (this.firstDataBlockRead && !seqData.containsKey(seqId))
718     {
719       throw new FileFormatException(
720               "Parse error: misplaced new sequence starting at " + dataLine);
721     }
722
723     String data = dataLine.substring(seqId.length() + 1).trim();
724
725     /*
726      * Do nothing if this line is _only_ a sequence id with no data following.
727      */
728     if (data != null && data.length() > 0)
729     {
730       data = addSequenceData(seqId, data);
731       setPositionsPerLine(Math.max(positionsPerLine, data.length()));
732       assertInterleaved(true, dataLine);
733     }
734   }
735
736   /**
737    * Remove spaces, and replace identity symbol, before appending the sequence
738    * data to the buffer for the sequence id. Returns the reformatted added data.
739    * Also updates a count of residues read for the sequence.
740    * 
741    * @param seqId
742    * @param data
743    * @return
744    */
745   protected String addSequenceData(String seqId, String data)
746   {
747     StringBuilder sb = getSequenceDataBuffer(seqId);
748     int len = sb.length();
749     String formatted = data.replace(SPACE, "");
750
751     /*
752      * If sequence contains '.' or other identity symbol; replace these with the
753      * same position from the first (reference) sequence
754      */
755     int nonGapped = 0;
756     StringBuilder referenceSequence = seqData.values().iterator().next();
757     StringBuilder sb1 = new StringBuilder(formatted.length());
758     for (int i = 0; i < formatted.length(); i++)
759     {
760       char nextChar = formatted.charAt(i);
761       if (nextChar != gapCharacter)
762       {
763         nonGapped++;
764       }
765       if (nextChar == identityCharacter
766               && len + i < referenceSequence.length())
767       {
768         sb1.append(referenceSequence.charAt(len + i));
769       }
770       else
771       {
772         sb1.append(nextChar);
773       }
774     }
775     formatted = sb1.toString();
776
777     data = formatted;
778     sb.append(data);
779
780     /*
781      * increment residue count for the sequence
782      */
783     if (nonGapped > 0)
784     {
785       Integer residueCount = residuesRead.get(seqId);
786       residuesRead.put(seqId, nonGapped
787               + (residueCount == null ? 0 : residueCount));
788     }
789
790     return data;
791   }
792
793   /**
794    * If the line begins with (e.g.) "#abcde " then returns "abcde" as the
795    * identifier. Else returns null.
796    * 
797    * @param dataLine
798    * @return
799    */
800   public static String getSequenceId(String dataLine)
801   {
802     // TODO refactor to a StringUtils type class
803     if (dataLine != null)
804     {
805       if (dataLine.startsWith(HASHSIGN))
806       {
807         int spacePos = dataLine.indexOf(" ");
808         return (spacePos == -1 ? dataLine.substring(1) : dataLine
809                 .substring(1, spacePos));
810       }
811     }
812     return null;
813   }
814
815   /**
816    * Read the #MEGA and Title/Format/Description header lines (if present).
817    * 
818    * Save as alignment properties in case useful.
819    * 
820    * @return the next non-blank line following the header lines.
821    * @throws IOException
822    */
823   protected String parseHeaderLines() throws IOException
824   {
825     String inputLine = null;
826     while ((inputLine = nextNonCommentLine()) != null)
827     {
828       inputLine = inputLine.trim();
829
830       /*
831        * skip blank lines
832        */
833       if (inputLine.length() == 0)
834       {
835         continue;
836       }
837
838       if (inputLine.toUpperCase().startsWith(MEGA_ID))
839       {
840         continue;
841       }
842
843       if (isTitle(inputLine))
844       {
845         this.title = getValue(inputLine);
846         setAlignmentProperty(PROP_TITLE, title);
847       }
848       else if (inputLine.startsWith(BANG + DESCRIPTION))
849       {
850         parseDescription(inputLine);
851       }
852
853       else if (inputLine.startsWith(BANG + FORMAT))
854       {
855         parseFormat(inputLine);
856       }
857       else if (!inputLine.toUpperCase().startsWith(MEGA_ID))
858       {
859
860         /*
861          * Return the first 'data line' i.e. one that is not blank, #MEGA or
862          * TITLE:
863          */
864         break;
865       }
866     }
867     return inputLine;
868   }
869
870   /**
871    * Parse a !Format statement. This may be multiline, and is ended by a
872    * semicolon.
873    * 
874    * @param inputLine
875    * @throws IOException
876    */
877   protected void parseFormat(String inputLine) throws IOException
878   {
879     while (inputLine != null)
880     {
881       parseFormatLine(inputLine);
882       if (inputLine.endsWith(SEMICOLON))
883       {
884         break;
885       }
886       inputLine = nextNonCommentLine();
887     }
888   }
889
890   /**
891    * Parse one line of a !Format statement. This may contain one or more
892    * keyword=value pairs.
893    * 
894    * @param inputLine
895    * @throws FileFormatException
896    */
897   protected void parseFormatLine(String inputLine)
898           throws FileFormatException
899   {
900     if (inputLine.startsWith(BANG + FORMAT))
901     {
902       inputLine = inputLine.substring((BANG + FORMAT).length());
903     }
904     if (inputLine.endsWith(SEMICOLON))
905     {
906       inputLine = inputLine.substring(0, inputLine.length() - 1);
907     }
908     if (inputLine.length() == 0)
909     {
910       return;
911     }
912     String[] tokens = inputLine.trim().split(WHITESPACE);
913     for (String token : tokens)
914     {
915       parseFormatKeyword(token);
916     }
917   }
918
919   /**
920    * Parse a Keyword=Value token. Possible keywords are
921    * <ul>
922    * <li>DataType= DNA, RNA, Nucleotide, Protein</li>
923    * <li>DataFormat= Interleaved, ?</li>
924    * <li>NSeqs= number of sequences (synonym NTaxa)</li>
925    * <li>NSites= number of bases / residues</li>
926    * <li>Property= Exon (or Coding), Intron (or Noncoding), End (of domain)</li>
927    * <li>Indel= gap character</li>
928    * <li>Identical= identity character (to first sequence) (synonym MatchChar)</li>
929    * <li>Missing= missing data character</li>
930    * <li>CodeTable= Standard, other (MEGA supports various)</li>
931    * </ul>
932    * 
933    * @param token
934    * @throws FileFormatException
935    *           if an unrecognised keyword or value is encountered
936    */
937   protected void parseFormatKeyword(String token)
938           throws FileFormatException
939   {
940     String msg = "Unrecognised Format command: " + token;
941     String[] bits = token.split(EQUALS);
942     if (bits.length != 2)
943     {
944       throw new FileFormatException(msg);
945     }
946     String keyword = bits[0];
947     String value = bits[1];
948
949     /*
950      * Jalview will work out whether nucleotide or not anyway
951      */
952     if (keyword.equalsIgnoreCase(DATATYPE))
953     {
954       if (value.equalsIgnoreCase("DNA") || value.equalsIgnoreCase("RNA")
955               || value.equalsIgnoreCase("Nucleotide"))
956       {
957         this.nucleotide = true;
958         // alignment computes whether or not it is nucleotide when created
959       }
960       else if (value.equalsIgnoreCase(PROTEIN))
961       {
962         this.nucleotide = false;
963       }
964       else
965       {
966         throw new FileFormatException(msg);
967       }
968       setAlignmentProperty(PROP_DATATYPE, value);
969     }
970
971     /*
972      * accept non-Standard code table but save in case we want to disable
973      * 'translate as cDNA'
974      */
975     else if (keyword.equalsIgnoreCase(CODETABLE))
976     {
977       setAlignmentProperty(PROP_CODETABLE, value);
978     }
979
980     /*
981      * save gap char to set later on alignment once created
982      */
983     else if (keyword.equalsIgnoreCase(INDEL))
984     {
985       this.gapCharacter = value.charAt(0);
986     }
987
988     else if (keyword.equalsIgnoreCase(IDENTICAL)
989             || keyword.equalsIgnoreCase("MatchChar"))
990     {
991       setAlignmentProperty(PROP_IDENTITY, value);
992       this.identityCharacter = value.charAt(0);
993       if (!".".equals(value))
994       {
995         System.err.println("Warning: " + token
996                 + " not supported, Jalview uses '.' for identity");
997       }
998     }
999
1000     else if (keyword.equalsIgnoreCase(MISSING))
1001     {
1002       setAlignmentProperty(PROP_MISSING, value);
1003       System.err.println("Warning: " + token + " not supported");
1004     }
1005
1006     else if (keyword.equalsIgnoreCase(PROPERTY))
1007     {
1008       // TODO: can Property appear in a Format command?
1009       // suspect this is a mistake in the manual
1010     }
1011
1012     else if (!keyword.equalsIgnoreCase(N_SEQS)
1013             && !keyword.equalsIgnoreCase("NTaxa")
1014             && !keyword.equalsIgnoreCase(N_SITES))
1015     {
1016       System.err.println("Warning: " + msg);
1017     }
1018   }
1019
1020   /**
1021    * Returns the trimmed data on the line following either whitespace or '=',
1022    * with any trailing semi-colon removed<br>
1023    * So
1024    * <ul>
1025    * <li>Hello World</li>
1026    * <li>!Hello: \tWorld;</li>
1027    * <li>!Hello=World</li>
1028    * <ul>
1029    * should all return "World"
1030    * 
1031    * @param inputLine
1032    * @return
1033    */
1034   protected static String getValue(String inputLine)
1035   {
1036     if (inputLine == null)
1037     {
1038       return null;
1039     }
1040     String value = null;
1041     String s = inputLine.replaceAll("\t", " ").trim();
1042
1043     /*
1044      * KEYWORD = VALUE should return VALUE
1045      */
1046     int equalsPos = s.indexOf("=");
1047     if (equalsPos >= 0)
1048     {
1049       value = s.substring(equalsPos + 1);
1050     }
1051     else
1052     {
1053       int spacePos = s.indexOf(' ');
1054       value = spacePos == -1 ? "" : s.substring(spacePos + 1);
1055     }
1056     value = value.trim();
1057     if (value.endsWith(SEMICOLON))
1058     {
1059       value = value.substring(0, value.length() - 1).trim();
1060     }
1061     return value;
1062   }
1063
1064   /**
1065    * Returns true if the input line starts with "TITLE" or "!TITLE" (not case
1066    * sensitive). The latter is the official format, some older data file
1067    * examples have it without the !.
1068    * 
1069    * @param inputLine
1070    * @return
1071    */
1072   protected static boolean isTitle(String inputLine)
1073   {
1074     if (inputLine == null)
1075     {
1076       return false;
1077     }
1078     String upper = inputLine.toUpperCase();
1079     return (upper.startsWith(TITLE.toUpperCase()) || upper.startsWith(BANG
1080             + TITLE.toUpperCase()));
1081   }
1082
1083   /**
1084    * Reads lines until terminated by semicolon, appending each to the
1085    * Description property value.
1086    * 
1087    * @throws IOException
1088    */
1089   protected void parseDescription(String firstDescriptionLine)
1090           throws IOException
1091   {
1092     StringBuilder desc = new StringBuilder(256);
1093     desc.append(getValue(firstDescriptionLine));
1094     if (!firstDescriptionLine.endsWith(SEMICOLON))
1095     {
1096       String line = nextNonCommentLine();
1097       while (line != null)
1098       {
1099         if (line.endsWith(SEMICOLON))
1100         {
1101           desc.append(line.substring(0, line.length() - 1));
1102           break;
1103         }
1104         else if (line.length() > 0)
1105         {
1106           desc.append(line).append(newline);
1107         }
1108         line = nextNonCommentLine();
1109       }
1110     }
1111     setAlignmentProperty(PROP_DESCRIPTION, desc.toString());
1112   }
1113
1114   /**
1115    * Returns the alignment sequences in Mega format.
1116    */
1117   @Override
1118   public String print()
1119   {
1120     return MEGA_ID + newline + print(getSeqsAsArray());
1121   }
1122
1123   /**
1124    * Write out the alignment sequences in Mega format - interleaved unless
1125    * explicitly noninterleaved.
1126    */
1127   protected String print(SequenceI[] s)
1128   {
1129     String result;
1130     if (this.interleaved != null && !this.interleaved)
1131     {
1132       result = printNonInterleaved(s);
1133     }
1134     else
1135     {
1136       result = printInterleaved(s);
1137     }
1138     return result;
1139   }
1140
1141   /**
1142    * Print to string in Interleaved format - blocks of next N characters of each
1143    * sequence in turn.
1144    * 
1145    * @param s
1146    */
1147   protected String printInterleaved(SequenceI[] s)
1148   {
1149     int maxIdLength = getMaxIdLength(s);
1150     int maxSequenceLength = getMaxSequenceLength(s);
1151     int numLines = maxSequenceLength / positionsPerLine + 3; // approx
1152
1153     int numDataBlocks = (maxSequenceLength - 1) / positionsPerLine + 1;
1154     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1155     int chunksPerLine = (positionsPerLine + spaceEvery - 1) / spaceEvery;
1156
1157     /*
1158      * Roughly size a buffer to hold the whole output
1159      */
1160     StringBuilder sb = new StringBuilder(numLines
1161             * (maxIdLength + positionsPerLine + chunksPerLine + 10));
1162
1163     /*
1164      * Output as: #Seqid CGT AGC ACT ... or blocks of 10 for peptide
1165      */
1166     int from = 0;
1167     for (int i = 0; i < numDataBlocks; i++)
1168     {
1169       sb.append(newline);
1170       boolean first = true;
1171       int advancedBy = 0;
1172       for (SequenceI seq : s)
1173       {
1174         int seqFrom = from;
1175         String seqId = String.format("#%-" + maxIdLength + "s",
1176                 seq.getName());
1177
1178         /*
1179          * output next line for this sequence
1180          */
1181         sb.append(seqId);
1182         int lastPos = seqFrom + positionsPerLine; // exclusive
1183         for (int j = 0; j < chunksPerLine; j++)
1184         {
1185           char[] subSequence = seq.getSequence(seqFrom,
1186                   Math.min(lastPos, seqFrom + spaceEvery));
1187           if (subSequence.length > 0)
1188           {
1189             sb.append(SPACE).append(subSequence);
1190           }
1191           seqFrom += subSequence.length;
1192           if (first)
1193           {
1194             // all sequences should be the same length in MEGA
1195             advancedBy += subSequence.length;
1196           }
1197         }
1198         // write last position as a comment
1199         if (writePositionNumbers)
1200         {
1201           sb.append(SPACE).append(COMMENT_START).append(from + advancedBy)
1202                   .append(COMMENT_END);
1203         }
1204         sb.append(newline);
1205         first = false;
1206       }
1207       from += advancedBy;
1208     }
1209
1210     return new String(sb);
1211   }
1212
1213   /**
1214    * Outputs to string the MEGA header and any other known and relevant
1215    * alignment properties
1216    * 
1217    * @param al
1218    */
1219   protected String printHeaders(AlignmentI al)
1220   {
1221     StringBuilder sb = new StringBuilder(128);
1222     sb.append(MEGA_ID).append(newline);
1223     String propertyValue = (String) al.getProperty(PROP_TITLE);
1224     if (propertyValue != null)
1225     {
1226       sb.append(BANG).append(TITLE).append(SPACE).append(propertyValue)
1227               .append(SEMICOLON).append(newline);
1228     }
1229     propertyValue = (String) al.getProperty(PROP_DESCRIPTION);
1230     if (propertyValue != null)
1231     {
1232       sb.append(BANG).append(DESCRIPTION).append(newline)
1233               .append(propertyValue).append(SEMICOLON)
1234               .append(newline);
1235     }
1236
1237     /*
1238      * !Format DataType CodeTable
1239      */
1240     sb.append(BANG).append(FORMAT).append(newline);
1241     String dataType = (String) al.getProperty(PROP_DATATYPE);
1242     if (dataType == null)
1243     {
1244       dataType = al.isNucleotide() ? NUCLEOTIDE : PROTEIN;
1245     }
1246     sb.append(INDENT).append(DATATYPE).append(EQUALS).append(dataType);
1247     String codeTable = (String) al.getProperty(PROP_CODETABLE);
1248     sb.append(SPACE).append(CODETABLE).append(EQUALS)
1249             .append(codeTable == null ? "Standard" : codeTable)
1250             .append(newline);
1251     
1252     /*
1253      * !Format NSeqs NSites (the length of sequences - they should all be the
1254      * same - including gaps)
1255      */
1256     sb.append(INDENT).append(N_SEQS).append(EQUALS).append(al.getHeight());
1257     sb.append(SPACE).append(N_SITES).append(EQUALS)
1258             .append(String.valueOf(al.getWidth()));
1259     sb.append(newline);
1260
1261     /*
1262      * !Format Indel Identical Missing
1263      */
1264     sb.append(INDENT);
1265     sb.append(INDEL).append(EQUALS).append(al.getGapCharacter());
1266     String identity = (String) al.getProperty(PROP_IDENTITY);
1267     if (identity != null)
1268     {
1269       sb.append(SPACE).append(IDENTICAL).append(EQUALS).append(identity);
1270     }
1271     String missing = (String) al.getProperty(PROP_MISSING);
1272     if (missing != null)
1273     {
1274       sb.append(SPACE).append(MISSING).append(EQUALS).append(missing);
1275     }
1276     sb.append(SEMICOLON).append(newline);
1277
1278     return sb.toString();
1279   }
1280
1281   /**
1282    * Get the longest sequence id (to allow aligned printout).
1283    * 
1284    * @param s
1285    * @return
1286    */
1287   protected static int getMaxIdLength(SequenceI[] s)
1288   {
1289     // TODO pull up for reuse
1290     int maxLength = 0;
1291     for (SequenceI seq : s)
1292     {
1293       int len = seq.getName().length();
1294       if (len > maxLength)
1295       {
1296         maxLength = len;
1297       }
1298     }
1299     return maxLength;
1300   }
1301
1302   /**
1303    * Get the longest sequence length
1304    * 
1305    * @param s
1306    * @return
1307    */
1308   protected static int getMaxSequenceLength(SequenceI[] s)
1309   {
1310     // TODO pull up for reuse
1311     int maxLength = 0;
1312     for (SequenceI seq : s)
1313     {
1314       int len = seq.getLength();
1315       if (len > maxLength)
1316       {
1317         maxLength = len;
1318       }
1319     }
1320     return maxLength;
1321   }
1322
1323   /**
1324    * Print to string in noninterleaved format - all of each sequence in turn, in
1325    * blocks of 50 characters.
1326    * 
1327    * @param s
1328    * @return
1329    */
1330   protected String printNonInterleaved(SequenceI[] s)
1331   {
1332     int maxSequenceLength = getMaxSequenceLength(s);
1333     // approx
1334     int numLines = maxSequenceLength / positionsPerLine + 2 + s.length;
1335
1336     /*
1337      * Roughly size a buffer to hold the whole output
1338      */
1339     StringBuilder sb = new StringBuilder(numLines * positionsPerLine);
1340
1341     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1342     int chunksPerLine = positionsPerLine / spaceEvery;
1343     for (SequenceI seq : s)
1344     {
1345       sb.append(newline);
1346       sb.append(HASHSIGN + seq.getName()).append(newline);
1347       int startPos = 0;
1348       while (startPos < seq.getLength())
1349       {
1350         boolean firstChunk = true;
1351         /*
1352          * print next line for this sequence
1353          */
1354         int lastPos = startPos + positionsPerLine; // exclusive
1355         for (int j = 0; j < chunksPerLine; j++)
1356         {
1357           char[] subSequence = seq.getSequence(startPos,
1358                   Math.min(lastPos, startPos + positionsPerLine));
1359           if (subSequence.length > 0)
1360           {
1361             if (!firstChunk)
1362             {
1363               sb.append(SPACE);
1364             }
1365             sb.append(subSequence);
1366             firstChunk = false;
1367           }
1368           startPos += subSequence.length;
1369         }
1370         sb.append(newline);
1371       }
1372     }
1373
1374     return new String(sb);
1375   }
1376
1377   /**
1378    * Flag this file as interleaved or not, based on data format. Throws an
1379    * exception if has previously been determined to be otherwise.
1380    * 
1381    * @param isIt
1382    * @param dataLine
1383    * @throws IOException
1384    */
1385   protected void assertInterleaved(boolean isIt, String dataLine)
1386           throws FileFormatException
1387   {
1388     if (this.interleaved != null && isIt != this.interleaved.booleanValue())
1389     {
1390       throw new FileFormatException(
1391               "Parse error: mix of interleaved and noninterleaved detected, at line: "
1392                       + dataLine);
1393     }
1394     this.interleaved = new Boolean(isIt);
1395     setAlignmentProperty(PROP_INTERLEAVED, interleaved.toString());
1396   }
1397
1398   public boolean isInterleaved()
1399   {
1400     return this.interleaved == null ? false : this.interleaved
1401             .booleanValue();
1402   }
1403
1404   /**
1405    * Adds saved parsed values either as alignment properties, or (in some cases)
1406    * as specific member fields of the alignment
1407    */
1408   @Override
1409   public void addProperties(AlignmentI al)
1410   {
1411     super.addProperties(al);
1412     al.setGapCharacter(gapCharacter);
1413     
1414     /*
1415      * warn if e.g. DataType=DNA but data is protein (or vice versa)
1416      */
1417     if (this.nucleotide != null && this.nucleotide != al.isNucleotide()) {
1418       System.err.println("Warning: " + this.title + " declared "
1419               + (nucleotide ? "" : " not ") + "nucleotide but it is"
1420               + (nucleotide ? " not" : ""));
1421     }
1422   }
1423
1424   /**
1425    * Print the given alignment in MEGA format. If the alignment was created by
1426    * parsing a MEGA file, it should have properties set (e.g. Title) which can
1427    * influence the output.
1428    */
1429   @Override
1430   public String print(AlignmentI al)
1431   {
1432     this.nucleotide = al.isNucleotide();
1433
1434     String lineLength = (String) al.getProperty(PROP_LINELENGTH);
1435     this.positionsPerLine = lineLength == null ? DEFAULT_LINE_LENGTH : Integer
1436             .parseInt(lineLength);
1437
1438     /*
1439      * round down to a multiple of 3 positions per line for nucleotide
1440      */
1441     if (nucleotide)
1442     {
1443       positionsPerLine = positionsPerLine - (positionsPerLine % 3);
1444     }
1445
1446     String interleave = (String) al.getProperty(PROP_INTERLEAVED);
1447     if (interleave != null)
1448     {
1449       this.interleaved = Boolean.valueOf(interleave);
1450     }
1451
1452     String headers = printHeaders(al);
1453     return headers + print(al.getSequencesArray());
1454   }
1455
1456   /**
1457    * Returns the number of sequence positions output per line
1458    * 
1459    * @return
1460    */
1461   public int getPositionsPerLine()
1462   {
1463     return positionsPerLine;
1464   }
1465
1466   /**
1467    * Sets the number of sequence positions output per line. Note these will be
1468    * formatted in blocks of 3 (nucleotide) or 10 (peptide).
1469    * 
1470    * @param p
1471    */
1472   public void setPositionsPerLine(int p)
1473   {
1474     this.positionsPerLine = p;
1475   }
1476 }