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