JAL-1499 output !Domain and !Gene statements from alignment annotations
[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         // NB use 'colour feature by label' to show up distinct instances of
778         // feature type 'Gene' or 'Domain' on the alignment
779         SequenceFeature sf = new SequenceFeature(featureType, featureValue,
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
1496     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1497     int chunksPerLine = (positionsPerLine + spaceEvery - 1) / spaceEvery;
1498
1499     /*
1500      * Roughly size a buffer to hold the whole output
1501      */
1502     int numLines = maxSequenceLength / positionsPerLine + 3; // approx
1503     StringBuilder sb = new StringBuilder(numLines
1504             * (maxIdLength + positionsPerLine + chunksPerLine + 10));
1505
1506     /*
1507      * Output as: #Seqid CGT AGC ACT ... or blocks of 10 for peptide
1508      */
1509     AlignmentAnnotation geneAnnotation = findAnnotation(MEGA_ANNOTATION_GENE);
1510     AlignmentAnnotation domainAnnotation = findAnnotation(MEGA_ANNOTATION_DOMAIN);
1511     currentGene = null;
1512     currentDomain = null;
1513     int from = 0;
1514
1515     while (from < maxSequenceLength)
1516     {
1517       printGeneOrDomainHeader(sb, from, geneAnnotation, domainAnnotation);
1518       int maxCol = from + positionsPerLine; // exclusive
1519       for (int col = from; col <= maxCol; col++)
1520       {
1521         if (geneOrDomainChanged(col, geneAnnotation, domainAnnotation)
1522                 || col == maxCol)
1523         {
1524           /*
1525            * print a block of sequences up to [col-1]
1526            */
1527           sb.append(newline);
1528           boolean first = true;
1529           int advancedBy = 0;
1530           for (SequenceI seq : s)
1531           {
1532             int seqFrom = from;
1533             String seqId = String.format("#%-" + maxIdLength + "s",
1534                     seq.getName());
1535
1536             /*
1537              * output next line for this sequence
1538              */
1539             sb.append(seqId);
1540             for (int j = 0; j < chunksPerLine; j++)
1541             {
1542               char[] subSequence = seq.getSequence(seqFrom,
1543                       Math.min(col, seqFrom + spaceEvery));
1544               if (subSequence.length > 0)
1545               {
1546                 sb.append(SPACE).append(subSequence);
1547               }
1548               seqFrom += subSequence.length;
1549               if (first)
1550               {
1551                 // all sequences should be the same length in MEGA
1552                 advancedBy += subSequence.length;
1553               }
1554             }
1555             // write last position as a comment
1556             if (writePositionNumbers)
1557             {
1558               sb.append(SPACE).append(COMMENT_START)
1559                       .append(String.valueOf(from + advancedBy))
1560                       .append(COMMENT_END);
1561             }
1562             sb.append(newline);
1563             first = false;
1564           }
1565           sb.append(printLabel(from, advancedBy, maxIdLength));
1566           from += advancedBy;
1567           break;
1568         }
1569       }
1570     }
1571
1572     return new String(sb);
1573   }
1574
1575   /**
1576    * Returns true if we detect a change of gene or domain at the given column
1577    * position. Currently done by inspecting any "MEGA Gene" or "MEGA Domain"
1578    * annotation for the column.
1579    * 
1580    * @param column
1581    * @param geneAnnotation
1582    * @param domainAnnotation
1583    * @return
1584    */
1585   protected boolean geneOrDomainChanged(int column,
1586           AlignmentAnnotation geneAnnotation,
1587           AlignmentAnnotation domainAnnotation)
1588   {
1589     String gene = getGeneFromAnnotation(column, geneAnnotation);
1590     String domain = getDomainFromAnnotation(column, domainAnnotation);
1591     boolean domainEnd = domain != null
1592             && domain.toUpperCase().indexOf(DOMAINEND.toUpperCase()) != -1;
1593
1594     if ((domainEnd || gene == null) && currentGene != null)
1595     {
1596       return true;
1597     }
1598     else if (gene != null && !gene.equals(currentGene))
1599     {
1600       return true;
1601     }
1602     // todo: Domain
1603     return false;
1604   }
1605
1606   /**
1607    * Extracts the name of a domain from MEGA Domain annotation at the given
1608    * column position, if any
1609    * 
1610    * @param column
1611    * @param domainAnnotation
1612    * @return
1613    */
1614   protected static String getDomainFromAnnotation(int column,
1615           AlignmentAnnotation domainAnnotation)
1616   {
1617     Annotation domainAnn = domainAnnotation == null ? null
1618             : (column >= domainAnnotation.annotations.length ? null
1619                     : domainAnnotation.annotations[column]);
1620     String domain = domainAnn == null ? null : domainAnn.description;
1621     if (domain != null && domain.indexOf("(") > 0)
1622     {
1623       domain = domain.substring(0, domain.indexOf("(")).trim();
1624     }
1625     return domain;
1626   }
1627
1628   /**
1629    * Extracts the name of a gene from MEGA Gene annotation at the given column
1630    * position, if any
1631    * 
1632    * @param column
1633    * @param geneAnnotation
1634    * @return
1635    */
1636   protected static String getGeneFromAnnotation(int column,
1637           AlignmentAnnotation geneAnnotation)
1638   {
1639     Annotation geneAnn = geneAnnotation == null ? null
1640             : (column >= geneAnnotation.annotations.length ? null
1641                     : geneAnnotation.annotations[column]);
1642     String gene = geneAnn == null ? null : geneAnn.description;
1643     return gene;
1644   }
1645
1646   /**
1647    * Appends a !Domain header on change of Domain (or Domain/Gene) annotation.
1648    * If it changes to a null value, appends a !Domain Property=domainend;
1649    * statement
1650    * 
1651    * @param sb
1652    *          buffer to append to
1653    * @param column
1654    *          of alignment being output
1655    * @param geneAnnotation
1656    *          "MEGA Gene" annotations
1657    * @param domainAnnotation
1658    *          "MEGA Domain" annotations
1659    */
1660   protected void printGeneOrDomainHeader(StringBuilder sb, int column,
1661           AlignmentAnnotation geneAnnotation,
1662           AlignmentAnnotation domainAnnotation)
1663   {
1664     String gene = getGeneFromAnnotation(column, geneAnnotation);
1665     String domain = getDomainFromAnnotation(column, domainAnnotation);
1666     String property = getPropertyFromAnnotation(column, domainAnnotation);
1667
1668     if ((gene == null && currentGene != null)
1669             || (gene == null && domain == null && currentDomain != null))
1670     {
1671       /*
1672        * end of Gene or Domain annotation
1673        */
1674       sb.append(newline).append(BANG).append(DOMAIN).append(EQUALS)
1675               .append(currentDomain).append(SPACE).append(PROPERTY)
1676               .append(EQUALS).append(DOMAINEND).append(SEMICOLON);
1677     }
1678     else if (gene != null && !gene.equals(currentGene))
1679     {
1680       /*
1681        * start of a new Gene; output as "!Domain Gene=..." if we have domain
1682        * annotation
1683        */
1684       if (domain != null)
1685       {
1686         sb.append(newline).append(BANG).append(DOMAIN).append(EQUALS)
1687                 .append(domain).append(SPACE).append(GENE).append(EQUALS)
1688                 .append(gene);
1689         if (property != null)
1690         {
1691           sb.append(SPACE).append(PROPERTY).append(EQUALS).append(property);
1692         }
1693         sb.append(SEMICOLON);
1694       }
1695       else
1696       {
1697         sb.append(newline).append(BANG).append(GENE).append(EQUALS)
1698               .append(gene).append(SEMICOLON);
1699       }
1700     }
1701     else if (domain != null && !domain.equals(currentGene))
1702     {
1703       /*
1704        * change of domain within same (or no) gene
1705        */
1706       sb.append(newline).append(BANG).append(DOMAIN).append(EQUALS)
1707               .append(domain);
1708       if (currentGene != null)
1709       {
1710         sb.append(SPACE).append(GENE).append(EQUALS).append(currentGene);
1711       }
1712       if (property != null)
1713       {
1714         sb.append(SPACE).append(PROPERTY).append(EQUALS).append(property);
1715       }
1716       sb.append(SEMICOLON);
1717     }
1718     currentGene = gene;
1719     currentDomain = domain;
1720   }
1721
1722   /**
1723    * Outputs to string the MEGA header and any other known and relevant
1724    * alignment properties
1725    * 
1726    * @param al
1727    */
1728   protected String printHeaders(AlignmentI al)
1729   {
1730     StringBuilder sb = new StringBuilder(128);
1731     sb.append(MEGA_ID).append(newline);
1732     String propertyValue = (String) al.getProperty(PROP_TITLE);
1733     if (propertyValue != null)
1734     {
1735       sb.append(BANG).append(TITLE).append(SPACE).append(propertyValue)
1736               .append(SEMICOLON).append(newline);
1737     }
1738     propertyValue = (String) al.getProperty(PROP_DESCRIPTION);
1739     if (propertyValue != null)
1740     {
1741       sb.append(BANG).append(DESCRIPTION).append(newline)
1742               .append(propertyValue).append(SEMICOLON)
1743               .append(newline);
1744     }
1745
1746     /*
1747      * !Format DataType CodeTable
1748      */
1749     sb.append(BANG).append(FORMAT).append(newline);
1750     String dataType = (String) al.getProperty(PROP_DATATYPE);
1751     if (dataType == null)
1752     {
1753       dataType = al.isNucleotide() ? NUCLEOTIDE : PROTEIN;
1754     }
1755     sb.append(INDENT).append(DATATYPE).append(EQUALS).append(dataType);
1756     String codeTable = (String) al.getProperty(PROP_CODETABLE);
1757     sb.append(SPACE).append(CODETABLE).append(EQUALS)
1758             .append(codeTable == null ? "Standard" : codeTable)
1759             .append(newline);
1760     
1761     /*
1762      * !Format NSeqs NSites (the length of sequences - they should all be the
1763      * same - including gaps)
1764      */
1765     sb.append(INDENT).append(N_SEQS).append(EQUALS).append(al.getHeight());
1766     sb.append(SPACE).append(N_SITES).append(EQUALS)
1767             .append(String.valueOf(al.getWidth()));
1768     sb.append(newline);
1769
1770     /*
1771      * !Format Indel Identical Missing
1772      */
1773     sb.append(INDENT);
1774     sb.append(INDEL).append(EQUALS).append(al.getGapCharacter());
1775     String identity = (String) al.getProperty(PROP_IDENTITY);
1776     if (identity != null)
1777     {
1778       sb.append(SPACE).append(IDENTICAL).append(EQUALS).append(identity);
1779     }
1780     String missing = (String) al.getProperty(PROP_MISSING);
1781     if (missing != null)
1782     {
1783       sb.append(SPACE).append(MISSING).append(EQUALS).append(missing);
1784     }
1785     sb.append(SEMICOLON).append(newline);
1786
1787     return sb.toString();
1788   }
1789
1790   /**
1791    * Get the longest sequence id (to allow aligned printout).
1792    * 
1793    * @param s
1794    * @return
1795    */
1796   protected static int getMaxIdLength(SequenceI[] s)
1797   {
1798     // TODO pull up for reuse
1799     int maxLength = 0;
1800     for (SequenceI seq : s)
1801     {
1802       int len = seq.getName().length();
1803       if (len > maxLength)
1804       {
1805         maxLength = len;
1806       }
1807     }
1808     return maxLength;
1809   }
1810
1811   /**
1812    * Get the longest sequence length (including gaps)
1813    * 
1814    * @param s
1815    * @return
1816    */
1817   protected static int getMaxSequenceLength(SequenceI[] s)
1818   {
1819     // TODO pull up for reuse
1820     int maxLength = 0;
1821     for (SequenceI seq : s)
1822     {
1823       int len = seq.getLength();
1824       if (len > maxLength)
1825       {
1826         maxLength = len;
1827       }
1828     }
1829     return maxLength;
1830   }
1831
1832   /**
1833    * Print to string in noninterleaved format - all of each sequence in turn, in
1834    * blocks of 50 characters.
1835    * 
1836    * @param s
1837    * @return
1838    */
1839   protected String printNonInterleaved(SequenceI[] s)
1840   {
1841     int maxSequenceLength = getMaxSequenceLength(s);
1842     // approx
1843     int numLines = maxSequenceLength / positionsPerLine + 2 + s.length;
1844
1845     /*
1846      * Roughly size a buffer to hold the whole output
1847      */
1848     StringBuilder sb = new StringBuilder(numLines * positionsPerLine);
1849
1850     for (SequenceI seq : s)
1851     {
1852       printSequence(sb, seq);
1853     }
1854
1855     return new String(sb);
1856   }
1857
1858   /**
1859    * Append a formatted complete sequence to the string buffer
1860    * 
1861    * @param sb
1862    * @param seq
1863    */
1864   protected void printSequence(StringBuilder sb, SequenceI seq)
1865   {
1866     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1867     // round down to output a whole number of spaced blocks
1868     int chunksPerLine = positionsPerLine / spaceEvery;
1869
1870     sb.append(newline);
1871     sb.append(HASHSIGN + seq.getName()).append(newline);
1872     int startPos = 0;
1873     while (startPos < seq.getLength())
1874     {
1875       /*
1876        * print next line for this sequence
1877        */
1878       boolean firstChunk = true;
1879       int lastPos = startPos + positionsPerLine; // exclusive
1880       for (int j = 0; j < chunksPerLine; j++)
1881       {
1882         char[] subSequence = seq.getSequence(startPos,
1883                 Math.min(lastPos, startPos + spaceEvery));
1884         if (subSequence.length > 0)
1885         {
1886           if (!firstChunk)
1887           {
1888             sb.append(SPACE);
1889           }
1890           sb.append(subSequence);
1891           firstChunk = false;
1892         }
1893         startPos += subSequence.length;
1894       }
1895       // line end position (base 1) as a comment
1896       sb.append(SPACE).append(COMMENT_START).append(startPos)
1897               .append(COMMENT_END);
1898       sb.append(newline);
1899     }
1900   }
1901
1902   /**
1903    * Returns a formatted string like <br>
1904    * !Label aa_b_ ab_b_ <br>
1905    * where underscore represents no annotation, any other character a MEGA label
1906    * character <br>
1907    * Returns an empty string if there is no non-null annotation in the given
1908    * alignment range
1909    * 
1910    * @param fromPos
1911    *          start column of the alignment (base 0)
1912    * @param positions
1913    *          number of positions to output
1914    * @param labelWidth
1915    *          padded width of !Label statement to output
1916    * @return
1917    */
1918   protected String printLabel(int fromPos, int positions, int labelWidth)
1919   {
1920     int spaceEvery = this.nucleotide != null && this.nucleotide ? 3 : 10;
1921     String none = "";
1922     AlignmentAnnotation ann = findAnnotation(MEGA_ANNOTATION_LABEL);
1923     if (ann == null)
1924     {
1925       return none;
1926     }
1927
1928     StringBuilder sb = new StringBuilder(positions + 20);
1929     sb.append(String.format("%-" + labelWidth + "s ", BANG + LABEL));
1930     Annotation[] anns = annotations.get(0).annotations;
1931     int blockCharCount = 0;
1932     boolean annotationFound = false;
1933
1934     for (int i = fromPos; i < fromPos + positions; i++)
1935     {
1936       String label = String.valueOf(UNDERSCORE);
1937       if (i < anns.length && anns[i] != null)
1938       {
1939         label = anns[i].displayCharacter;
1940       }
1941       sb.append(label);
1942       if (label.charAt(0) != UNDERSCORE)
1943       {
1944         annotationFound = true;
1945       }
1946       // add a space after each block except the last
1947       if (++blockCharCount % spaceEvery == 0
1948               && (i < fromPos + positions - 1))
1949       {
1950         sb.append(SPACE);
1951       }
1952     }
1953     sb.append(SEMICOLON).append(newline);
1954
1955     return annotationFound ? sb.toString() : none;
1956   }
1957
1958   /**
1959    * Returns the first stored annotation found with the given label, or null
1960    * 
1961    * @param annotationLabel
1962    * @return
1963    */
1964   protected AlignmentAnnotation findAnnotation(String annotationLabel)
1965   {
1966     if (annotations != null)
1967     {
1968       for (AlignmentAnnotation ann : annotations)
1969       {
1970         if (annotationLabel.equals(ann.label))
1971         {
1972           return ann;
1973         }
1974       }
1975     }
1976     return null;
1977   }
1978
1979   /**
1980    * Flag this file as interleaved or not, based on data format. Throws an
1981    * exception if has previously been determined to be otherwise.
1982    * 
1983    * @param isIt
1984    * @param dataLine
1985    * @throws FileFormatException
1986    */
1987   protected void assertInterleaved(boolean isIt, String dataLine)
1988           throws FileFormatException
1989   {
1990     if (this.interleaved != null && isIt != this.interleaved.booleanValue())
1991     {
1992       throw new FileFormatException("Parse error: interleaved was " + !isIt
1993               + " but now seems to be " + isIt + ", at line: " + dataLine);
1994     }
1995     this.interleaved = new Boolean(isIt);
1996     setAlignmentProperty(PROP_INTERLEAVED, interleaved.toString());
1997   }
1998
1999   public boolean isInterleaved()
2000   {
2001     return this.interleaved == null ? false : this.interleaved
2002             .booleanValue();
2003   }
2004
2005   /**
2006    * Adds saved parsed values either as alignment properties, or (in some cases)
2007    * as specific member fields of the alignment
2008    */
2009   @Override
2010   public void addProperties(AlignmentI al)
2011   {
2012     super.addProperties(al);
2013     
2014     /*
2015      * record gap character specified, but convert to '-' if not one we support
2016      */
2017     al.setGapCharacter(Comparison.isGap(gapCharacter) ? gapCharacter
2018             : DEFAULT_GAP);
2019
2020     /*
2021      * warn if e.g. DataType=DNA but data is protein (or vice versa)
2022      */
2023     if (this.nucleotide != null && this.nucleotide != al.isNucleotide()) {
2024       System.err.println("Warning: " + this.title + " declared "
2025               + (nucleotide ? "" : " not ") + "nucleotide but it is"
2026               + (nucleotide ? " not" : ""));
2027     }
2028   }
2029
2030   /**
2031    * Print the given alignment in MEGA format. If the alignment was created by
2032    * parsing a MEGA file, it should have properties set (e.g. Title) which can
2033    * surface in the output.
2034    */
2035   @Override
2036   public String print(AlignmentI al)
2037   {
2038     this.nucleotide = al.isNucleotide();
2039
2040     /*
2041      * "MEGA Gene", "MEGA Domain" or "MEGA Label" annotations can be output
2042      */
2043     AlignmentAnnotation[] anns = al.getAlignmentAnnotation();
2044     if (anns != null)
2045     {
2046       this.annotations = new Vector<AlignmentAnnotation>();
2047       for (AlignmentAnnotation ann : anns)
2048       {
2049         annotations.add(ann);
2050       }
2051     }
2052
2053     String lineLength = (String) al.getProperty(PROP_LINELENGTH);
2054     this.positionsPerLine = lineLength == null ? DEFAULT_LINE_LENGTH : Integer
2055             .parseInt(lineLength);
2056
2057     /*
2058      * round down to a multiple of 3 positions per line for nucleotide
2059      */
2060     if (nucleotide)
2061     {
2062       positionsPerLine = positionsPerLine - (positionsPerLine % 3);
2063     }
2064
2065     String interleave = (String) al.getProperty(PROP_INTERLEAVED);
2066     if (interleave != null)
2067     {
2068       this.interleaved = Boolean.valueOf(interleave);
2069     }
2070
2071     String headers = printHeaders(al);
2072     return headers + print(al.getSequencesArray());
2073   }
2074
2075   /**
2076    * Returns the number of sequence positions output per line
2077    * 
2078    * @return
2079    */
2080   public int getPositionsPerLine()
2081   {
2082     return positionsPerLine;
2083   }
2084
2085   /**
2086    * Sets the number of sequence positions output per line. Note these will be
2087    * formatted in blocks of 3 (nucleotide) or 10 (peptide).
2088    * 
2089    * @param p
2090    */
2091   public void setPositionsPerLine(int p)
2092   {
2093     this.positionsPerLine = p;
2094   }
2095
2096   /**
2097    * Extracts the Coding / Noncoding property of a domain from MEGA Domain
2098    * annotation at the given column position, if present
2099    * 
2100    * @param column
2101    * @param domainAnnotation
2102    * @return
2103    */
2104   protected static String getPropertyFromAnnotation(int column,
2105           AlignmentAnnotation domainAnnotation)
2106   {
2107     /*
2108      * currently depends on parsing "Exon1 (Aspx Coding)" or similar
2109      */
2110     String property = null;
2111     Annotation domainAnn = domainAnnotation == null ? null
2112             : (column >= domainAnnotation.annotations.length ? null
2113                     : domainAnnotation.annotations[column]);
2114     String domain = domainAnn == null ? null : domainAnn.description;
2115     if (domain != null && domain.indexOf("(") > 0)
2116     {
2117       domain = domain.substring(domain.indexOf("(") +1);
2118       if (domain.indexOf(SPACE) > -1 && domain.endsWith(")")) {
2119         property = domain.substring(domain.indexOf(SPACE) + 1,
2120                 domain.length() - 1);
2121       }
2122     }
2123     return property;
2124   }
2125 }