JAL-2490 tidy code and tests for feature export in Jalview format
[jalview.git] / src / jalview / io / FeaturesFile.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ 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
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.io;
22
23 import jalview.analysis.AlignmentUtils;
24 import jalview.analysis.SequenceIdMatcher;
25 import jalview.api.AlignViewportI;
26 import jalview.api.FeatureColourI;
27 import jalview.api.FeaturesSourceI;
28 import jalview.datamodel.AlignedCodonFrame;
29 import jalview.datamodel.Alignment;
30 import jalview.datamodel.AlignmentI;
31 import jalview.datamodel.SequenceDummy;
32 import jalview.datamodel.SequenceFeature;
33 import jalview.datamodel.SequenceI;
34 import jalview.io.gff.GffHelperBase;
35 import jalview.io.gff.GffHelperFactory;
36 import jalview.io.gff.GffHelperI;
37 import jalview.schemes.FeatureColour;
38 import jalview.util.ColorUtils;
39 import jalview.util.MapList;
40 import jalview.util.ParseHtmlBodyAndLinks;
41 import jalview.util.StringUtils;
42
43 import java.awt.Color;
44 import java.io.IOException;
45 import java.util.ArrayList;
46 import java.util.Arrays;
47 import java.util.Collections;
48 import java.util.Comparator;
49 import java.util.HashMap;
50 import java.util.HashSet;
51 import java.util.List;
52 import java.util.Map;
53 import java.util.Map.Entry;
54 import java.util.Set;
55
56 /**
57  * Parses and writes features files, which may be in Jalview, GFF2 or GFF3
58  * format. These are tab-delimited formats but with differences in the use of
59  * columns.
60  * 
61  * A Jalview feature file may define feature colours and then declare that the
62  * remainder of the file is in GFF format with the line 'GFF'.
63  * 
64  * GFF3 files may include alignment mappings for features, which Jalview will
65  * attempt to model, and may include sequence data following a ##FASTA line.
66  * 
67  * 
68  * @author AMW
69  * @author jbprocter
70  * @author gmcarstairs
71  */
72 public class FeaturesFile extends AlignFile implements FeaturesSourceI
73 {
74   private static final String ID_NOT_SPECIFIED = "ID_NOT_SPECIFIED";
75
76   private static final String NOTE = "Note";
77
78   protected static final String TAB = "\t";
79
80   protected static final String GFF_VERSION = "##gff-version";
81
82   private static final Comparator<String> SORT_NULL_LAST = new Comparator<String>()
83   {
84     @Override
85     public int compare(String o1, String o2)
86     {
87       if (o1 == null)
88       {
89         return o2 == null ? 0 : 1;
90       }
91       return (o2 == null ? -1 : o1.compareTo(o2));
92     }
93   };
94
95   private AlignmentI lastmatchedAl = null;
96
97   private SequenceIdMatcher matcher = null;
98
99   protected AlignmentI dataset;
100
101   protected int gffVersion;
102
103   /**
104    * Creates a new FeaturesFile object.
105    */
106   public FeaturesFile()
107   {
108   }
109
110   /**
111    * Constructor which does not parse the file immediately
112    * 
113    * @param inFile
114    * @param paste
115    * @throws IOException
116    */
117   public FeaturesFile(String inFile, DataSourceType paste)
118           throws IOException
119   {
120     super(false, inFile, paste);
121   }
122
123   /**
124    * @param source
125    * @throws IOException
126    */
127   public FeaturesFile(FileParse source) throws IOException
128   {
129     super(source);
130   }
131
132   /**
133    * Constructor that optionally parses the file immediately
134    * 
135    * @param parseImmediately
136    * @param inFile
137    * @param type
138    * @throws IOException
139    */
140   public FeaturesFile(boolean parseImmediately, String inFile,
141           DataSourceType type)
142           throws IOException
143   {
144     super(parseImmediately, inFile, type);
145   }
146
147   /**
148    * Parse GFF or sequence features file using case-independent matching,
149    * discarding URLs
150    * 
151    * @param align
152    *          - alignment/dataset containing sequences that are to be annotated
153    * @param colours
154    *          - hashtable to store feature colour definitions
155    * @param removeHTML
156    *          - process html strings into plain text
157    * @return true if features were added
158    */
159   public boolean parse(AlignmentI align,
160           Map<String, FeatureColourI> colours, boolean removeHTML)
161   {
162     return parse(align, colours, removeHTML, false);
163   }
164
165   /**
166    * Extends the default addProperties by also adding peptide-to-cDNA mappings
167    * (if any) derived while parsing a GFF file
168    */
169   @Override
170   public void addProperties(AlignmentI al)
171   {
172     super.addProperties(al);
173     if (dataset != null && dataset.getCodonFrames() != null)
174     {
175       AlignmentI ds = (al.getDataset() == null) ? al : al.getDataset();
176       for (AlignedCodonFrame codons : dataset.getCodonFrames())
177       {
178         ds.addCodonFrame(codons);
179       }
180     }
181   }
182
183   /**
184    * Parse GFF or Jalview format sequence features file
185    * 
186    * @param align
187    *          - alignment/dataset containing sequences that are to be annotated
188    * @param colours
189    *          - hashtable to store feature colour definitions
190    * @param removeHTML
191    *          - process html strings into plain text
192    * @param relaxedIdmatching
193    *          - when true, ID matches to compound sequence IDs are allowed
194    * @return true if features were added
195    */
196   public boolean parse(AlignmentI align,
197           Map<String, FeatureColourI> colours, boolean removeHTML,
198           boolean relaxedIdmatching)
199   {
200     Map<String, String> gffProps = new HashMap<String, String>();
201     /*
202      * keep track of any sequences we try to create from the data
203      */
204     List<SequenceI> newseqs = new ArrayList<SequenceI>();
205
206     String line = null;
207     try
208     {
209       String[] gffColumns;
210       String featureGroup = null;
211
212       while ((line = nextLine()) != null)
213       {
214         // skip comments/process pragmas
215         if (line.length() == 0 || line.startsWith("#"))
216         {
217           if (line.toLowerCase().startsWith("##"))
218           {
219             processGffPragma(line, gffProps, align, newseqs);
220           }
221           continue;
222         }
223
224         gffColumns = line.split("\\t"); // tab as regex
225         if (gffColumns.length == 1)
226         {
227           if (line.trim().equalsIgnoreCase("GFF"))
228           {
229             /*
230              * Jalview features file with appended GFF
231              * assume GFF2 (though it may declare ##gff-version 3)
232              */
233             gffVersion = 2;
234             continue;
235           }
236         }
237
238         if (gffColumns.length > 1 && gffColumns.length < 4)
239         {
240           /*
241            * if 2 or 3 tokens, we anticipate either 'startgroup', 'endgroup' or
242            * a feature type colour specification
243            */
244           String ft = gffColumns[0];
245           if (ft.equalsIgnoreCase("startgroup"))
246           {
247             featureGroup = gffColumns[1];
248           }
249           else if (ft.equalsIgnoreCase("endgroup"))
250           {
251             // We should check whether this is the current group,
252             // but at present there's no way of showing more than 1 group
253             featureGroup = null;
254           }
255           else
256           {
257             String colscheme = gffColumns[1];
258             FeatureColourI colour = FeatureColour
259                     .parseJalviewFeatureColour(colscheme);
260             if (colour != null)
261             {
262               colours.put(ft, colour);
263             }
264           }
265           continue;
266         }
267
268         /*
269          * if not a comment, GFF pragma, startgroup, endgroup or feature
270          * colour specification, that just leaves a feature details line
271          * in either Jalview or GFF format
272          */
273         if (gffVersion == 0)
274         {
275           parseJalviewFeature(line, gffColumns, align, colours, removeHTML,
276                   relaxedIdmatching, featureGroup);
277         }
278         else
279         {
280           parseGff(gffColumns, align, relaxedIdmatching, newseqs);
281         }
282       }
283       resetMatcher();
284     } catch (Exception ex)
285     {
286       // should report somewhere useful for UI if necessary
287       warningMessage = ((warningMessage == null) ? "" : warningMessage)
288               + "Parsing error at\n" + line;
289       System.out.println("Error parsing feature file: " + ex + "\n" + line);
290       ex.printStackTrace(System.err);
291       resetMatcher();
292       return false;
293     }
294
295     /*
296      * experimental - add any dummy sequences with features to the alignment
297      * - we need them for Ensembl feature extraction - though maybe not otherwise
298      */
299     for (SequenceI newseq : newseqs)
300     {
301       if (newseq.getFeatures().hasFeatures())
302       {
303         align.addSequence(newseq);
304       }
305     }
306     return true;
307   }
308
309   /**
310    * Try to parse a Jalview format feature specification and add it as a
311    * sequence feature to any matching sequences in the alignment. Returns true
312    * if successful (a feature was added), or false if not.
313    * 
314    * @param line
315    * @param gffColumns
316    * @param alignment
317    * @param featureColours
318    * @param removeHTML
319    * @param relaxedIdmatching
320    * @param featureGroup
321    */
322   protected boolean parseJalviewFeature(String line, String[] gffColumns,
323           AlignmentI alignment, Map<String, FeatureColourI> featureColours,
324           boolean removeHTML, boolean relaxedIdMatching, String featureGroup)
325   {
326     /*
327      * tokens: description seqid seqIndex start end type [score]
328      */
329     if (gffColumns.length < 6)
330     {
331       System.err.println("Ignoring feature line '" + line
332               + "' with too few columns (" + gffColumns.length + ")");
333       return false;
334     }
335     String desc = gffColumns[0];
336     String seqId = gffColumns[1];
337     SequenceI seq = findSequence(seqId, alignment, null, relaxedIdMatching);
338
339     if (!ID_NOT_SPECIFIED.equals(seqId))
340     {
341       seq = findSequence(seqId, alignment, null, relaxedIdMatching);
342     }
343     else
344     {
345       seqId = null;
346       seq = null;
347       String seqIndex = gffColumns[2];
348       try
349       {
350         int idx = Integer.parseInt(seqIndex);
351         seq = alignment.getSequenceAt(idx);
352       } catch (NumberFormatException ex)
353       {
354         System.err.println("Invalid sequence index: " + seqIndex);
355       }
356     }
357
358     if (seq == null)
359     {
360       System.out.println("Sequence not found: " + line);
361       return false;
362     }
363
364     int startPos = Integer.parseInt(gffColumns[3]);
365     int endPos = Integer.parseInt(gffColumns[4]);
366
367     String ft = gffColumns[5];
368
369     if (!featureColours.containsKey(ft))
370     {
371       /* 
372        * Perhaps an old style groups file with no colours -
373        * synthesize a colour from the feature type
374        */
375       Color colour = ColorUtils.createColourFromName(ft);
376       featureColours.put(ft, new FeatureColour(colour));
377     }
378     SequenceFeature sf = new SequenceFeature(ft, desc, "", startPos,
379             endPos, featureGroup);
380     if (gffColumns.length > 6)
381     {
382       float score = Float.NaN;
383       try
384       {
385         score = new Float(gffColumns[6]).floatValue();
386         // update colourgradient bounds if allowed to
387       } catch (NumberFormatException ex)
388       {
389         // leave as NaN
390       }
391       sf.setScore(score);
392     }
393
394     parseDescriptionHTML(sf, removeHTML);
395
396     seq.addSequenceFeature(sf);
397
398     while (seqId != null
399             && (seq = alignment.findName(seq, seqId, false)) != null)
400     {
401       seq.addSequenceFeature(new SequenceFeature(sf));
402     }
403     return true;
404   }
405
406   /**
407    * clear any temporary handles used to speed up ID matching
408    */
409   protected void resetMatcher()
410   {
411     lastmatchedAl = null;
412     matcher = null;
413   }
414
415   /**
416    * Returns a sequence matching the given id, as follows
417    * <ul>
418    * <li>strict matching is on exact sequence name</li>
419    * <li>relaxed matching allows matching on a token within the sequence name,
420    * or a dbxref</li>
421    * <li>first tries to find a match in the alignment sequences</li>
422    * <li>else tries to find a match in the new sequences already generated while
423    * parsing the features file</li>
424    * <li>else creates a new placeholder sequence, adds it to the new sequences
425    * list, and returns it</li>
426    * </ul>
427    * 
428    * @param seqId
429    * @param align
430    * @param newseqs
431    * @param relaxedIdMatching
432    * 
433    * @return
434    */
435   protected SequenceI findSequence(String seqId, AlignmentI align,
436           List<SequenceI> newseqs, boolean relaxedIdMatching)
437   {
438     // TODO encapsulate in SequenceIdMatcher, share the matcher
439     // with the GffHelper (removing code duplication)
440     SequenceI match = null;
441     if (relaxedIdMatching)
442     {
443       if (lastmatchedAl != align)
444       {
445         lastmatchedAl = align;
446         matcher = new SequenceIdMatcher(align.getSequencesArray());
447         if (newseqs != null)
448         {
449           matcher.addAll(newseqs);
450         }
451       }
452       match = matcher.findIdMatch(seqId);
453     }
454     else
455     {
456       match = align.findName(seqId, true);
457       if (match == null && newseqs != null)
458       {
459         for (SequenceI m : newseqs)
460         {
461           if (seqId.equals(m.getName()))
462           {
463             return m;
464           }
465         }
466       }
467
468     }
469     if (match == null && newseqs != null)
470     {
471       match = new SequenceDummy(seqId);
472       if (relaxedIdMatching)
473       {
474         matcher.addAll(Arrays.asList(new SequenceI[] { match }));
475       }
476       // add dummy sequence to the newseqs list
477       newseqs.add(match);
478     }
479     return match;
480   }
481
482   public void parseDescriptionHTML(SequenceFeature sf, boolean removeHTML)
483   {
484     if (sf.getDescription() == null)
485     {
486       return;
487     }
488     ParseHtmlBodyAndLinks parsed = new ParseHtmlBodyAndLinks(
489             sf.getDescription(), removeHTML, newline);
490
491     if (removeHTML)
492     {
493       sf.setDescription(parsed.getNonHtmlContent());
494     }
495
496     for (String link : parsed.getLinks())
497     {
498       sf.addLink(link);
499     }
500   }
501
502   /**
503    * Returns contents of a Jalview format features file
504    * 
505    * @param sequences
506    *          source of features
507    * @param visible
508    *          map of colour for each visible feature type
509    * @param includeNonPositional
510    *          if true, include non-positional features (regardless of group or
511    *          type)
512    * @return
513    */
514   public String printJalviewFormat(SequenceI[] sequences,
515           Map<String, FeatureColourI> visible, boolean includeNonPositional)
516   {
517     if (!includeNonPositional && (visible == null || visible.isEmpty()))
518     {
519       // no point continuing.
520       return "No Features Visible";
521     }
522
523     /*
524      * write out feature colours (if we know them)
525      */
526     // TODO: decide if feature links should also be written here ?
527     StringBuilder out = new StringBuilder(256);
528     if (visible != null)
529     {
530       for (Entry<String, FeatureColourI> featureColour : visible.entrySet())
531       {
532         FeatureColourI colour = featureColour.getValue();
533         out.append(colour.toJalviewFormat(featureColour.getKey())).append(
534                 newline);
535       }
536     }
537
538     // Work out which groups are both present and visible
539     Set<String> groups = new HashSet<String>();
540     String[] types = visible == null ? new String[0] : visible.keySet()
541             .toArray(new String[visible.keySet().size()]);
542
543     for (int i = 0; i < sequences.length; i++)
544     {
545       groups.addAll(sequences[i].getFeatures()
546               .getFeatureGroups(true, types));
547       if (includeNonPositional)
548       {
549         groups.addAll(sequences[i].getFeatures().getFeatureGroups(false,
550                 types));
551       }
552     }
553
554     /*
555      * sort distinct groups so null group is output last
556      */
557     List<String> sortedGroups = new ArrayList<String>(groups);
558     Collections.sort(sortedGroups, SORT_NULL_LAST);
559
560     // TODO check where null group should be output
561     boolean foundSome = false;
562     for (String group : sortedGroups)
563     {
564       if (group != null)
565       {
566         out.append(newline);
567         out.append("STARTGROUP").append(TAB);
568         out.append(group);
569         out.append(newline);
570       }
571
572       /*
573        * output features within groups (non-positional first if wanted)
574        */
575       for (int i = 0; i < sequences.length; i++)
576       {
577         List<SequenceFeature> features = new ArrayList<SequenceFeature>();
578         if (includeNonPositional)
579         {
580           features.addAll(sequences[i].getFeatures().getFeaturesForGroup(
581                   false, group, types));
582         }
583         if (types.length > 0)
584         {
585           features.addAll(sequences[i].getFeatures().getFeaturesForGroup(
586                   true, group, types));
587         }
588
589         for (SequenceFeature sequenceFeature : features)
590         {
591           // we have features to output
592           foundSome = true;
593           if (sequenceFeature.description == null
594                   || sequenceFeature.description.equals(""))
595           {
596             out.append(sequenceFeature.type).append(TAB);
597           }
598           else
599           {
600             if (sequenceFeature.links != null
601                     && sequenceFeature.getDescription().indexOf("<html>") == -1)
602             {
603               out.append("<html>");
604             }
605
606             out.append(sequenceFeature.description);
607             if (sequenceFeature.links != null)
608             {
609               for (int l = 0; l < sequenceFeature.links.size(); l++)
610               {
611                 String label = sequenceFeature.links.elementAt(l);
612                 String href = label.substring(label.indexOf("|") + 1);
613                 label = label.substring(0, label.indexOf("|"));
614
615                 if (sequenceFeature.description.indexOf(href) == -1)
616                 {
617                   out.append(" <a href=\"" + href + "\">" + label + "</a>");
618                 }
619               }
620
621               if (sequenceFeature.getDescription().indexOf("</html>") == -1)
622               {
623                 out.append("</html>");
624               }
625             }
626
627             out.append(TAB);
628           }
629           out.append(sequences[i].getName());
630           out.append("\t-1\t");
631           out.append(sequenceFeature.begin);
632           out.append(TAB);
633           out.append(sequenceFeature.end);
634           out.append(TAB);
635           out.append(sequenceFeature.type);
636           if (!Float.isNaN(sequenceFeature.score))
637           {
638             out.append(TAB);
639             out.append(sequenceFeature.score);
640           }
641           out.append(newline);
642         }
643       }
644
645       if (group != null)
646       {
647         out.append("ENDGROUP").append(TAB);
648         out.append(group);
649         out.append(newline);
650       }
651     }
652
653     return foundSome ? out.toString() : "No Features Visible";
654   }
655
656   /**
657    * Parse method that is called when a GFF file is dragged to the desktop
658    */
659   @Override
660   public void parse()
661   {
662     AlignViewportI av = getViewport();
663     if (av != null)
664     {
665       if (av.getAlignment() != null)
666       {
667         dataset = av.getAlignment().getDataset();
668       }
669       if (dataset == null)
670       {
671         // working in the applet context ?
672         dataset = av.getAlignment();
673       }
674     }
675     else
676     {
677       dataset = new Alignment(new SequenceI[] {});
678     }
679
680     Map<String, FeatureColourI> featureColours = new HashMap<String, FeatureColourI>();
681     boolean parseResult = parse(dataset, featureColours, false, true);
682     if (!parseResult)
683     {
684       // pass error up somehow
685     }
686     if (av != null)
687     {
688       // update viewport with the dataset data ?
689     }
690     else
691     {
692       setSeqs(dataset.getSequencesArray());
693     }
694   }
695
696   /**
697    * Implementation of unused abstract method
698    * 
699    * @return error message
700    */
701   @Override
702   public String print(SequenceI[] sqs, boolean jvsuffix)
703   {
704     System.out.println("Use printGffFormat() or printJalviewFormat()");
705     return null;
706   }
707
708   /**
709    * Returns features output in GFF2 format
710    * 
711    * @param sequences
712    *          the sequences whose features are to be output
713    * @param visible
714    *          a map whose keys are the type names of visible features
715    * @param includeNonPositionalFeatures
716    * @return
717    */
718   public String printGffFormat(SequenceI[] sequences,
719           Map<String, FeatureColourI> visible,
720           boolean includeNonPositionalFeatures)
721   {
722     StringBuilder out = new StringBuilder(256);
723     int version = gffVersion == 0 ? 2 : gffVersion;
724     out.append(String.format("%s %d\n", GFF_VERSION, version));
725     String source;
726     boolean isnonpos;
727     for (SequenceI seq : sequences)
728     {
729       SequenceFeature[] features = seq.getSequenceFeatures();
730       if (features != null)
731       {
732         for (SequenceFeature sf : features)
733         {
734           isnonpos = sf.begin == 0 && sf.end == 0;
735           if (!includeNonPositionalFeatures && isnonpos)
736           {
737             /*
738              * ignore non-positional features if not wanted
739              */
740             continue;
741           }
742           if (!isnonpos && !visible.containsKey(sf.type))
743           {
744             /*
745              * ignore not visible features if not wanted
746              */
747             continue;
748           }
749
750           source = sf.featureGroup;
751           if (source == null)
752           {
753             source = sf.getDescription();
754           }
755
756           out.append(seq.getName());
757           out.append(TAB);
758           out.append(source);
759           out.append(TAB);
760           out.append(sf.type);
761           out.append(TAB);
762           out.append(sf.begin);
763           out.append(TAB);
764           out.append(sf.end);
765           out.append(TAB);
766           out.append(sf.score);
767           out.append(TAB);
768
769           int strand = sf.getStrand();
770           out.append(strand == 1 ? "+" : (strand == -1 ? "-" : "."));
771           out.append(TAB);
772
773           String phase = sf.getPhase();
774           out.append(phase == null ? "." : phase);
775
776           // miscellaneous key-values (GFF column 9)
777           String attributes = sf.getAttributes();
778           if (attributes != null)
779           {
780             out.append(TAB).append(attributes);
781           }
782
783           out.append(newline);
784         }
785       }
786     }
787
788     return out.toString();
789   }
790
791   /**
792    * Returns a mapping given list of one or more Align descriptors (exonerate
793    * format)
794    * 
795    * @param alignedRegions
796    *          a list of "Align fromStart toStart fromCount"
797    * @param mapIsFromCdna
798    *          if true, 'from' is dna, else 'from' is protein
799    * @param strand
800    *          either 1 (forward) or -1 (reverse)
801    * @return
802    * @throws IOException
803    */
804   protected MapList constructCodonMappingFromAlign(
805           List<String> alignedRegions, boolean mapIsFromCdna, int strand)
806           throws IOException
807   {
808     if (strand == 0)
809     {
810       throw new IOException(
811               "Invalid strand for a codon mapping (cannot be 0)");
812     }
813     int regions = alignedRegions.size();
814     // arrays to hold [start, end] for each aligned region
815     int[] fromRanges = new int[regions * 2]; // from dna
816     int[] toRanges = new int[regions * 2]; // to protein
817     int fromRangesIndex = 0;
818     int toRangesIndex = 0;
819
820     for (String range : alignedRegions)
821     {
822       /* 
823        * Align mapFromStart mapToStart mapFromCount
824        * e.g. if mapIsFromCdna
825        *     Align 11270 143 120
826        * means:
827        *     120 bases from pos 11270 align to pos 143 in peptide
828        * if !mapIsFromCdna this would instead be
829        *     Align 143 11270 40 
830        */
831       String[] tokens = range.split(" ");
832       if (tokens.length != 3)
833       {
834         throw new IOException("Wrong number of fields for Align");
835       }
836       int fromStart = 0;
837       int toStart = 0;
838       int fromCount = 0;
839       try
840       {
841         fromStart = Integer.parseInt(tokens[0]);
842         toStart = Integer.parseInt(tokens[1]);
843         fromCount = Integer.parseInt(tokens[2]);
844       } catch (NumberFormatException nfe)
845       {
846         throw new IOException("Invalid number in Align field: "
847                 + nfe.getMessage());
848       }
849
850       /*
851        * Jalview always models from dna to protein, so adjust values if the
852        * GFF mapping is from protein to dna
853        */
854       if (!mapIsFromCdna)
855       {
856         fromCount *= 3;
857         int temp = fromStart;
858         fromStart = toStart;
859         toStart = temp;
860       }
861       fromRanges[fromRangesIndex++] = fromStart;
862       fromRanges[fromRangesIndex++] = fromStart + strand * (fromCount - 1);
863
864       /*
865        * If a codon has an intron gap, there will be contiguous 'toRanges';
866        * this is handled for us by the MapList constructor. 
867        * (It is not clear that exonerate ever generates this case)  
868        */
869       toRanges[toRangesIndex++] = toStart;
870       toRanges[toRangesIndex++] = toStart + (fromCount - 1) / 3;
871     }
872
873     return new MapList(fromRanges, toRanges, 3, 1);
874   }
875
876   /**
877    * Parse a GFF format feature. This may include creating a 'dummy' sequence to
878    * hold the feature, or for its mapped sequence, or both, to be resolved
879    * either later in the GFF file (##FASTA section), or when the user loads
880    * additional sequences.
881    * 
882    * @param gffColumns
883    * @param alignment
884    * @param relaxedIdMatching
885    * @param newseqs
886    * @return
887    */
888   protected SequenceI parseGff(String[] gffColumns, AlignmentI alignment,
889           boolean relaxedIdMatching, List<SequenceI> newseqs)
890   {
891     /*
892      * GFF: seqid source type start end score strand phase [attributes]
893      */
894     if (gffColumns.length < 5)
895     {
896       System.err.println("Ignoring GFF feature line with too few columns ("
897               + gffColumns.length + ")");
898       return null;
899     }
900
901     /*
902      * locate referenced sequence in alignment _or_ 
903      * as a forward or external reference (SequenceDummy)
904      */
905     String seqId = gffColumns[0];
906     SequenceI seq = findSequence(seqId, alignment, newseqs,
907             relaxedIdMatching);
908
909     SequenceFeature sf = null;
910     GffHelperI helper = GffHelperFactory.getHelper(gffColumns);
911     if (helper != null)
912     {
913       try
914       {
915         sf = helper.processGff(seq, gffColumns, alignment, newseqs,
916                 relaxedIdMatching);
917         if (sf != null)
918         {
919           seq.addSequenceFeature(sf);
920           while ((seq = alignment.findName(seq, seqId, true)) != null)
921           {
922             seq.addSequenceFeature(new SequenceFeature(sf));
923           }
924         }
925       } catch (IOException e)
926       {
927         System.err.println("GFF parsing failed with: " + e.getMessage());
928         return null;
929       }
930     }
931
932     return seq;
933   }
934
935   /**
936    * Process the 'column 9' data of the GFF file. This is less formally defined,
937    * and its interpretation will vary depending on the tool that has generated
938    * it.
939    * 
940    * @param attributes
941    * @param sf
942    */
943   protected void processGffColumnNine(String attributes, SequenceFeature sf)
944   {
945     sf.setAttributes(attributes);
946
947     /*
948      * Parse attributes in column 9 and add them to the sequence feature's 
949      * 'otherData' table; use Note as a best proxy for description
950      */
951     char nameValueSeparator = gffVersion == 3 ? '=' : ' ';
952     // TODO check we don't break GFF2 values which include commas here
953     Map<String, List<String>> nameValues = GffHelperBase
954             .parseNameValuePairs(attributes, ";", nameValueSeparator, ",");
955     for (Entry<String, List<String>> attr : nameValues.entrySet())
956     {
957       String values = StringUtils.listToDelimitedString(attr.getValue(),
958               "; ");
959       sf.setValue(attr.getKey(), values);
960       if (NOTE.equals(attr.getKey()))
961       {
962         sf.setDescription(values);
963       }
964     }
965   }
966
967   /**
968    * After encountering ##fasta in a GFF3 file, process the remainder of the
969    * file as FAST sequence data. Any placeholder sequences created during
970    * feature parsing are updated with the actual sequences.
971    * 
972    * @param align
973    * @param newseqs
974    * @throws IOException
975    */
976   protected void processAsFasta(AlignmentI align, List<SequenceI> newseqs)
977           throws IOException
978   {
979     try
980     {
981       mark();
982     } catch (IOException q)
983     {
984     }
985     FastaFile parser = new FastaFile(this);
986     List<SequenceI> includedseqs = parser.getSeqs();
987
988     SequenceIdMatcher smatcher = new SequenceIdMatcher(newseqs);
989
990     /*
991      * iterate over includedseqs, and replacing matching ones with newseqs
992      * sequences. Generic iterator not used here because we modify
993      * includedseqs as we go
994      */
995     for (int p = 0, pSize = includedseqs.size(); p < pSize; p++)
996     {
997       // search for any dummy seqs that this sequence can be used to update
998       SequenceI includedSeq = includedseqs.get(p);
999       SequenceI dummyseq = smatcher.findIdMatch(includedSeq);
1000       if (dummyseq != null && dummyseq instanceof SequenceDummy)
1001       {
1002         // probably have the pattern wrong
1003         // idea is that a flyweight proxy for a sequence ID can be created for
1004         // 1. stable reference creation
1005         // 2. addition of annotation
1006         // 3. future replacement by a real sequence
1007         // current pattern is to create SequenceDummy objects - a convenience
1008         // constructor for a Sequence.
1009         // problem is that when promoted to a real sequence, all references
1010         // need to be updated somehow. We avoid that by keeping the same object.
1011         ((SequenceDummy) dummyseq).become(includedSeq);
1012         dummyseq.createDatasetSequence();
1013
1014         /*
1015          * Update mappings so they are now to the dataset sequence
1016          */
1017         for (AlignedCodonFrame mapping : align.getCodonFrames())
1018         {
1019           mapping.updateToDataset(dummyseq);
1020         }
1021
1022         /*
1023          * replace parsed sequence with the realised forward reference
1024          */
1025         includedseqs.set(p, dummyseq);
1026
1027         /*
1028          * and remove from the newseqs list
1029          */
1030         newseqs.remove(dummyseq);
1031       }
1032     }
1033
1034     /*
1035      * finally add sequences to the dataset
1036      */
1037     for (SequenceI seq : includedseqs)
1038     {
1039       // experimental: mapping-based 'alignment' to query sequence
1040       AlignmentUtils.alignSequenceAs(seq, align,
1041               String.valueOf(align.getGapCharacter()), false, true);
1042
1043       // rename sequences if GFF handler requested this
1044       // TODO a more elegant way e.g. gffHelper.postProcess(newseqs) ?
1045       SequenceFeature[] sfs = seq.getSequenceFeatures();
1046       if (sfs != null)
1047       {
1048         String newName = (String) sfs[0].getValue(GffHelperI.RENAME_TOKEN);
1049         if (newName != null)
1050         {
1051           seq.setName(newName);
1052         }
1053       }
1054       align.addSequence(seq);
1055     }
1056   }
1057
1058   /**
1059    * Process a ## directive
1060    * 
1061    * @param line
1062    * @param gffProps
1063    * @param align
1064    * @param newseqs
1065    * @throws IOException
1066    */
1067   protected void processGffPragma(String line,
1068           Map<String, String> gffProps, AlignmentI align,
1069           List<SequenceI> newseqs) throws IOException
1070   {
1071     line = line.trim();
1072     if ("###".equals(line))
1073     {
1074       // close off any open 'forward references'
1075       return;
1076     }
1077
1078     String[] tokens = line.substring(2).split(" ");
1079     String pragma = tokens[0];
1080     String value = tokens.length == 1 ? null : tokens[1];
1081
1082     if ("gff-version".equalsIgnoreCase(pragma))
1083     {
1084       if (value != null)
1085       {
1086         try
1087         {
1088           // value may be e.g. "3.1.2"
1089           gffVersion = Integer.parseInt(value.split("\\.")[0]);
1090         } catch (NumberFormatException e)
1091         {
1092           // ignore
1093         }
1094       }
1095     }
1096     else if ("sequence-region".equalsIgnoreCase(pragma))
1097     {
1098       // could capture <seqid start end> if wanted here
1099     }
1100     else if ("feature-ontology".equalsIgnoreCase(pragma))
1101     {
1102       // should resolve against the specified feature ontology URI
1103     }
1104     else if ("attribute-ontology".equalsIgnoreCase(pragma))
1105     {
1106       // URI of attribute ontology - not currently used in GFF3
1107     }
1108     else if ("source-ontology".equalsIgnoreCase(pragma))
1109     {
1110       // URI of source ontology - not currently used in GFF3
1111     }
1112     else if ("species-build".equalsIgnoreCase(pragma))
1113     {
1114       // save URI of specific NCBI taxon version of annotations
1115       gffProps.put("species-build", value);
1116     }
1117     else if ("fasta".equalsIgnoreCase(pragma))
1118     {
1119       // process the rest of the file as a fasta file and replace any dummy
1120       // sequence IDs
1121       processAsFasta(align, newseqs);
1122     }
1123     else
1124     {
1125       System.err.println("Ignoring unknown pragma: " + line);
1126     }
1127   }
1128 }