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