JAL-1499 Gene and Domain parsed to AlignmentAnnotation (currently as
[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.AlignmentAnnotation;
22 import jalview.datamodel.AlignmentI;
23 import jalview.datamodel.Annotation;
24 import jalview.datamodel.Sequence;
25 import jalview.datamodel.SequenceFeature;
26 import jalview.datamodel.SequenceI;
27 import jalview.schemes.UserColourScheme;
28 import jalview.util.Comparison;
29
30 import java.awt.Color;
31 import java.io.IOException;
32 import java.util.ArrayList;
33 import java.util.HashMap;
34 import java.util.Iterator;
35 import java.util.LinkedHashMap;
36 import java.util.List;
37 import java.util.Map;
38 import java.util.Map.Entry;
39 import java.util.Set;
40 import java.util.Vector;
41
42 /**
43  * A parser for input or output of MEGA format files. <br>
44  * <br>
45  * Tamura K, Stecher G, Peterson D, Filipski A, and Kumar S (2013) MEGA6:
46  * Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology and
47  * Evolution 30: 2725-2729. <br>
48  * <br>
49  * 
50  * MEGA file format is supported as described in
51  * http://www.megasoftware.net/manual.pdf <br>
52  * Limitations:
53  * <ul>
54  * <li>any comments (delimited by [ ]) are ignored and not preserved</li>
55  * </ul>
56  * 
57  * @see http://www.megasoftware.net/
58  */
59 public class MegaFile extends AlignFile
60 {
61   private static final String MEGA = "MEGA";
62
63   private static final String MEGA_ANNOTATION_LABEL = MEGA + " Label";
64
65   private static final String MEGA_ANNOTATION_GENE = MEGA + " Gene";
66
67   private static final String MEGA_ANNOTATION_DOMAIN = MEGA + " Domain";
68
69   private static final char UNDERSCORE = '_';
70
71   private static final String WHITESPACE = "\\s+";
72
73   private static final int DEFAULT_LINE_LENGTH = 60;
74
75   private static final String INDENT = "    ";
76
77   private static final String N_SITES = "NSites";
78
79   private static final String N_SEQS = "NSeqs";
80
81   private static final String MISSING = "Missing";
82
83   private static final String IDENTICAL = "Identical";
84
85   private static final String INDEL = "Indel";
86
87   private static final String CODETABLE = "CodeTable";
88
89   private static final String PROTEIN = "Protein";
90
91   private static final String NUCLEOTIDE = "Nucleotide";
92
93   private static final String DATATYPE = "DataType";
94
95   private static final char COMMENT_START = '[';
96
97   private static final char COMMENT_END = ']';
98
99   private static final String HASHSIGN = "#";
100
101   private static final String SEMICOLON = ";";
102
103   private static final String BANG = "!";
104
105   private static final String EQUALS = "=";
106
107   private static final String MEGA_ID = HASHSIGN + MEGA;
108
109   private static final String TITLE = "Title";
110
111   private static final String FORMAT = "Format";
112
113   private static final String DESCRIPTION = "Description";
114
115   private static final String GENE = "Gene";
116
117   private static final String DOMAIN = "Domain";
118
119   private static final String PROPERTY = "Property";
120
121   private static final String CODONSTART = "CodonStart";
122
123   private static final String DOMAINEND = "domainend";
124
125   private static final String LABEL = "Label";
126
127   /*
128    * names of properties to save to the alignment (may affect eventual output
129    * format)
130    */
131   static final String PROP_TITLE = "MEGA_TITLE";
132
133   static final String PROP_INTERLEAVED = "MEGA_INTERLEAVED";
134
135   static final String PROP_DESCRIPTION = "MEGA_DESCRIPTION";
136
137   static final String PROP_CODETABLE = "MEGA_CODETABLE";
138
139   static final String PROP_IDENTITY = "MEGA_IDENTITY";
140
141   static final String PROP_MISSING = "MEGA_MISSING";
142
143   static final String PROP_DATATYPE = "MEGA_DATATYPE";
144
145   // number of bases per line of file (value is inferred)
146   static final String PROP_LINELENGTH = "MEGA_LINELENGTH";
147
148   // TODO: need a controlled name for Gene as a feature if we want to be able to
149   // output the MEGA file with !Gene headers
150   // WTF do we do if the sequences get realigned?
151
152   // initial size for sequence data buffer
153   private static final int SEQBUFFERSIZE = 256;
154
155   private static final String SPACE = " ";
156
157   private static final String TAB = "\t";
158
159   private static final char DEFAULT_GAP = '-';
160
161   /*
162    * number of sequence positions output per line
163    */
164   private int positionsPerLine;
165
166   private String title;
167
168   // gap character may be explicitly declared, default is -
169   private char gapCharacter = DEFAULT_GAP;
170
171   // identity character if declared
172   private char identityCharacter = 0;
173
174   // this can be True, False or null (meaning not asserted in file)
175   private Boolean nucleotide;
176
177   // set once we have seen one block of interleaved data
178   private boolean firstDataBlockRead = false;
179
180   // this can be True, False or null (meaning we don't know yet)
181   private Boolean interleaved;
182
183   // write end of line positions as a comment
184   private boolean writePositionNumbers = true;
185
186   // id of sequence being processed
187   private String currentSequenceId;
188
189   /*
190    * Temporary store of {sequenceId, positionData} while parsing interleaved
191    * sequences; sequences are maintained in the order in which they are added
192    * i.e. read in the file
193    */
194   Map<String, StringBuilder> seqData;
195   
196   // number of residues read (so far) per sequence
197   Map<String, Integer> residuesRead;
198   
199   // current Gene if any we are parsing
200   private String currentGene;
201
202   // start position in alignment (base 0) of current Gene
203   private int currentGeneStartCol;
204
205   // start residue (base 1) per sequence of current Gene
206   private Map<String, Integer> geneStart;
207
208   // current Domain if any we are parsing
209   private String currentDomain;
210
211   // start position in alignment (base 0) of current Domain
212   private int currentDomainStartCol;
213
214   // start residue (base 1) per sequence of current Domain
215   private Map<String, Integer> domainStart;
216
217   // map of SequenceFeature's by sequence id
218   private Map<String, List<SequenceFeature>> sequenceFeatures;
219
220   // each !Label line character becomes an Annotation (except underscores)
221   List<Annotation> labelAnnotations;
222
223   // records any declared Gene positions (including null values)
224   List<Annotation> geneAnnotations;
225
226   // records any declared Domain positions (including null values)
227   List<Annotation> domainAnnotations;
228
229   public MegaFile()
230   {
231   }
232
233   public MegaFile(String inFile, String type) throws IOException
234   {
235     super(inFile, type);
236   }
237
238   public MegaFile(FileParse source) throws IOException
239   {
240     super(source);
241   }
242
243   /**
244    * Parse the input stream.
245    */
246   @Override
247   public void parse() throws IOException
248   {
249     gapCharacter = DEFAULT_GAP;
250     sequenceFeatures = new HashMap<String, List<SequenceFeature>>();
251     geneStart = new HashMap<String, Integer>();
252     domainStart = new HashMap<String, Integer>();
253     residuesRead = new HashMap<String, Integer>();
254     labelAnnotations = new ArrayList<Annotation>();
255     geneAnnotations = new ArrayList<Annotation>();
256     domainAnnotations = new ArrayList<Annotation>();
257     currentDomainStartCol = -1;
258     currentGeneStartCol = -1;
259
260     /*
261      * Read and process MEGA and Title/Format/Description headers if present.
262      * Returns the first data line following the headers.
263      */
264     String dataLine = parseHeaderLines();
265
266     /*
267      * order-preserving map to hold sequences by id as they are built up during
268      * parsing
269      */
270     seqData = new LinkedHashMap<String, StringBuilder>();
271
272     /*
273      * The id of the sequence being read (for non-interleaved)
274      */
275     currentSequenceId = "";
276
277     while (dataLine != null)
278     {
279       dataLine = dataLine.trim();
280       if (dataLine.length() > 0)
281       {
282         dataLine = dataLine.replace(TAB, SPACE);
283         String upperCased = dataLine.toUpperCase();
284         if (upperCased.startsWith(BANG + GENE.toUpperCase())
285                 || upperCased.startsWith(BANG + DOMAIN.toUpperCase()))
286         {
287           parseGeneOrDomain(dataLine);
288         }
289         else if (upperCased.startsWith(BANG + LABEL.toUpperCase()))
290         {
291           parseLabel(dataLine);
292         }
293         else
294         {
295           currentSequenceId = parseDataLine(dataLine);
296         }
297       }
298       else if (!seqData.isEmpty())
299       {
300         /*
301          * Blank line after processing some data...
302          */
303         endOfDataBlock();
304       }
305       dataLine = nextNonCommentLine();
306     }
307
308     /*
309      * close off any features currently being parsed
310      */
311     createFeature(GENE, currentGene, geneStart);
312     createFeature(DOMAIN, currentDomain, domainStart);
313
314     extendAnnotation(geneAnnotations, currentGene, currentGeneStartCol);
315     extendAnnotation(domainAnnotations, currentDomain,
316             currentDomainStartCol);
317
318     // remember the (longest) line length read in, so we can output the same
319     setAlignmentProperty(PROP_LINELENGTH, String.valueOf(positionsPerLine));
320
321     deriveSequencesAndFeatures();
322
323     deriveAnnotations();
324   }
325
326   /**
327    * Create AlignmentAnnotation for Label, Gene and Domain (provided at least
328    * one non-null annotation is present)
329    */
330   protected void deriveAnnotations()
331   {
332     deriveAnnotation(this.labelAnnotations, MEGA_ANNOTATION_LABEL);
333     deriveAnnotation(this.geneAnnotations, MEGA_ANNOTATION_GENE);
334     deriveAnnotation(this.domainAnnotations, MEGA_ANNOTATION_DOMAIN);
335   }
336
337   /**
338    * Create and ad an AlignmentAnnotation (provided at least one non-null
339    * annotation is present)
340    * 
341    * @param anns
342    * @param label
343    */
344   protected void deriveAnnotation(List<Annotation> anns, String label)
345   {
346     if (anns.size() > 0 && hasNonNullEntry(anns))
347     {
348       Annotation[] annotationArray = anns.toArray(new Annotation[anns
349               .size()]);
350       AlignmentAnnotation aa = new AlignmentAnnotation(label, "",
351               annotationArray);
352       this.annotations.add(aa);
353     }
354   }
355
356   protected static boolean hasNonNullEntry(List<? extends Object> l)
357   {
358     for (Object o : l)
359     {
360       if (o != null)
361       {
362         return true;
363       }
364     }
365     return false;
366   }
367
368   /**
369    * Parse a !Label line. This contains a single character per position (column)
370    * of the alignment block above. An underscore character represents no label.
371    * Labels are assembled into an AlignmentAnnotation object.
372    * 
373    * @param dataLine
374    * @return true if any non-null annotation was created
375    * @throws FileFormatException
376    */
377   protected boolean parseLabel(String dataLine) throws FileFormatException
378   {
379     // strip off leading !Label and following spaces
380     dataLine = dataLine.substring(LABEL.length() + 1).trim();
381
382     // remove internal spacing and any leading tab
383     String labels = dataLine.replace(SPACE, "");
384     if (labels.endsWith(SEMICOLON))
385     {
386       labels = labels.substring(0, labels.length() - 1);
387     }
388     else
389     {
390       System.err.println("Warning: '" + dataLine
391               + "' should end with semi-colon");
392     }
393     boolean added = false;
394     for (char c : labels.toCharArray())
395     {
396       if (c == UNDERSCORE)
397       {
398         this.labelAnnotations.add(null);
399       }
400       else
401       {
402         this.labelAnnotations.add(new Annotation(String.valueOf(c), "",
403                 ' ', 0f));
404         added = true;
405       }
406     }
407     return added;
408   }
409
410   /**
411    * Post-processing after reading one block of interleaved data
412    */
413   protected void endOfDataBlock()
414   {
415     this.firstDataBlockRead = true;
416
417     padAnnotations(labelAnnotations);
418   }
419
420   /**
421    * Append null annotations to keep the annotations list the same length as the
422    * sequences. This ensures that when the list is converted to an array it is
423    * correctly aligned with the alignment columns. It is needed when there are
424    * gaps in declared 'annotations' in a MEGA file, such as lines with no !Label
425    * statement, or regions between marked genes or domains.
426    * 
427    * @param anns
428    */
429   protected void padAnnotations(List<Annotation> anns)
430   {
431     addNullAnnotations(anns, getAlignmentWidth());
432   }
433
434   /**
435    * Append null annotations for positions up to (and excluding) the given end
436    * column (base 0)
437    * 
438    * @param anns
439    * @param upTo
440    */
441   protected void addNullAnnotations(List<Annotation> anns, int upTo)
442   {
443     int annotationsToAdd = upTo - anns.size();
444     for (int i = 0; i < annotationsToAdd; i++)
445     {
446       anns.add(null);
447     }
448   }
449
450   /**
451    * Parse a !Gene or !Domain command line. MEGA accepts
452    * <ul>
453    * <li>!Gene=name;</li>
454    * <li>!Gene=name Property=Coding/Noncoding CodonStart=1/2/3;</li>
455    * <li>!Gene=genename Domain=domainname Property= etc</li>
456    * <li>!Domain=domainname Gene=genename Property= etc</li>
457    * <li>!Domain=domainname Property= etc</li>
458    * <li>!domain=domainname property=domainend</li>
459    * </ul>
460    * Properly, a Gene should be composed of Domain segments, but MEGA accepts
461    * without. Note that keywords don't seem to be case sensitive.
462    * 
463    * @param dataLine
464    * @throws FileFormatException
465    */
466   protected void parseGeneOrDomain(String dataLine)
467           throws FileFormatException
468   {
469     String domain = null;
470     String gene = null;
471     String property = null;
472     String codonStart = null;
473     String errorMsg = "Unrecognized format: " + dataLine;
474
475     if (!dataLine.startsWith(BANG) || !dataLine.endsWith(SEMICOLON))
476     {
477       throw new FileFormatException(errorMsg);
478     }
479     String trimmed = dataLine.substring(1, dataLine.length() - 1).trim();
480     String[] tokens = trimmed.split(WHITESPACE);
481     for (String token : tokens)
482     {
483       String[] keyValue = token.split("=");
484       if (keyValue.length != 2)
485       {
486         throw new FileFormatException(errorMsg);
487       }
488       String key = keyValue[0];
489       if (GENE.equalsIgnoreCase(key))
490       {
491         gene = keyValue[1];
492       }
493       else if (DOMAIN.equalsIgnoreCase(key))
494       {
495         domain = keyValue[1];
496       }
497       else if (PROPERTY.equalsIgnoreCase(key))
498       {
499         property = keyValue[1];
500       }
501       else if (CODONSTART.equalsIgnoreCase(key))
502       {
503         codonStart = keyValue[1];
504       }
505       else
506       {
507         System.err.println("Unrecognised token: '" + key + "; in "
508                 + dataLine);
509       }
510     }
511
512     processGeneOrDomain(gene, domain, property, codonStart);
513   }
514
515   /**
516    * Process a statement containing one or both of Gene and Domain, and
517    * optionally Property or CodonStart commands.
518    * 
519    * @param gene
520    *          the Gene name if specified, else null
521    * @param domain
522    *          the Domain name if specified, else null
523    * @param property
524    *          the Property value if specified, else null
525    * @param codonStart
526    *          the CodonStart value if specified, else null
527    */
528   protected void processGeneOrDomain(String gene, String domain,
529           String property, String codonStart)
530   {
531     /*
532      * the order of processing below ensures that we correctly handle a domain
533      * in the context of an enclosing gene
534      */
535     processDomainEnd(domain, gene, property);
536
537     processGeneEnd(gene);
538
539     processGeneStart(gene);
540
541     processDomainStart(domain, property);
542
543     // TODO save codonStart if we plan to involve it in 'translate as cDNA'
544   }
545
546   /**
547    * If we have declared a domain, and it is not continuing, start a sequence
548    * feature for it
549    * 
550    * @param domain
551    * @param property
552    */
553   protected void processDomainStart(String domain, String property)
554   {
555     if (DOMAINEND.equalsIgnoreCase(property))
556     {
557       currentDomain = null;
558       return;
559     }
560
561     if (domain != null && !domain.equals(currentDomain))
562     {
563       String verboseDomain = makeVerboseDomainName(domain, property);
564       startSequenceFeature(domainStart);
565       currentDomainStartCol = getAlignmentWidth();
566
567       currentDomain = verboseDomain;
568     }
569   }
570
571   /**
572    * Returns the width of alignment parsed so far. Note we assume (as does MEGA)
573    * that all sequences are the same length, so we can just take the length of
574    * the first one.
575    * 
576    * @return
577    */
578   protected int getAlignmentWidth()
579   {
580     return seqData.isEmpty() ? 0 : seqData.values()
581             .iterator().next().length();
582   }
583
584   /**
585    * If we have declared a gene, and it is not continuing, start a sequence
586    * feature for it
587    * 
588    * @param gene
589    */
590   protected void processGeneStart(String gene)
591   {
592     if (gene != null && !gene.equals(currentGene))
593     {
594       startSequenceFeature(geneStart);
595       currentGeneStartCol = getAlignmentWidth();
596     }
597     currentGene = gene;
598   }
599
600   /**
601    * If we have been processing a domain, and it is not being continued, then
602    * make a sequence feature for the domain just ended. Criteria for the domain
603    * not being continued are either an explicit new domain or gene name, or a
604    * 'Property=domainend' statement
605    * 
606    * @param domain
607    * @param gene
608    * @param property
609    * @return true if a feature is created, else false
610    */
611   protected boolean processDomainEnd(String domain, String gene,
612           String property)
613   {
614     boolean newGene = (gene != null && !gene.equals(currentGene));
615
616     String verboseDomain = makeVerboseDomainName(domain, property);
617
618     if (this.currentDomain != null)
619     {
620       boolean newDomain = !this.currentDomain.equals(verboseDomain);
621       boolean domainEnded = DOMAINEND.equalsIgnoreCase(property);
622       if (newDomain || newGene || domainEnded)
623       {
624         createFeature(DOMAIN, currentDomain, domainStart);
625         // and/or... create annnotations for domain
626         extendAnnotation(domainAnnotations, currentDomain,
627                 currentDomainStartCol);
628         currentDomain = null;
629         currentDomainStartCol = -1;
630         return true;
631       }
632     }
633     return false;
634   }
635
636   /**
637    * If we have been processing a gene, and it is not being continued, then make
638    * a sequence feature for the gene just ended
639    * 
640    * @param gene
641    * @return true if a feature is created, else false
642    */
643   protected boolean processGeneEnd(String gene)
644   {
645     boolean created = false;
646     /*
647      * If we were processing a gene and now have either another, or none, create
648      * a sequence feature for that gene
649      */
650     if (this.currentGene != null && !this.currentGene.equals(gene))
651     {
652       createFeature(GENE, currentGene, geneStart);
653       // and/or... add annotations for Gene
654       extendAnnotation(geneAnnotations, currentGene, currentGeneStartCol);
655       currentGene = null;
656       currentGeneStartCol = -1;
657       created = true;
658     }
659
660     return created;
661   }
662
663   /**
664    * Helper method to add Annotation elements, with the given description and
665    * starting at the given start column, up to the end of the sequence length
666    * parsed so far. Null elements are inserted for any skipped columns since the
667    * last annotation (if any), i.e. positions with no annotation of this type.
668    * 
669    * @param anns
670    * @param description
671    * @param startColumn
672    *          the start column of the annotations to add, or -1 if nothing to
673    *          add
674    */
675   protected void extendAnnotation(List<Annotation> anns,
676           String description, int startColumn)
677   {
678     int alignmentWidth = getAlignmentWidth();
679     addNullAnnotations(anns, startColumn == -1 ? alignmentWidth
680             : startColumn);
681
682     int numberToAdd = alignmentWidth - anns.size();
683     if (numberToAdd > 0)
684     {
685       Color col = description == null ? Color.black : UserColourScheme
686               .createColourFromName(description);
687       for (int i = 0; i < numberToAdd; i++)
688       {
689         anns.add(new Annotation("X", description, ' ', 0f, col));
690       }
691     }
692   }
693
694   /**
695    * Makes an expanded descriptive name for Domain if possible e.g.
696    * "Intron1 (Adh Coding)". Currently incorporates the current gene name (if
697    * any) and the Coding/Noncoding property value (if given).
698    * 
699    * @param domain
700    * @param property
701    * @return
702    */
703   protected String makeVerboseDomainName(String domain, String property)
704   {
705     String verboseDomain = domain;
706     if (domain != null)
707     {
708       String coding = "";
709       if ("Exon".equalsIgnoreCase(property)
710               || "Coding".equalsIgnoreCase(property))
711       {
712         coding = " Coding";
713       }
714       else if ("Intron".equalsIgnoreCase(property)
715               || "Noncoding".equalsIgnoreCase(property))
716       {
717         coding = " Noncoding";
718       }
719       verboseDomain = domain
720               + (currentGene == null ? "" : " (" + currentGene + coding
721                       + ")");
722     }
723     return verboseDomain;
724   }
725
726   /**
727    * Start processing a new feature
728    * 
729    * @param startPositions
730    */
731   protected void startSequenceFeature(Map<String, Integer> startPositions)
732   {
733     /*
734      * If the feature declaration precedes all sequences, we will know in
735      * createFeature that it started with residue 1; otherwise note now where it
736      * starts in each sequence
737      */
738     if (!residuesRead.isEmpty())
739     {
740       for (Entry<String, Integer> entry : residuesRead.entrySet())
741       {
742         String seqId = entry.getKey();
743         Integer nextResidue = entry.getValue() + 1;
744         startPositions.put(seqId, nextResidue);
745       }
746     }
747   }
748
749   /**
750    * Add a SequenceFeature to each sequence, using the given start/end values
751    * per sequence
752    * 
753    * @param featureType
754    * @param featureValue
755    * @param featureStartResidues
756    */
757   protected void createFeature(String featureType, String featureValue,
758           Map<String, Integer> featureStartResidues)
759   {
760     if (featureValue == null)
761     {
762       return;
763     }
764
765     Iterator<String> seqids = this.seqData.keySet().iterator();
766     while (seqids.hasNext())
767     {
768       String seqid = seqids.next();
769       Integer startAt = featureStartResidues.get(seqid);
770       int sfstart = startAt == null ? 1 : startAt.intValue();
771       int sfend = residuesRead.get(seqid);
772       if (sfend >= sfstart)
773       {
774         /*
775          * don't add feature if entirely gapped in the sequence
776          */
777         // TODO: type="Gene" (but then all coloured the same) or
778         // type="GeneName"?
779         SequenceFeature sf = new SequenceFeature(featureValue, featureType,
780                 sfstart, sfend, 0f, null);
781         sequenceFeatures.get(seqid).add(sf);
782       }
783     }
784   }
785
786   /**
787    * Returns the next line that is not a comment, or null at end of file.
788    * Comments in MEGA are within [ ] brackets, and may be nested.
789    * 
790    * @return
791    * @throws FileFormatException
792    */
793   protected String nextNonCommentLine() throws FileFormatException
794   {
795     return nextNonCommentLine(0);
796   }
797
798   /**
799    * Returns the next non-comment line (or part line), or null at end of file.
800    * Comments in MEGA are within [ ] brackets, and may be nested. They may occur
801    * anywhere within a line (for example at the end with position numbers); this
802    * method returns the line with any comments removed.
803    * 
804    * @param depth
805    *          current depth of nesting of comments while parsing
806    * @return
807    * @throws FileFormatException
808    */
809   protected String nextNonCommentLine(final int depth)
810           throws FileFormatException
811   {
812     String data = null;
813     try
814     {
815       data = nextLine();
816     } catch (IOException e)
817     {
818       throw new FileFormatException(e.getMessage());
819     }
820     if (data == null)
821     {
822       if (depth > 0)
823       {
824         System.err.println("Warning: unterminated comment in data file");
825       }
826       return data;
827     }
828
829     /*
830      * If we are in a (possibly nested) comment after parsing this line, keep
831      * reading recursively until the comment has unwound
832      */
833     int newDepth = commentDepth(data, depth);
834     if (newDepth > 0)
835     {
836       return nextNonCommentLine(newDepth);
837     }
838     else
839     {
840       /*
841        * not in a comment by end of this line; return what is left
842        */
843       String nonCommentPart = getNonCommentContent(data, depth);
844       return nonCommentPart;
845     }
846   }
847
848   /**
849    * Returns what is left of the input data after removing any comments, whether
850    * 'in progress' from preceding lines, or embedded in the current line
851    * 
852    * @param data
853    *          input data
854    * @param depth
855    *          nested depth of comments pending termination
856    * @return
857    * @throws FileFormatException
858    */
859   protected static String getNonCommentContent(String data, int depth)
860           throws FileFormatException
861   {
862     int len = data.length();
863     StringBuilder result = new StringBuilder(len);
864     for (int i = 0; i < len; i++)
865     {
866       char c = data.charAt(i);
867       switch (c)
868       {
869       case COMMENT_START:
870         depth++;
871         break;
872
873       case COMMENT_END:
874         if (depth > 0)
875         {
876           depth--;
877         }
878         else
879         {
880           result.append(c);
881         }
882         break;
883
884       default:
885         if (depth == 0)
886         {
887           result.append(c);
888         }
889       }
890     }
891     return result.toString();
892   }
893
894   /**
895    * Calculates new depth of comment after parsing an input line i.e. the excess
896    * of opening '[' over closing ']' characters. Any excess ']' are ignored (not
897    * treated as comment delimiters).
898    * 
899    * @param data
900    *          input line
901    * @param depth
902    *          current comment nested depth before parsing the line
903    * @return new depth after parsing the line
904    */
905   protected static int commentDepth(CharSequence data, int depth)
906   {
907     int newDepth = depth;
908     int len = data.length();
909     for (int i = 0; i < len; i++)
910     {
911       char c = data.charAt(i);
912       if (c == COMMENT_START)
913       {
914         newDepth++;
915       }
916       else if (c == COMMENT_END && newDepth > 0)
917       {
918         newDepth--;
919       }
920     }
921     return newDepth;
922   }
923
924   /**
925    * Convert the parsed sequence strings to objects and store them in the model.
926    */
927   protected void deriveSequencesAndFeatures()
928   {
929     Set<Entry<String, StringBuilder>> datasets = seqData.entrySet();
930
931     for (Entry<String, StringBuilder> dataset : datasets)
932     {
933       String sequenceId = dataset.getKey();
934       StringBuilder characters = dataset.getValue();
935       SequenceI s = new Sequence(sequenceId, new String(characters));
936       this.seqs.addElement(s);
937
938       /*
939        * and add any derived sequence features to the sequence
940        */
941       for (SequenceFeature sf : sequenceFeatures.get(sequenceId))
942       {
943         s.addSequenceFeature(sf);
944       }
945     }
946   }
947
948   /**
949    * Process one line of sequence data. If it has no sequence identifier, append
950    * to the current id's sequence. Else parse out the sequence id and append the
951    * data (if any) to that id's sequence. Returns the sequence id (implicit or
952    * explicit) for this line.
953    * 
954    * @param dataLine
955    * @return
956    * @throws FileFormatException
957    */
958   protected String parseDataLine(String dataLine)
959           throws FileFormatException
960   {
961     String seqId = getSequenceId(dataLine);
962     if (seqId == null)
963     {
964       /*
965        * Just character data
966        */
967       parseNoninterleavedDataLine(dataLine);
968       return currentSequenceId;
969     }
970     else if ((HASHSIGN + seqId).trim().equals(dataLine.trim()))
971     {
972       /*
973        * Sequence id only - header line for noninterleaved data
974        */
975       return seqId;
976     }
977     else
978     {
979       /*
980        * Sequence id followed by data
981        */
982       parseInterleavedDataLine(dataLine, seqId);
983       return seqId;
984     }
985   }
986
987   /**
988    * Add a line of sequence data to the buffer for the given sequence id. Start
989    * a new one if we haven't seen it before.
990    * 
991    * @param dataLine
992    * @throws FileFormatException
993    */
994   protected void parseNoninterleavedDataLine(String dataLine)
995           throws FileFormatException
996   {
997     if (currentSequenceId == null)
998     {
999       /*
1000        * Oops. Data but no sequence id context.
1001        */
1002       throw new FileFormatException("No sequence id context at: "
1003               + dataLine);
1004     }
1005
1006     assertInterleaved(false, dataLine);
1007
1008     dataLine = addSequenceData(currentSequenceId, dataLine);
1009
1010     setPositionsPerLine(Math.max(positionsPerLine, dataLine.length()));
1011   }
1012
1013   /**
1014    * Get the sequence data for this sequence id, starting a new one if
1015    * necessary.
1016    * 
1017    * @param currentId
1018    * @return
1019    */
1020   protected StringBuilder getSequenceDataBuffer(String currentId)
1021   {
1022     StringBuilder sb = seqData.get(currentId);
1023     if (sb == null)
1024     {
1025       // first data met for this sequence id, start a new buffer
1026       sb = new StringBuilder(SEQBUFFERSIZE);
1027       seqData.put(currentId, sb);
1028
1029       // and a placeholder for any SequenceFeature found
1030       sequenceFeatures.put(currentId, new ArrayList<SequenceFeature>());
1031     }
1032     return sb;
1033   }
1034
1035   /**
1036    * Parse one line of interleaved data e.g.
1037    * 
1038    * <pre>
1039    * #TheSeqId CGATCGCATGCA
1040    * </pre>
1041    * 
1042    * @param dataLine
1043    * @param seqId
1044    * @throws FileFormatException
1045    */
1046   protected void parseInterleavedDataLine(String dataLine, String seqId)
1047           throws FileFormatException
1048   {
1049     /*
1050      * New sequence found in second or later data block - error.
1051      */
1052     if (this.firstDataBlockRead && !seqData.containsKey(seqId))
1053     {
1054       throw new FileFormatException(
1055               "Parse error: misplaced new sequence starting at " + dataLine);
1056     }
1057
1058     String data = dataLine.substring(seqId.length() + 1).trim();
1059
1060     /*
1061      * Do nothing if this line is _only_ a sequence id with no data following.
1062      */
1063     if (data != null && data.length() > 0)
1064     {
1065       data = addSequenceData(seqId, data);
1066       setPositionsPerLine(Math.max(positionsPerLine, data.length()));
1067       assertInterleaved(true, dataLine);
1068     }
1069   }
1070
1071   /**
1072    * Remove spaces, and replace identity symbol, before appending the sequence
1073    * data to the buffer for the sequence id. Returns the reformatted added data.
1074    * Also updates a count of residues read for the sequence.
1075    * 
1076    * @param seqId
1077    * @param data
1078    * @return
1079    */
1080   protected String addSequenceData(String seqId, String data)
1081   {
1082     StringBuilder sb = getSequenceDataBuffer(seqId);
1083     int len = sb.length();
1084     String formatted = data.replace(SPACE, "");
1085
1086     /*
1087      * If sequence contains '.' or other identity symbol; replace these with the
1088      * same position from the first (reference) sequence
1089      */
1090     int nonGapped = 0;
1091     StringBuilder referenceSequence = seqData.values().iterator().next();
1092     StringBuilder sb1 = new StringBuilder(formatted.length());
1093     for (int i = 0; i < formatted.length(); i++)
1094     {
1095       char nextChar = formatted.charAt(i);
1096       if (nextChar == gapCharacter)
1097       {
1098         sb1.append(Comparison.isGap(nextChar) ? nextChar : DEFAULT_GAP);
1099       }
1100       else
1101       {
1102         nonGapped++;
1103         if (nextChar == identityCharacter
1104                 && len + i < referenceSequence.length())
1105         {
1106           sb1.append(referenceSequence.charAt(len + i));
1107         }
1108         else
1109         {
1110           sb1.append(nextChar);
1111         }
1112       }
1113     }
1114     formatted = sb1.toString();
1115
1116     data = formatted;
1117     sb.append(data);
1118
1119     /*
1120      * increment residue count for the sequence
1121      */
1122     if (nonGapped > 0)
1123     {
1124       Integer residueCount = residuesRead.get(seqId);
1125       residuesRead.put(seqId, nonGapped
1126               + (residueCount == null ? 0 : residueCount));
1127     }
1128
1129     return data;
1130   }
1131
1132   /**
1133    * If the line begins with (e.g.) "#abcde " then returns "abcde" as the
1134    * identifier. Else returns null.
1135    * 
1136    * @param dataLine
1137    * @return
1138    */
1139   public static String getSequenceId(String dataLine)
1140   {
1141     // TODO refactor to a StringUtils type class
1142     if (dataLine != null)
1143     {
1144       if (dataLine.startsWith(HASHSIGN))
1145       {
1146         int spacePos = dataLine.indexOf(" ");
1147         return (spacePos == -1 ? dataLine.substring(1) : dataLine
1148                 .substring(1, spacePos));
1149       }
1150     }
1151     return null;
1152   }
1153
1154   /**
1155    * Read the #MEGA and Title/Format/Description header lines (if present).
1156    * 
1157    * Save as alignment properties in case useful.
1158    * 
1159    * @return the next non-blank line following the header lines.
1160    * @throws FileFormatException
1161    */
1162   protected String parseHeaderLines() throws FileFormatException
1163   {
1164     String inputLine = null;
1165     while ((inputLine = nextNonCommentLine()) != null)
1166     {
1167       inputLine = inputLine.trim();
1168
1169       /*
1170        * skip blank lines
1171        */
1172       if (inputLine.length() == 0)
1173       {
1174         continue;
1175       }
1176
1177       if (inputLine.toUpperCase().startsWith(MEGA_ID))
1178       {
1179         continue;
1180       }
1181
1182       if (isTitle(inputLine))
1183       {
1184         this.title = getValue(inputLine);
1185         setAlignmentProperty(PROP_TITLE, title);
1186       }
1187       else if (inputLine.startsWith(BANG + DESCRIPTION))
1188       {
1189         parseDescription(inputLine);
1190       }
1191
1192       else if (inputLine.startsWith(BANG + FORMAT))
1193       {
1194         parseFormat(inputLine);
1195       }
1196       else if (!inputLine.toUpperCase().startsWith(MEGA_ID))
1197       {
1198
1199         /*
1200          * Return the first 'data line' i.e. one that is not blank, #MEGA or
1201          * TITLE:
1202          */
1203         break;
1204       }
1205     }
1206     return inputLine;
1207   }
1208
1209   /**
1210    * Parse a !Format statement. This may be multiline, and is ended by a
1211    * semicolon.
1212    * 
1213    * @param inputLine
1214    * @throws FileFormatException
1215    */
1216   protected void parseFormat(String inputLine) throws FileFormatException
1217   {
1218     while (inputLine != null)
1219     {
1220       parseFormatLine(inputLine);
1221       if (inputLine.endsWith(SEMICOLON))
1222       {
1223         break;
1224       }
1225       inputLine = nextNonCommentLine();
1226     }
1227   }
1228
1229   /**
1230    * Parse one line of a !Format statement. This may contain one or more
1231    * keyword=value pairs.
1232    * 
1233    * @param inputLine
1234    * @throws FileFormatException
1235    */
1236   protected void parseFormatLine(String inputLine)
1237           throws FileFormatException
1238   {
1239     if (inputLine.startsWith(BANG + FORMAT))
1240     {
1241       inputLine = inputLine.substring((BANG + FORMAT).length());
1242     }
1243     if (inputLine.endsWith(SEMICOLON))
1244     {
1245       inputLine = inputLine.substring(0, inputLine.length() - 1);
1246     }
1247     if (inputLine.length() == 0)
1248     {
1249       return;
1250     }
1251     String[] tokens = inputLine.trim().split(WHITESPACE);
1252     for (String token : tokens)
1253     {
1254       parseFormatKeyword(token);
1255     }
1256   }
1257
1258   /**
1259    * Parse a Keyword=Value token. Possible keywords are
1260    * <ul>
1261    * <li>DataType= DNA, RNA, Nucleotide, Protein</li>
1262    * <li>DataFormat= Interleaved, ?</li>
1263    * <li>NSeqs= number of sequences (synonym NTaxa)</li>
1264    * <li>NSites= number of bases / residues</li>
1265    * <li>Property= Exon (or Coding), Intron (or Noncoding), End (of domain)</li>
1266    * <li>Indel= gap character</li>
1267    * <li>Identical= identity character (to first sequence) (synonym MatchChar)</li>
1268    * <li>Missing= missing data character</li>
1269    * <li>CodeTable= Standard, other (MEGA supports various)</li>
1270    * </ul>
1271    * 
1272    * @param token
1273    * @throws FileFormatException
1274    *           if an unrecognised keyword or value is encountered
1275    */
1276   protected void parseFormatKeyword(String token)
1277           throws FileFormatException
1278   {
1279     String msg = "Unrecognised Format command: " + token;
1280     String[] bits = token.split(EQUALS);
1281     if (bits.length != 2)
1282     {
1283       throw new FileFormatException(msg);
1284     }
1285     String keyword = bits[0];
1286     String value = bits[1];
1287
1288     /*
1289      * Jalview will work out whether nucleotide or not anyway
1290      */
1291     if (keyword.equalsIgnoreCase(DATATYPE))
1292     {
1293       if (value.equalsIgnoreCase("DNA") || value.equalsIgnoreCase("RNA")
1294               || value.equalsIgnoreCase("Nucleotide"))
1295       {
1296         this.nucleotide = true;
1297         // alignment computes whether or not it is nucleotide when created
1298       }
1299       else if (value.equalsIgnoreCase(PROTEIN))
1300       {
1301         this.nucleotide = false;
1302       }
1303       else
1304       {
1305         throw new FileFormatException(msg);
1306       }
1307       setAlignmentProperty(PROP_DATATYPE, value);
1308     }
1309
1310     /*
1311      * accept non-Standard code table but save in case we want to disable
1312      * 'translate as cDNA'
1313      */
1314     else if (keyword.equalsIgnoreCase(CODETABLE))
1315     {
1316       setAlignmentProperty(PROP_CODETABLE, value);
1317     }
1318
1319     /*
1320      * save gap char to set later on alignment once created
1321      */
1322     else if (keyword.equalsIgnoreCase(INDEL))
1323     {
1324       this.gapCharacter = value.charAt(0);
1325       if (!Comparison.isGap(gapCharacter))
1326       {
1327         System.err.println("Jalview doesn't support '" + gapCharacter
1328                 + "' for gaps, will be converted to '" + DEFAULT_GAP + "'");
1329       }
1330     }
1331
1332     else if (keyword.equalsIgnoreCase(IDENTICAL)
1333             || keyword.equalsIgnoreCase("MatchChar"))
1334     {
1335       setAlignmentProperty(PROP_IDENTITY, value);
1336       this.identityCharacter = value.charAt(0);
1337       if (!".".equals(value))
1338       {
1339         System.err.println("Warning: " + token
1340                 + " not supported, Jalview uses '.' for identity");
1341       }
1342     }
1343
1344     else if (keyword.equalsIgnoreCase(MISSING))
1345     {
1346       setAlignmentProperty(PROP_MISSING, value);
1347       System.err.println("Warning: " + token + " not supported");
1348     }
1349
1350     else if (keyword.equalsIgnoreCase(PROPERTY))
1351     {
1352       // TODO: can Property appear in a Format command?
1353       // suspect this is a mistake in the manual
1354     }
1355
1356     else if (!keyword.equalsIgnoreCase(N_SEQS)
1357             && !keyword.equalsIgnoreCase("NTaxa")
1358             && !keyword.equalsIgnoreCase(N_SITES))
1359     {
1360       System.err.println("Warning: " + msg);
1361     }
1362   }
1363
1364   /**
1365    * Returns the trimmed data on the line following either whitespace or '=',
1366    * with any trailing semi-colon removed<br>
1367    * So
1368    * <ul>
1369    * <li>Hello World</li>
1370    * <li>!Hello: \tWorld;</li>
1371    * <li>!Hello=World</li>
1372    * <ul>
1373    * should all return "World"
1374    * 
1375    * @param inputLine
1376    * @return
1377    */
1378   protected static String getValue(String inputLine)
1379   {
1380     if (inputLine == null)
1381     {
1382       return null;
1383     }
1384     String value = null;
1385     String s = inputLine.replaceAll("\t", " ").trim();
1386
1387     /*
1388      * KEYWORD = VALUE should return VALUE
1389      */
1390     int equalsPos = s.indexOf("=");
1391     if (equalsPos >= 0)
1392     {
1393       value = s.substring(equalsPos + 1);
1394     }
1395     else
1396     {
1397       int spacePos = s.indexOf(' ');
1398       value = spacePos == -1 ? "" : s.substring(spacePos + 1);
1399     }
1400     value = value.trim();
1401     if (value.endsWith(SEMICOLON))
1402     {
1403       value = value.substring(0, value.length() - 1).trim();
1404     }
1405     return value;
1406   }
1407
1408   /**
1409    * Returns true if the input line starts with "TITLE" or "!TITLE" (not case
1410    * sensitive). The latter is the official format, some older data file
1411    * examples have it without the !.
1412    * 
1413    * @param inputLine
1414    * @return
1415    */
1416   protected static boolean isTitle(String inputLine)
1417   {
1418     if (inputLine == null)
1419     {
1420       return false;
1421     }
1422     String upper = inputLine.toUpperCase();
1423     return (upper.startsWith(TITLE.toUpperCase()) || upper.startsWith(BANG
1424             + TITLE.toUpperCase()));
1425   }
1426
1427   /**
1428    * Reads lines until terminated by semicolon, appending each to the
1429    * Description property value.
1430    * 
1431    * @throws FileFormatException
1432    */
1433   protected void parseDescription(String firstDescriptionLine)
1434           throws FileFormatException
1435   {
1436     StringBuilder desc = new StringBuilder(256);
1437     desc.append(getValue(firstDescriptionLine));
1438     if (!firstDescriptionLine.endsWith(SEMICOLON))
1439     {
1440       String line = nextNonCommentLine();
1441       while (line != null)
1442       {
1443         if (line.endsWith(SEMICOLON))
1444         {
1445           desc.append(line.substring(0, line.length() - 1));
1446           break;
1447         }
1448         else if (line.length() > 0)
1449         {
1450           desc.append(line).append(newline);
1451         }
1452         line = nextNonCommentLine();
1453       }
1454     }
1455     setAlignmentProperty(PROP_DESCRIPTION, desc.toString());
1456   }
1457
1458   /**
1459    * Returns the alignment sequences in Mega format.
1460    */
1461   @Override
1462   public String print()
1463   {
1464     return MEGA_ID + newline + print(getSeqsAsArray());
1465   }
1466
1467   /**
1468    * Write out the alignment sequences in Mega format - interleaved unless
1469    * explicitly noninterleaved.
1470    */
1471   protected String print(SequenceI[] s)
1472   {
1473     String result;
1474     if (this.interleaved != null && !this.interleaved)
1475     {
1476       result = printNonInterleaved(s);
1477     }
1478     else
1479     {
1480       result = printInterleaved(s);
1481     }
1482     return result;
1483   }
1484
1485   /**
1486    * Print to string in Interleaved format - blocks of next N characters of each
1487    * sequence in turn.
1488    * 
1489    * @param s
1490    */
1491   protected String printInterleaved(SequenceI[] s)
1492   {
1493     int maxIdLength = getMaxIdLength(s);
1494     int maxSequenceLength = getMaxSequenceLength(s);
1495     int numLines = maxSequenceLength / positionsPerLine + 3; // approx
1496
1497     int numDataBlocks = (maxSequenceLength - 1) / positionsPerLine + 1;
1498     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1499     int chunksPerLine = (positionsPerLine + spaceEvery - 1) / spaceEvery;
1500
1501     /*
1502      * Roughly size a buffer to hold the whole output
1503      */
1504     StringBuilder sb = new StringBuilder(numLines
1505             * (maxIdLength + positionsPerLine + chunksPerLine + 10));
1506
1507     /*
1508      * Output as: #Seqid CGT AGC ACT ... or blocks of 10 for peptide
1509      */
1510     int from = 0;
1511     for (int i = 0; i < numDataBlocks; i++)
1512     {
1513       sb.append(newline);
1514       boolean first = true;
1515       int advancedBy = 0;
1516       for (SequenceI seq : s)
1517       {
1518         int seqFrom = from;
1519         String seqId = String.format("#%-" + maxIdLength + "s",
1520                 seq.getName());
1521
1522         /*
1523          * output next line for this sequence
1524          */
1525         sb.append(seqId);
1526         int lastPos = seqFrom + positionsPerLine; // exclusive
1527         for (int j = 0; j < chunksPerLine; j++)
1528         {
1529           char[] subSequence = seq.getSequence(seqFrom,
1530                   Math.min(lastPos, seqFrom + spaceEvery));
1531           if (subSequence.length > 0)
1532           {
1533             sb.append(SPACE).append(subSequence);
1534           }
1535           seqFrom += subSequence.length;
1536           if (first)
1537           {
1538             // all sequences should be the same length in MEGA
1539             advancedBy += subSequence.length;
1540           }
1541         }
1542         // write last position as a comment
1543         if (writePositionNumbers)
1544         {
1545           sb.append(SPACE).append(COMMENT_START).append(from + advancedBy)
1546                   .append(COMMENT_END);
1547         }
1548         sb.append(newline);
1549         first = false;
1550       }
1551       sb.append(printLabel(from, advancedBy, maxIdLength));
1552       from += advancedBy;
1553     }
1554
1555     return new String(sb);
1556   }
1557
1558   /**
1559    * Outputs to string the MEGA header and any other known and relevant
1560    * alignment properties
1561    * 
1562    * @param al
1563    */
1564   protected String printHeaders(AlignmentI al)
1565   {
1566     StringBuilder sb = new StringBuilder(128);
1567     sb.append(MEGA_ID).append(newline);
1568     String propertyValue = (String) al.getProperty(PROP_TITLE);
1569     if (propertyValue != null)
1570     {
1571       sb.append(BANG).append(TITLE).append(SPACE).append(propertyValue)
1572               .append(SEMICOLON).append(newline);
1573     }
1574     propertyValue = (String) al.getProperty(PROP_DESCRIPTION);
1575     if (propertyValue != null)
1576     {
1577       sb.append(BANG).append(DESCRIPTION).append(newline)
1578               .append(propertyValue).append(SEMICOLON)
1579               .append(newline);
1580     }
1581
1582     /*
1583      * !Format DataType CodeTable
1584      */
1585     sb.append(BANG).append(FORMAT).append(newline);
1586     String dataType = (String) al.getProperty(PROP_DATATYPE);
1587     if (dataType == null)
1588     {
1589       dataType = al.isNucleotide() ? NUCLEOTIDE : PROTEIN;
1590     }
1591     sb.append(INDENT).append(DATATYPE).append(EQUALS).append(dataType);
1592     String codeTable = (String) al.getProperty(PROP_CODETABLE);
1593     sb.append(SPACE).append(CODETABLE).append(EQUALS)
1594             .append(codeTable == null ? "Standard" : codeTable)
1595             .append(newline);
1596     
1597     /*
1598      * !Format NSeqs NSites (the length of sequences - they should all be the
1599      * same - including gaps)
1600      */
1601     sb.append(INDENT).append(N_SEQS).append(EQUALS).append(al.getHeight());
1602     sb.append(SPACE).append(N_SITES).append(EQUALS)
1603             .append(String.valueOf(al.getWidth()));
1604     sb.append(newline);
1605
1606     /*
1607      * !Format Indel Identical Missing
1608      */
1609     sb.append(INDENT);
1610     sb.append(INDEL).append(EQUALS).append(al.getGapCharacter());
1611     String identity = (String) al.getProperty(PROP_IDENTITY);
1612     if (identity != null)
1613     {
1614       sb.append(SPACE).append(IDENTICAL).append(EQUALS).append(identity);
1615     }
1616     String missing = (String) al.getProperty(PROP_MISSING);
1617     if (missing != null)
1618     {
1619       sb.append(SPACE).append(MISSING).append(EQUALS).append(missing);
1620     }
1621     sb.append(SEMICOLON).append(newline);
1622
1623     return sb.toString();
1624   }
1625
1626   /**
1627    * Get the longest sequence id (to allow aligned printout).
1628    * 
1629    * @param s
1630    * @return
1631    */
1632   protected static int getMaxIdLength(SequenceI[] s)
1633   {
1634     // TODO pull up for reuse
1635     int maxLength = 0;
1636     for (SequenceI seq : s)
1637     {
1638       int len = seq.getName().length();
1639       if (len > maxLength)
1640       {
1641         maxLength = len;
1642       }
1643     }
1644     return maxLength;
1645   }
1646
1647   /**
1648    * Get the longest sequence length
1649    * 
1650    * @param s
1651    * @return
1652    */
1653   protected static int getMaxSequenceLength(SequenceI[] s)
1654   {
1655     // TODO pull up for reuse
1656     int maxLength = 0;
1657     for (SequenceI seq : s)
1658     {
1659       int len = seq.getLength();
1660       if (len > maxLength)
1661       {
1662         maxLength = len;
1663       }
1664     }
1665     return maxLength;
1666   }
1667
1668   /**
1669    * Print to string in noninterleaved format - all of each sequence in turn, in
1670    * blocks of 50 characters.
1671    * 
1672    * @param s
1673    * @return
1674    */
1675   protected String printNonInterleaved(SequenceI[] s)
1676   {
1677     int maxSequenceLength = getMaxSequenceLength(s);
1678     // approx
1679     int numLines = maxSequenceLength / positionsPerLine + 2 + s.length;
1680
1681     /*
1682      * Roughly size a buffer to hold the whole output
1683      */
1684     StringBuilder sb = new StringBuilder(numLines * positionsPerLine);
1685
1686     for (SequenceI seq : s)
1687     {
1688       printSequence(sb, seq);
1689     }
1690
1691     return new String(sb);
1692   }
1693
1694   /**
1695    * Append a formatted complete sequence to the string buffer
1696    * 
1697    * @param sb
1698    * @param seq
1699    */
1700   protected void printSequence(StringBuilder sb, SequenceI seq)
1701   {
1702     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1703     // round down to output a whole number of spaced blocks
1704     int chunksPerLine = positionsPerLine / spaceEvery;
1705
1706     sb.append(newline);
1707     sb.append(HASHSIGN + seq.getName()).append(newline);
1708     int startPos = 0;
1709     while (startPos < seq.getLength())
1710     {
1711       /*
1712        * print next line for this sequence
1713        */
1714       boolean firstChunk = true;
1715       int lastPos = startPos + positionsPerLine; // exclusive
1716       for (int j = 0; j < chunksPerLine; j++)
1717       {
1718         char[] subSequence = seq.getSequence(startPos,
1719                 Math.min(lastPos, startPos + spaceEvery));
1720         if (subSequence.length > 0)
1721         {
1722           if (!firstChunk)
1723           {
1724             sb.append(SPACE);
1725           }
1726           sb.append(subSequence);
1727           firstChunk = false;
1728         }
1729         startPos += subSequence.length;
1730       }
1731       // line end position (base 1) as a comment
1732       sb.append(SPACE).append(COMMENT_START).append(startPos)
1733               .append(COMMENT_END);
1734       sb.append(newline);
1735     }
1736   }
1737
1738   /**
1739    * Returns a formatted string like <br>
1740    * !Label aa_b_ ab_b_ <br>
1741    * where underscore represents no annotation, any other character a MEGA label
1742    * character <br>
1743    * Returns an empty string if there is no non-null annotation in the given
1744    * alignment range
1745    * 
1746    * @param fromPos
1747    *          start column of the alignment (base 0)
1748    * @param positions
1749    *          number of positions to output
1750    * @param labelWidth
1751    *          padded width of !Label statement to output
1752    * @return
1753    */
1754   protected String printLabel(int fromPos, int positions, int labelWidth)
1755   {
1756     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1757     String none = "";
1758     AlignmentAnnotation ann = findAnnotation(MEGA_ANNOTATION_LABEL);
1759     if (ann == null)
1760     {
1761       return none;
1762     }
1763
1764     StringBuilder sb = new StringBuilder(positions + 20);
1765     sb.append(String.format("%-" + labelWidth + "s ", BANG + LABEL));
1766     Annotation[] anns = annotations.get(0).annotations;
1767     int blockCharCount = 0;
1768     boolean annotationFound = false;
1769
1770     for (int i = fromPos; i < fromPos + positions; i++)
1771     {
1772       String label = String.valueOf(UNDERSCORE);
1773       if (i < anns.length && anns[i] != null)
1774       {
1775         label = anns[i].displayCharacter;
1776       }
1777       sb.append(label);
1778       if (label.charAt(0) != UNDERSCORE)
1779       {
1780         annotationFound = true;
1781       }
1782       // add a space after each block except the last
1783       if (++blockCharCount % spaceEvery == 0
1784               && (i < fromPos + positions - 1))
1785       {
1786         sb.append(SPACE);
1787       }
1788     }
1789     sb.append(SEMICOLON).append(newline);
1790
1791     return annotationFound ? sb.toString() : none;
1792   }
1793
1794   /**
1795    * Returns the first stored annotation found with the given label, or null
1796    * 
1797    * @param annotationLabel
1798    * @return
1799    */
1800   protected AlignmentAnnotation findAnnotation(String annotationLabel)
1801   {
1802     if (annotations != null)
1803     {
1804       for (AlignmentAnnotation ann : annotations)
1805       {
1806         if (annotationLabel.equals(ann.label))
1807         {
1808           return ann;
1809         }
1810       }
1811     }
1812     return null;
1813   }
1814
1815   /**
1816    * Flag this file as interleaved or not, based on data format. Throws an
1817    * exception if has previously been determined to be otherwise.
1818    * 
1819    * @param isIt
1820    * @param dataLine
1821    * @throws FileFormatException
1822    */
1823   protected void assertInterleaved(boolean isIt, String dataLine)
1824           throws FileFormatException
1825   {
1826     if (this.interleaved != null && isIt != this.interleaved.booleanValue())
1827     {
1828       throw new FileFormatException("Parse error: interleaved was " + !isIt
1829               + " but now seems to be " + isIt + ", at line: " + dataLine);
1830     }
1831     this.interleaved = new Boolean(isIt);
1832     setAlignmentProperty(PROP_INTERLEAVED, interleaved.toString());
1833   }
1834
1835   public boolean isInterleaved()
1836   {
1837     return this.interleaved == null ? false : this.interleaved
1838             .booleanValue();
1839   }
1840
1841   /**
1842    * Adds saved parsed values either as alignment properties, or (in some cases)
1843    * as specific member fields of the alignment
1844    */
1845   @Override
1846   public void addProperties(AlignmentI al)
1847   {
1848     super.addProperties(al);
1849     
1850     /*
1851      * record gap character specified, but convert to '-' if not one we support
1852      */
1853     al.setGapCharacter(Comparison.isGap(gapCharacter) ? gapCharacter
1854             : DEFAULT_GAP);
1855
1856     /*
1857      * warn if e.g. DataType=DNA but data is protein (or vice versa)
1858      */
1859     if (this.nucleotide != null && this.nucleotide != al.isNucleotide()) {
1860       System.err.println("Warning: " + this.title + " declared "
1861               + (nucleotide ? "" : " not ") + "nucleotide but it is"
1862               + (nucleotide ? " not" : ""));
1863     }
1864   }
1865
1866   /**
1867    * Print the given alignment in MEGA format. If the alignment was created by
1868    * parsing a MEGA file, it should have properties set (e.g. Title) which can
1869    * surface in the output.
1870    */
1871   @Override
1872   public String print(AlignmentI al)
1873   {
1874     this.nucleotide = al.isNucleotide();
1875
1876     /*
1877      * if the alignment has a "MEGA" annotation, we'll output its values as
1878      * !Label statements; MEGA only supports one of these
1879      */
1880     AlignmentAnnotation[] anns = al.getAlignmentAnnotation();
1881     if (anns != null)
1882     {
1883       for (AlignmentAnnotation ann : anns)
1884       {
1885         if (MEGA_ANNOTATION_LABEL.equals(ann.label))
1886         {
1887           this.annotations = new Vector<AlignmentAnnotation>();
1888           annotations.add(ann);
1889           break;
1890         }
1891       }
1892     }
1893
1894     String lineLength = (String) al.getProperty(PROP_LINELENGTH);
1895     this.positionsPerLine = lineLength == null ? DEFAULT_LINE_LENGTH : Integer
1896             .parseInt(lineLength);
1897
1898     /*
1899      * round down to a multiple of 3 positions per line for nucleotide
1900      */
1901     if (nucleotide)
1902     {
1903       positionsPerLine = positionsPerLine - (positionsPerLine % 3);
1904     }
1905
1906     String interleave = (String) al.getProperty(PROP_INTERLEAVED);
1907     if (interleave != null)
1908     {
1909       this.interleaved = Boolean.valueOf(interleave);
1910     }
1911
1912     String headers = printHeaders(al);
1913     return headers + print(al.getSequencesArray());
1914   }
1915
1916   /**
1917    * Returns the number of sequence positions output per line
1918    * 
1919    * @return
1920    */
1921   public int getPositionsPerLine()
1922   {
1923     return positionsPerLine;
1924   }
1925
1926   /**
1927    * Sets the number of sequence positions output per line. Note these will be
1928    * formatted in blocks of 3 (nucleotide) or 10 (peptide).
1929    * 
1930    * @param p
1931    */
1932   public void setPositionsPerLine(int p)
1933   {
1934     this.positionsPerLine = p;
1935   }
1936 }