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.FeatureRenderer;
  28 import jalview.api.FeaturesSourceI;
  29 import jalview.datamodel.AlignedCodonFrame;
  30 import jalview.datamodel.Alignment;
  31 import jalview.datamodel.AlignmentI;
  32 import jalview.datamodel.MappedFeatures;
  33 import jalview.datamodel.SequenceDummy;
  34 import jalview.datamodel.SequenceFeature;
  35 import jalview.datamodel.SequenceI;
  36 import jalview.datamodel.features.FeatureMatcherSet;
  37 import jalview.datamodel.features.FeatureMatcherSetI;
  38 import jalview.gui.Desktop;
  39 import jalview.io.gff.GffHelperFactory;
  40 import jalview.io.gff.GffHelperI;
  41 import jalview.schemes.FeatureColour;
  42 import jalview.util.ColorUtils;
  43 import jalview.util.MapList;
  44 import jalview.util.ParseHtmlBodyAndLinks;
  45 import jalview.util.StringUtils;
  46
  47 import java.awt.Color;
  48 import java.io.IOException;
  49 import java.util.ArrayList;
  50 import java.util.Arrays;
  51 import java.util.Collections;
  52 import java.util.HashMap;
  53 import java.util.LinkedHashMap;
  54 import java.util.List;
  55 import java.util.Map;
  56 import java.util.Map.Entry;
  57 import java.util.TreeMap;
  58
  59 /**
  60  * Parses and writes features files, which may be in Jalview, GFF2 or GFF3
  61  * format. These are tab-delimited formats but with differences in the use of
  62  * columns.
  63  *
  64  * A Jalview feature file may define feature colours and then declare that the
  65  * remainder of the file is in GFF format with the line 'GFF'.
  66  *
  67  * GFF3 files may include alignment mappings for features, which Jalview will
  68  * attempt to model, and may include sequence data following a ##FASTA line.
  69  *
  70  *
  71  * @author AMW
  72  * @author jbprocter
  73  * @author gmcarstairs
  74  */
  75 public class FeaturesFile extends AlignFile implements FeaturesSourceI
  76 {
  77   private static final String EQUALS = "=";
  78
  79   private static final String TAB_REGEX = "\\t";
  80
  81   private static final String STARTGROUP = "STARTGROUP";
  82
  83   private static final String ENDGROUP = "ENDGROUP";
  84
  85   private static final String STARTFILTERS = "STARTFILTERS";
  86
  87   private static final String ENDFILTERS = "ENDFILTERS";
  88
  89   private static final String ID_NOT_SPECIFIED = "ID_NOT_SPECIFIED";
  90
  91   protected static final String GFF_VERSION = "##gff-version";
  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 file File or String filename
 112    * @param paste
 113    * @throws IOException
 114    */
 115   public FeaturesFile(Object file, DataSourceType paste)
 116           throws IOException
 117   {
 118     super(false, file, 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 file
 135    * @param type
 136    * @throws IOException
 137    */
 138   public FeaturesFile(boolean parseImmediately, Object file,
 139           DataSourceType type) throws IOException
 140   {
 141     super(parseImmediately, file, type);
 142   }
 143
 144   /**
 145    * Parse GFF or sequence features file using case-independent matching,
 146    * discarding URLs
 147    *
 148    * @param align
 149    *          - alignment/dataset containing sequences that are to be annotated
 150    * @param colours
 151    *          - hashtable to store feature colour definitions
 152    * @param removeHTML
 153    *          - process html strings into plain text
 154    * @return true if features were added
 155    */
 156   public boolean parse(AlignmentI align,
 157           Map<String, FeatureColourI> colours, boolean removeHTML)
 158   {
 159     return parse(align, colours, removeHTML, false);
 160   }
 161
 162   /**
 163    * Extends the default addProperties by also adding peptide-to-cDNA mappings
 164    * (if any) derived while parsing a GFF file
 165    */
 166   @Override
 167   public void addProperties(AlignmentI al)
 168   {
 169     super.addProperties(al);
 170     if (dataset != null && dataset.getCodonFrames() != null)
 171     {
 172       AlignmentI ds = (al.getDataset() == null) ? al : al.getDataset();
 173       for (AlignedCodonFrame codons : dataset.getCodonFrames())
 174       {
 175         ds.addCodonFrame(codons);
 176       }
 177     }
 178   }
 179
 180   /**
 181    * Parse GFF or Jalview format sequence features file
 182    *
 183    * @param align
 184    *          - alignment/dataset containing sequences that are to be annotated
 185    * @param colours
 186    *          - map to store feature colour definitions
 187    * @param removeHTML
 188    *          - process html strings into plain text
 189    * @param relaxedIdmatching
 190    *          - when true, ID matches to compound sequence IDs are allowed
 191    * @return true if features were added
 192    */
 193   public boolean parse(AlignmentI align,
 194           Map<String, FeatureColourI> colours, boolean removeHTML,
 195           boolean relaxedIdmatching)
 196   {
 197     return parse(align, colours, null, removeHTML, relaxedIdmatching);
 198   }
 199
 200   /**
 201    * Parse GFF or Jalview format sequence features file
 202    *
 203    * @param align
 204    *          - alignment/dataset containing sequences that are to be annotated
 205    * @param colours
 206    *          - map to store feature colour definitions
 207    * @param filters
 208    *          - map to store feature filter definitions
 209    * @param removeHTML
 210    *          - process html strings into plain text
 211    * @param relaxedIdmatching
 212    *          - when true, ID matches to compound sequence IDs are allowed
 213    * @return true if features were added
 214    */
 215   public boolean parse(AlignmentI align,
 216           Map<String, FeatureColourI> colours,
 217           Map<String, FeatureMatcherSetI> filters, boolean removeHTML,
 218           boolean relaxedIdmatching)
 219   {
 220     Map<String, String> gffProps = new HashMap<>();
 221     /*
 222      * keep track of any sequences we try to create from the data
 223      */
 224     List<SequenceI> newseqs = new ArrayList<>();
 225
 226     String line = null;
 227     try
 228     {
 229       String[] gffColumns;
 230       String featureGroup = null;
 231
 232       while ((line = nextLine()) != null)
 233       {
 234         // skip comments/process pragmas
 235         if (line.length() == 0 || line.startsWith("#"))
 236         {
 237           if (line.toLowerCase().startsWith("##"))
 238           {
 239             processGffPragma(line, gffProps, align, newseqs);
 240           }
 241           continue;
 242         }
 243
 244         gffColumns = line.split(TAB_REGEX);
 245         if (gffColumns.length == 1)
 246         {
 247           if (line.trim().equalsIgnoreCase("GFF"))
 248           {
 249             /*
 250              * Jalview features file with appended GFF
 251              * assume GFF2 (though it may declare ##gff-version 3)
 252              */
 253             gffVersion = 2;
 254             continue;
 255           }
 256         }
 257
 258         if (gffColumns.length > 0 && gffColumns.length < 4)
 259         {
 260           /*
 261            * if 2 or 3 tokens, we anticipate either 'startgroup', 'endgroup' or
 262            * a feature type colour specification
 263            */
 264           String ft = gffColumns[0];
 265           if (ft.equalsIgnoreCase(STARTFILTERS))
 266           {
 267             parseFilters(filters);
 268             continue;
 269           }
 270           if (ft.equalsIgnoreCase(STARTGROUP))
 271           {
 272             featureGroup = gffColumns[1];
 273           }
 274           else if (ft.equalsIgnoreCase(ENDGROUP))
 275           {
 276             // We should check whether this is the current group,
 277             // but at present there's no way of showing more than 1 group
 278             featureGroup = null;
 279           }
 280           else
 281           {
 282             String colscheme = gffColumns[1];
 283             FeatureColourI colour = FeatureColour
 284                     .parseJalviewFeatureColour(colscheme);
 285             if (colour != null)
 286             {
 287               colours.put(ft, colour);
 288             }
 289           }
 290           continue;
 291         }
 292
 293         /*
 294          * if not a comment, GFF pragma, startgroup, endgroup or feature
 295          * colour specification, that just leaves a feature details line
 296          * in either Jalview or GFF format
 297          */
 298         if (gffVersion == 0)
 299         {
 300           parseJalviewFeature(line, gffColumns, align, colours, removeHTML,
 301                   relaxedIdmatching, featureGroup);
 302         }
 303         else
 304         {
 305           parseGff(gffColumns, align, relaxedIdmatching, newseqs);
 306         }
 307       }
 308       resetMatcher();
 309     } catch (Exception ex)
 310     {
 311       // should report somewhere useful for UI if necessary
 312       warningMessage = ((warningMessage == null) ? "" : warningMessage)
 313               + "Parsing error at\n" + line;
 314       System.out.println("Error parsing feature file: " + ex + "\n" + line);
 315       ex.printStackTrace(System.err);
 316       resetMatcher();
 317       return false;
 318     }
 319
 320     /*
 321      * experimental - add any dummy sequences with features to the alignment
 322      * - we need them for Ensembl feature extraction - though maybe not otherwise
 323      */
 324     for (SequenceI newseq : newseqs)
 325     {
 326       if (newseq.getFeatures().hasFeatures())
 327       {
 328         align.addSequence(newseq);
 329       }
 330     }
 331     return true;
 332   }
 333
 334   /**
 335    * Reads input lines from STARTFILTERS to ENDFILTERS and adds a feature type
 336    * filter to the map for each line parsed. After exit from this method,
 337    * nextLine() should return the line after ENDFILTERS (or we are already at
 338    * end of file if ENDFILTERS was missing).
 339    *
 340    * @param filters
 341    * @throws IOException
 342    */
 343   protected void parseFilters(Map<String, FeatureMatcherSetI> filters)
 344           throws IOException
 345   {
 346     String line;
 347     while ((line = nextLine()) != null)
 348     {
 349       if (line.toUpperCase().startsWith(ENDFILTERS))
 350       {
 351         return;
 352       }
 353       String[] tokens = line.split(TAB_REGEX);
 354       if (tokens.length != 2)
 355       {
 356         System.err.println(String.format("Invalid token count %d for %d",
 357                 tokens.length, line));
 358       }
 359       else
 360       {
 361         String featureType = tokens[0];
 362         FeatureMatcherSetI fm = FeatureMatcherSet.fromString(tokens[1]);
 363         if (fm != null && filters != null)
 364         {
 365           filters.put(featureType, fm);
 366         }
 367       }
 368     }
 369   }
 370
 371   /**
 372    * Try to parse a Jalview format feature specification and add it as a
 373    * sequence feature to any matching sequences in the alignment. Returns true
 374    * if successful (a feature was added), or false if not.
 375    *
 376    * @param line
 377    * @param gffColumns
 378    * @param alignment
 379    * @param featureColours
 380    * @param removeHTML
 381    * @param relaxedIdmatching
 382    * @param featureGroup
 383    */
 384   protected boolean parseJalviewFeature(String line, String[] gffColumns,
 385           AlignmentI alignment, Map<String, FeatureColourI> featureColours,
 386           boolean removeHTML, boolean relaxedIdMatching,
 387           String featureGroup)
 388   {
 389     /*
 390      * tokens: description seqid seqIndex start end type [score]
 391      */
 392     if (gffColumns.length < 6)
 393     {
 394       System.err.println("Ignoring feature line '" + line
 395               + "' with too few columns (" + gffColumns.length + ")");
 396       return false;
 397     }
 398     String desc = gffColumns[0];
 399     String seqId = gffColumns[1];
 400     SequenceI seq = findSequence(seqId, alignment, null, relaxedIdMatching);
 401
 402     if (!ID_NOT_SPECIFIED.equals(seqId))
 403     {
 404       seq = findSequence(seqId, alignment, null, relaxedIdMatching);
 405     }
 406     else
 407     {
 408       seqId = null;
 409       seq = null;
 410       String seqIndex = gffColumns[2];
 411       try
 412       {
 413         int idx = Integer.parseInt(seqIndex);
 414         seq = alignment.getSequenceAt(idx);
 415       } catch (NumberFormatException ex)
 416       {
 417         System.err.println("Invalid sequence index: " + seqIndex);
 418       }
 419     }
 420
 421     if (seq == null)
 422     {
 423       System.out.println("Sequence not found: " + line);
 424       return false;
 425     }
 426
 427     int startPos = Integer.parseInt(gffColumns[3]);
 428     int endPos = Integer.parseInt(gffColumns[4]);
 429
 430     String ft = gffColumns[5];
 431
 432     if (!featureColours.containsKey(ft))
 433     {
 434       /*
 435        * Perhaps an old style groups file with no colours -
 436        * synthesize a colour from the feature type
 437        */
 438       Color colour = ColorUtils.createColourFromName(ft);
 439       featureColours.put(ft, new FeatureColour(colour));
 440     }
 441     SequenceFeature sf = null;
 442     if (gffColumns.length > 6)
 443     {
 444       float score = Float.NaN;
 445       try
 446       {
 447         score = Float.valueOf(gffColumns[6]).floatValue();
 448       } catch (NumberFormatException ex)
 449       {
 450         sf = new SequenceFeature(ft, desc, startPos, endPos, featureGroup);
 451       }
 452       sf = new SequenceFeature(ft, desc, startPos, endPos, score,
 453               featureGroup);
 454     }
 455     else
 456     {
 457       sf = new SequenceFeature(ft, desc, startPos, endPos, featureGroup);
 458     }
 459
 460     parseDescriptionHTML(sf, removeHTML);
 461
 462     seq.addSequenceFeature(sf);
 463
 464     while (seqId != null
 465             && (seq = alignment.findName(seq, seqId, false)) != null)
 466     {
 467       seq.addSequenceFeature(new SequenceFeature(sf));
 468     }
 469     return true;
 470   }
 471
 472   /**
 473    * clear any temporary handles used to speed up ID matching
 474    */
 475   protected void resetMatcher()
 476   {
 477     lastmatchedAl = null;
 478     matcher = null;
 479   }
 480
 481   /**
 482    * Returns a sequence matching the given id, as follows
 483    * <ul>
 484    * <li>strict matching is on exact sequence name</li>
 485    * <li>relaxed matching allows matching on a token within the sequence name,
 486    * or a dbxref</li>
 487    * <li>first tries to find a match in the alignment sequences</li>
 488    * <li>else tries to find a match in the new sequences already generated while
 489    * parsing the features file</li>
 490    * <li>else creates a new placeholder sequence, adds it to the new sequences
 491    * list, and returns it</li>
 492    * </ul>
 493    *
 494    * @param seqId
 495    * @param align
 496    * @param newseqs
 497    * @param relaxedIdMatching
 498    *
 499    * @return
 500    */
 501   protected SequenceI findSequence(String seqId, AlignmentI align,
 502           List<SequenceI> newseqs, boolean relaxedIdMatching)
 503   {
 504     // TODO encapsulate in SequenceIdMatcher, share the matcher
 505     // with the GffHelper (removing code duplication)
 506     SequenceI match = null;
 507     if (relaxedIdMatching)
 508     {
 509       if (lastmatchedAl != align)
 510       {
 511         lastmatchedAl = align;
 512         matcher = new SequenceIdMatcher(align.getSequencesArray());
 513         if (newseqs != null)
 514         {
 515           matcher.addAll(newseqs);
 516         }
 517       }
 518       match = matcher.findIdMatch(seqId);
 519     }
 520     else
 521     {
 522       match = align.findName(seqId, true);
 523       if (match == null && newseqs != null)
 524       {
 525         for (SequenceI m : newseqs)
 526         {
 527           if (seqId.equals(m.getName()))
 528           {
 529             return m;
 530           }
 531         }
 532       }
 533
 534     }
 535     if (match == null && newseqs != null)
 536     {
 537       match = new SequenceDummy(seqId);
 538       if (relaxedIdMatching)
 539       {
 540         matcher.addAll(Arrays.asList(new SequenceI[] { match }));
 541       }
 542       // add dummy sequence to the newseqs list
 543       newseqs.add(match);
 544     }
 545     return match;
 546   }
 547
 548   public void parseDescriptionHTML(SequenceFeature sf, boolean removeHTML)
 549   {
 550     if (sf.getDescription() == null)
 551     {
 552       return;
 553     }
 554     ParseHtmlBodyAndLinks parsed = new ParseHtmlBodyAndLinks(
 555             sf.getDescription(), removeHTML, newline);
 556
 557     if (removeHTML)
 558     {
 559       sf.setDescription(parsed.getNonHtmlContent());
 560     }
 561
 562     for (String link : parsed.getLinks())
 563     {
 564       sf.addLink(link);
 565     }
 566   }
 567
 568   /**
 569    * Returns contents of a Jalview format features file, for visible features, as
 570    * filtered by type and group. Features with a null group are displayed if their
 571    * feature type is visible. Non-positional features may optionally be included
 572    * (with no check on type or group).
 573    *
 574    * @param sequences
 575    * @param fr
 576    * @param includeNonPositional
 577    *                               if true, include non-positional features
 578    *                               (regardless of group or type)
 579    * @param includeComplement
 580    *                               if true, include visible complementary
 581    *                               (CDS/protein) positional features, with
 582    *                               locations converted to local sequence
 583    *                               coordinates
 584    * @return
 585    */
 586   public String printJalviewFormat(SequenceI[] sequences,
 587           FeatureRenderer fr, boolean includeNonPositional,
 588           boolean includeComplement)
 589   {
 590     Map<String, FeatureColourI> visibleColours = fr
 591             .getDisplayedFeatureCols();
 592     Map<String, FeatureMatcherSetI> featureFilters = fr.getFeatureFilters();
 593
 594     /*
 595      * write out feature colours (if we know them)
 596      */
 597     // TODO: decide if feature links should also be written here ?
 598     StringBuilder out = new StringBuilder(256);
 599     if (visibleColours != null)
 600     {
 601       for (Entry<String, FeatureColourI> featureColour : visibleColours
 602               .entrySet())
 603       {
 604         FeatureColourI colour = featureColour.getValue();
 605         out.append(colour.toJalviewFormat(featureColour.getKey())).append(
 606                 newline);
 607       }
 608     }
 609
 610     String[] types = visibleColours == null ? new String[0]
 611             : visibleColours.keySet()
 612                     .toArray(new String[visibleColours.keySet().size()]);
 613
 614     /*
 615      * feature filters if any
 616      */
 617     outputFeatureFilters(out, visibleColours, featureFilters);
 618
 619     /*
 620      * output features within groups
 621      */
 622     int count = outputFeaturesByGroup(out, fr, types, sequences,
 623             includeNonPositional);
 624
 625     if (includeComplement)
 626     {
 627       count += outputComplementFeatures(out, fr, sequences);
 628     }
 629
 630     return count > 0 ? out.toString() : "No Features Visible";
 631   }
 632
 633   /**
 634    * Outputs any visible complementary (CDS/peptide) positional features as
 635    * Jalview format, within feature group. The coordinates of the linked features
 636    * are converted to the corresponding positions of the local sequences.
 637    *
 638    * @param out
 639    * @param fr
 640    * @param sequences
 641    * @return
 642    */
 643   private int outputComplementFeatures(StringBuilder out,
 644           FeatureRenderer fr, SequenceI[] sequences)
 645   {
 646     AlignViewportI comp = fr.getViewport().getCodingComplement();
 647     FeatureRenderer fr2 = Desktop.getAlignFrameFor(comp)
 648             .getFeatureRenderer();
 649
 650     /*
 651      * bin features by feature group and sequence
 652      */
 653     Map<String, Map<String, List<SequenceFeature>>> map = new TreeMap<>(
 654             String.CASE_INSENSITIVE_ORDER);
 655     int count = 0;
 656
 657     for (SequenceI seq : sequences)
 658     {
 659       /*
 660        * find complementary features
 661        */
 662       List<SequenceFeature> complementary = findComplementaryFeatures(seq,
 663               fr2);
 664       String seqName = seq.getName();
 665
 666       for (SequenceFeature sf : complementary)
 667       {
 668         String group = sf.getFeatureGroup();
 669         if (!map.containsKey(group))
 670         {
 671           map.put(group, new LinkedHashMap<>()); // preserves sequence order
 672         }
 673         Map<String, List<SequenceFeature>> groupFeatures = map.get(group);
 674         if (!groupFeatures.containsKey(seqName))
 675         {
 676           groupFeatures.put(seqName, new ArrayList<>());
 677         }
 678         List<SequenceFeature> foundFeatures = groupFeatures.get(seqName);
 679         foundFeatures.add(sf);
 680         count++;
 681       }
 682     }
 683
 684     /*
 685      * output features by group
 686      */
 687     for (Entry<String, Map<String, List<SequenceFeature>>> groupFeatures : map.entrySet())
 688     {
 689       out.append(newline);
 690       String group = groupFeatures.getKey();
 691       if (!"".equals(group))
 692       {
 693         out.append(STARTGROUP).append(TAB).append(group).append(newline);
 694       }
 695       Map<String, List<SequenceFeature>> seqFeaturesMap = groupFeatures
 696               .getValue();
 697       for (Entry<String, List<SequenceFeature>> seqFeatures : seqFeaturesMap
 698               .entrySet())
 699       {
 700         String sequenceName = seqFeatures.getKey();
 701         for (SequenceFeature sf : seqFeatures.getValue())
 702         {
 703           formatJalviewFeature(out, sequenceName, sf);
 704         }
 705       }
 706       if (!"".equals(group))
 707       {
 708         out.append(ENDGROUP).append(TAB).append(group).append(newline);
 709       }
 710     }
 711
 712     return count;
 713   }
 714
 715   /**
 716    * Answers a list of mapped features visible in the (CDS/protein) complement,
 717    * with feature positions translated to local sequence coordinates
 718    *
 719    * @param seq
 720    * @param fr2
 721    * @return
 722    */
 723   protected List<SequenceFeature> findComplementaryFeatures(SequenceI seq,
 724           FeatureRenderer fr2)
 725   {
 726     /*
 727      * avoid duplication of features (e.g. peptide feature
 728      * at all 3 mapped codon positions)
 729      */
 730     List<SequenceFeature> found = new ArrayList<>();
 731     List<SequenceFeature> complementary = new ArrayList<>();
 732
 733     for (int pos = seq.getStart(); pos <= seq.getEnd(); pos++)
 734     {
 735       MappedFeatures mf = fr2.findComplementFeaturesAtResidue(seq, pos);
 736
 737       if (mf != null)
 738       {
 739         MapList mapping = mf.mapping.getMap();
 740         for (SequenceFeature sf : mf.features)
 741         {
 742           /*
 743            * make a virtual feature with local coordinates
 744            */
 745           if (!found.contains(sf))
 746           {
 747             String group = sf.getFeatureGroup();
 748             if (group == null)
 749             {
 750               group = "";
 751             }
 752             found.add(sf);
 753             int begin = sf.getBegin();
 754             int end = sf.getEnd();
 755             int[] range = mf.mapping.getTo() == seq.getDatasetSequence()
 756                     ? mapping.locateInTo(begin, end)
 757                     : mapping.locateInFrom(begin, end);
 758             SequenceFeature sf2 = new SequenceFeature(sf, range[0],
 759                     range[1], group, sf.getScore());
 760             complementary.add(sf2);
 761           }
 762         }
 763       }
 764     }
 765
 766     return complementary;
 767   }
 768
 769   /**
 770    * Outputs any feature filters defined for visible feature types, sandwiched by
 771    * STARTFILTERS and ENDFILTERS lines
 772    *
 773    * @param out
 774    * @param visible
 775    * @param featureFilters
 776    */
 777   void outputFeatureFilters(StringBuilder out,
 778           Map<String, FeatureColourI> visible,
 779           Map<String, FeatureMatcherSetI> featureFilters)
 780   {
 781     if (visible == null || featureFilters == null
 782             || featureFilters.isEmpty())
 783     {
 784       return;
 785     }
 786
 787     boolean first = true;
 788     for (String featureType : visible.keySet())
 789     {
 790       FeatureMatcherSetI filter = featureFilters.get(featureType);
 791       if (filter != null)
 792       {
 793         if (first)
 794         {
 795           first = false;
 796           out.append(newline).append(STARTFILTERS).append(newline);
 797         }
 798         out.append(featureType).append(TAB).append(filter.toStableString())
 799                 .append(newline);
 800       }
 801     }
 802     if (!first)
 803     {
 804       out.append(ENDFILTERS).append(newline);
 805     }
 806
 807   }
 808
 809   /**
 810    * Appends output of visible sequence features within feature groups to the
 811    * output buffer. Groups other than the null or empty group are sandwiched by
 812    * STARTGROUP and ENDGROUP lines. Answers the number of features written.
 813    *
 814    * @param out
 815    * @param fr
 816    * @param featureTypes
 817    * @param sequences
 818    * @param includeNonPositional
 819    * @return
 820    */
 821   private int outputFeaturesByGroup(StringBuilder out,
 822           FeatureRenderer fr, String[] featureTypes,
 823           SequenceI[] sequences, boolean includeNonPositional)
 824   {
 825     List<String> featureGroups = fr.getFeatureGroups();
 826
 827     /*
 828      * sort groups alphabetically, and ensure that features with a
 829      * null or empty group are output after those in named groups
 830      */
 831     List<String> sortedGroups = new ArrayList<>(featureGroups);
 832     sortedGroups.remove(null);
 833     sortedGroups.remove("");
 834     Collections.sort(sortedGroups);
 835     sortedGroups.add(null);
 836     sortedGroups.add("");
 837
 838     int count = 0;
 839     List<String> visibleGroups = fr.getDisplayedFeatureGroups();
 840
 841     /*
 842      * loop over all groups (may be visible or not);
 843      * non-positional features are output even if group is not visible
 844      */
 845     for (String group : sortedGroups)
 846     {
 847       boolean firstInGroup = true;
 848       boolean isNullGroup = group == null || "".equals(group);
 849
 850       for (int i = 0; i < sequences.length; i++)
 851       {
 852         String sequenceName = sequences[i].getName();
 853         List<SequenceFeature> features = new ArrayList<>();
 854
 855         /*
 856          * get any non-positional features in this group, if wanted
 857          * (for any feature type, whether visible or not)
 858          */
 859         if (includeNonPositional)
 860         {
 861           features.addAll(sequences[i].getFeatures()
 862                   .getFeaturesForGroup(false, group));
 863         }
 864
 865         /*
 866          * add positional features for visible feature types, but
 867          * (for named groups) only if feature group is visible
 868          */
 869         if (featureTypes.length > 0
 870                 && (isNullGroup || visibleGroups.contains(group)))
 871         {
 872           features.addAll(sequences[i].getFeatures().getFeaturesForGroup(
 873                   true, group, featureTypes));
 874         }
 875
 876         for (SequenceFeature sf : features)
 877         {
 878           if (sf.isNonPositional() || fr.isVisible(sf))
 879           {
 880             count++;
 881             if (firstInGroup)
 882             {
 883               out.append(newline);
 884               if (!isNullGroup)
 885               {
 886                 out.append(STARTGROUP).append(TAB).append(group)
 887                         .append(newline);
 888               }
 889             }
 890             firstInGroup = false;
 891             formatJalviewFeature(out, sequenceName, sf);
 892           }
 893         }
 894       }
 895
 896       if (!isNullGroup && !firstInGroup)
 897       {
 898         out.append(ENDGROUP).append(TAB).append(group).append(newline);
 899       }
 900     }
 901     return count;
 902   }
 903
 904   /**
 905    * Formats one feature in Jalview format and appends to the string buffer
 906    *
 907    * @param out
 908    * @param sequenceName
 909    * @param sequenceFeature
 910    */
 911   protected void formatJalviewFeature(
 912           StringBuilder out, String sequenceName,
 913           SequenceFeature sequenceFeature)
 914   {
 915     if (sequenceFeature.description == null
 916             || sequenceFeature.description.equals(""))
 917     {
 918       out.append(sequenceFeature.type).append(TAB);
 919     }
 920     else
 921     {
 922       if (sequenceFeature.links != null
 923               && sequenceFeature.getDescription().indexOf("<html>") == -1)
 924       {
 925         out.append("<html>");
 926       }
 927
 928       out.append(sequenceFeature.description);
 929       if (sequenceFeature.links != null)
 930       {
 931         for (int l = 0; l < sequenceFeature.links.size(); l++)
 932         {
 933           String label = sequenceFeature.links.elementAt(l);
 934           String href = label.substring(label.indexOf("|") + 1);
 935           label = label.substring(0, label.indexOf("|"));
 936
 937           if (sequenceFeature.description.indexOf(href) == -1)
 938           {
 939             out.append(" <a href=\"").append(href).append("\">")
 940                     .append(label).append("</a>");
 941           }
 942         }
 943
 944         if (sequenceFeature.getDescription().indexOf("</html>") == -1)
 945         {
 946           out.append("</html>");
 947         }
 948       }
 949
 950       out.append(TAB);
 951     }
 952     out.append(sequenceName);
 953     out.append("\t-1\t");
 954     out.append(sequenceFeature.begin);
 955     out.append(TAB);
 956     out.append(sequenceFeature.end);
 957     out.append(TAB);
 958     out.append(sequenceFeature.type);
 959     if (!Float.isNaN(sequenceFeature.score))
 960     {
 961       out.append(TAB);
 962       out.append(sequenceFeature.score);
 963     }
 964     out.append(newline);
 965   }
 966
 967   /**
 968    * Parse method that is called when a GFF file is dragged to the desktop
 969    */
 970   @Override
 971   public void parse()
 972   {
 973     AlignViewportI av = getViewport();
 974     if (av != null)
 975     {
 976       if (av.getAlignment() != null)
 977       {
 978         dataset = av.getAlignment().getDataset();
 979       }
 980       if (dataset == null)
 981       {
 982         // working in the applet context ?
 983         dataset = av.getAlignment();
 984       }
 985     }
 986     else
 987     {
 988       dataset = new Alignment(new SequenceI[] {});
 989     }
 990
 991     Map<String, FeatureColourI> featureColours = new HashMap<>();
 992     boolean parseResult = parse(dataset, featureColours, false, true);
 993     if (!parseResult)
 994     {
 995       // pass error up somehow
 996     }
 997     if (av != null)
 998     {
 999       // update viewport with the dataset data ?
1000     }
1001     else
1002     {
1003       setSeqs(dataset.getSequencesArray());
1004     }
1005   }
1006
1007   /**
1008    * Implementation of unused abstract method
1009    *
1010    * @return error message
1011    */
1012   @Override
1013   public String print(SequenceI[] sqs, boolean jvsuffix)
1014   {
1015     System.out.println("Use printGffFormat() or printJalviewFormat()");
1016     return null;
1017   }
1018
1019   /**
1020    * Returns features output in GFF2 format
1021    *
1022    * @param sequences
1023    *                                       the sequences whose features are to be
1024    *                                       output
1025    * @param visible
1026    *                                       a map whose keys are the type names of
1027    *                                       visible features
1028    * @param visibleFeatureGroups
1029    * @param includeNonPositionalFeatures
1030    * @param includeComplement
1031    * @return
1032    */
1033   public String printGffFormat(SequenceI[] sequences,
1034           FeatureRenderer fr, boolean includeNonPositionalFeatures,
1035           boolean includeComplement)
1036   {
1037     FeatureRenderer fr2 = null;
1038     if (includeComplement)
1039     {
1040       AlignViewportI comp = fr.getViewport().getCodingComplement();
1041       fr2 = Desktop.getAlignFrameFor(comp).getFeatureRenderer();
1042     }
1043
1044     Map<String, FeatureColourI> visibleColours = fr.getDisplayedFeatureCols();
1045
1046     StringBuilder out = new StringBuilder(256);
1047
1048     out.append(String.format("%s %d\n", GFF_VERSION, gffVersion == 0 ? 2 : gffVersion));
1049
1050     String[] types = visibleColours == null ? new String[0]
1051             : visibleColours.keySet()
1052                     .toArray(new String[visibleColours.keySet().size()]);
1053
1054     for (SequenceI seq : sequences)
1055     {
1056       List<SequenceFeature> seqFeatures = new ArrayList<>();
1057       List<SequenceFeature> features = new ArrayList<>();
1058       if (includeNonPositionalFeatures)
1059       {
1060         features.addAll(seq.getFeatures().getNonPositionalFeatures());
1061       }
1062       if (visibleColours != null && !visibleColours.isEmpty())
1063       {
1064         features.addAll(seq.getFeatures().getPositionalFeatures(types));
1065       }
1066       for (SequenceFeature sf : features)
1067       {
1068         if (sf.isNonPositional() || fr.isVisible(sf))
1069         {
1070           /*
1071            * drop features hidden by group visibility, colour threshold,
1072            * or feature filter condition
1073            */
1074           seqFeatures.add(sf);
1075         }
1076       }
1077
1078       if (includeComplement)
1079       {
1080         seqFeatures.addAll(findComplementaryFeatures(seq, fr2));
1081       }
1082
1083       /*
1084        * sort features here if wanted
1085        */
1086       for (SequenceFeature sf : seqFeatures)
1087       {
1088         formatGffFeature(out, seq, sf);
1089         out.append(newline);
1090       }
1091     }
1092
1093     return out.toString();
1094   }
1095
1096   /**
1097    * Formats one feature as GFF and appends to the string buffer
1098    */
1099   private void formatGffFeature(StringBuilder out, SequenceI seq,
1100           SequenceFeature sf)
1101   {
1102     String source = sf.featureGroup;
1103     if (source == null)
1104     {
1105       source = sf.getDescription();
1106     }
1107
1108     out.append(seq.getName());
1109     out.append(TAB);
1110     out.append(source);
1111     out.append(TAB);
1112     out.append(sf.type);
1113     out.append(TAB);
1114     out.append(sf.begin);
1115     out.append(TAB);
1116     out.append(sf.end);
1117     out.append(TAB);
1118     out.append(sf.score);
1119     out.append(TAB);
1120
1121     int strand = sf.getStrand();
1122     out.append(strand == 1 ? "+" : (strand == -1 ? "-" : "."));
1123     out.append(TAB);
1124
1125     String phase = sf.getPhase();
1126     out.append(phase == null ? "." : phase);
1127
1128     if (sf.otherDetails != null && !sf.otherDetails.isEmpty())
1129     {
1130       Map<String, Object> map = sf.otherDetails;
1131       formatAttributes(out, map);
1132     }
1133   }
1134
1135   /**
1136    * A helper method that outputs attributes stored in the map as
1137    * semicolon-delimited values e.g.
1138    *
1139    * <pre>
1140    * AC_Male=0;AF_NFE=0.00000e 00;Hom_FIN=0;GQ_MEDIAN=9
1141    * </pre>
1142    *
1143    * A map-valued attribute is formatted as a comma-delimited list within braces,
1144    * for example
1145    *
1146    * <pre>
1147    * jvmap_CSQ={ALLELE_NUM=1,UNIPARC=UPI0002841053,Feature=ENST00000585561}
1148    * </pre>
1149    *
1150    * The {@code jvmap_} prefix designates a values map and is removed if the value
1151    * is parsed when read in. (The GFF3 specification allows 'semi-structured data'
1152    * to be represented provided the attribute name begins with a lower case
1153    * letter.)
1154    *
1155    * @param sb
1156    * @param map
1157    * @see http://gmod.org/wiki/GFF3#GFF3_Format
1158    */
1159   void formatAttributes(StringBuilder sb, Map<String, Object> map)
1160   {
1161     sb.append(TAB);
1162     boolean first = true;
1163     for (String key : map.keySet())
1164     {
1165       if (SequenceFeature.STRAND.equals(key)
1166               || SequenceFeature.PHASE.equals(key))
1167       {
1168         /*
1169          * values stashed in map but output to their own columns
1170          */
1171         continue;
1172       }
1173       {
1174         if (!first)
1175         {
1176           sb.append(";");
1177         }
1178       }
1179       first = false;
1180       Object value = map.get(key);
1181       if (value instanceof Map<?, ?>)
1182       {
1183         formatMapAttribute(sb, key, (Map<?, ?>) value);
1184       }
1185       else
1186       {
1187         String formatted = StringUtils.urlEncode(value.toString(),
1188                 GffHelperI.GFF_ENCODABLE);
1189         sb.append(key).append(EQUALS).append(formatted);
1190       }
1191     }
1192   }
1193
1194   /**
1195    * Formats the map entries as
1196    *
1197    * <pre>
1198    * key=key1=value1,key2=value2,...
1199    * </pre>
1200    *
1201    * and appends this to the string buffer
1202    *
1203    * @param sb
1204    * @param key
1205    * @param map
1206    */
1207   private void formatMapAttribute(StringBuilder sb, String key,
1208           Map<?, ?> map)
1209   {
1210     if (map == null || map.isEmpty())
1211     {
1212       return;
1213     }
1214
1215     /*
1216      * AbstractMap.toString would be a shortcut here, but more reliable
1217      * to code the required format in case toString changes in future
1218      */
1219     sb.append(key).append(EQUALS);
1220     boolean first = true;
1221     for (Entry<?, ?> entry : map.entrySet())
1222     {
1223       if (!first)
1224       {
1225         sb.append(",");
1226       }
1227       first = false;
1228       sb.append(entry.getKey().toString()).append(EQUALS);
1229       String formatted = StringUtils.urlEncode(entry.getValue().toString(),
1230               GffHelperI.GFF_ENCODABLE);
1231       sb.append(formatted);
1232     }
1233   }
1234
1235   /**
1236    * Returns a mapping given list of one or more Align descriptors (exonerate
1237    * format)
1238    *
1239    * @param alignedRegions
1240    *                         a list of "Align fromStart toStart fromCount"
1241    * @param mapIsFromCdna
1242    *                         if true, 'from' is dna, else 'from' is protein
1243    * @param strand
1244    *                         either 1 (forward) or -1 (reverse)
1245    * @return
1246    * @throws IOException
1247    */
1248   protected MapList constructCodonMappingFromAlign(
1249           List<String> alignedRegions, boolean mapIsFromCdna, int strand)
1250           throws IOException
1251   {
1252     if (strand == 0)
1253     {
1254       throw new IOException(
1255               "Invalid strand for a codon mapping (cannot be 0)");
1256     }
1257     int regions = alignedRegions.size();
1258     // arrays to hold [start, end] for each aligned region
1259     int[] fromRanges = new int[regions * 2]; // from dna
1260     int[] toRanges = new int[regions * 2]; // to protein
1261     int fromRangesIndex = 0;
1262     int toRangesIndex = 0;
1263
1264     for (String range : alignedRegions)
1265     {
1266       /*
1267        * Align mapFromStart mapToStart mapFromCount
1268        * e.g. if mapIsFromCdna
1269        *     Align 11270 143 120
1270        * means:
1271        *     120 bases from pos 11270 align to pos 143 in peptide
1272        * if !mapIsFromCdna this would instead be
1273        *     Align 143 11270 40
1274        */
1275       String[] tokens = range.split(" ");
1276       if (tokens.length != 3)
1277       {
1278         throw new IOException("Wrong number of fields for Align");
1279       }
1280       int fromStart = 0;
1281       int toStart = 0;
1282       int fromCount = 0;
1283       try
1284       {
1285         fromStart = Integer.parseInt(tokens[0]);
1286         toStart = Integer.parseInt(tokens[1]);
1287         fromCount = Integer.parseInt(tokens[2]);
1288       } catch (NumberFormatException nfe)
1289       {
1290         throw new IOException(
1291                 "Invalid number in Align field: " + nfe.getMessage());
1292       }
1293
1294       /*
1295        * Jalview always models from dna to protein, so adjust values if the
1296        * GFF mapping is from protein to dna
1297        */
1298       if (!mapIsFromCdna)
1299       {
1300         fromCount *= 3;
1301         int temp = fromStart;
1302         fromStart = toStart;
1303         toStart = temp;
1304       }
1305       fromRanges[fromRangesIndex++] = fromStart;
1306       fromRanges[fromRangesIndex++] = fromStart + strand * (fromCount - 1);
1307
1308       /*
1309        * If a codon has an intron gap, there will be contiguous 'toRanges';
1310        * this is handled for us by the MapList constructor.
1311        * (It is not clear that exonerate ever generates this case)
1312        */
1313       toRanges[toRangesIndex++] = toStart;
1314       toRanges[toRangesIndex++] = toStart + (fromCount - 1) / 3;
1315     }
1316
1317     return new MapList(fromRanges, toRanges, 3, 1);
1318   }
1319
1320   /**
1321    * Parse a GFF format feature. This may include creating a 'dummy' sequence to
1322    * hold the feature, or for its mapped sequence, or both, to be resolved
1323    * either later in the GFF file (##FASTA section), or when the user loads
1324    * additional sequences.
1325    *
1326    * @param gffColumns
1327    * @param alignment
1328    * @param relaxedIdMatching
1329    * @param newseqs
1330    * @return
1331    */
1332   protected SequenceI parseGff(String[] gffColumns, AlignmentI alignment,
1333           boolean relaxedIdMatching, List<SequenceI> newseqs)
1334   {
1335     /*
1336      * GFF: seqid source type start end score strand phase [attributes]
1337      */
1338     if (gffColumns.length < 5)
1339     {
1340       System.err.println("Ignoring GFF feature line with too few columns ("
1341               + gffColumns.length + ")");
1342       return null;
1343     }
1344
1345     /*
1346      * locate referenced sequence in alignment _or_
1347      * as a forward or external reference (SequenceDummy)
1348      */
1349     String seqId = gffColumns[0];
1350     SequenceI seq = findSequence(seqId, alignment, newseqs,
1351             relaxedIdMatching);
1352
1353     SequenceFeature sf = null;
1354     GffHelperI helper = GffHelperFactory.getHelper(gffColumns);
1355     if (helper != null)
1356     {
1357       try
1358       {
1359         sf = helper.processGff(seq, gffColumns, alignment, newseqs,
1360                 relaxedIdMatching);
1361         if (sf != null)
1362         {
1363           seq.addSequenceFeature(sf);
1364           while ((seq = alignment.findName(seq, seqId, true)) != null)
1365           {
1366             seq.addSequenceFeature(new SequenceFeature(sf));
1367           }
1368         }
1369       } catch (IOException e)
1370       {
1371         System.err.println("GFF parsing failed with: " + e.getMessage());
1372         return null;
1373       }
1374     }
1375
1376     return seq;
1377   }
1378
1379   /**
1380    * After encountering ##fasta in a GFF3 file, process the remainder of the
1381    * file as FAST sequence data. Any placeholder sequences created during
1382    * feature parsing are updated with the actual sequences.
1383    *
1384    * @param align
1385    * @param newseqs
1386    * @throws IOException
1387    */
1388   protected void processAsFasta(AlignmentI align, List<SequenceI> newseqs)
1389           throws IOException
1390   {
1391     try
1392     {
1393       mark();
1394     } catch (IOException q)
1395     {
1396     }
1397     FastaFile parser = new FastaFile(this);
1398     List<SequenceI> includedseqs = parser.getSeqs();
1399
1400     SequenceIdMatcher smatcher = new SequenceIdMatcher(newseqs);
1401
1402     /*
1403      * iterate over includedseqs, and replacing matching ones with newseqs
1404      * sequences. Generic iterator not used here because we modify
1405      * includedseqs as we go
1406      */
1407     for (int p = 0, pSize = includedseqs.size(); p < pSize; p++)
1408     {
1409       // search for any dummy seqs that this sequence can be used to update
1410       SequenceI includedSeq = includedseqs.get(p);
1411       SequenceI dummyseq = smatcher.findIdMatch(includedSeq);
1412       if (dummyseq != null && dummyseq instanceof SequenceDummy)
1413       {
1414         // probably have the pattern wrong
1415         // idea is that a flyweight proxy for a sequence ID can be created for
1416         // 1. stable reference creation
1417         // 2. addition of annotation
1418         // 3. future replacement by a real sequence
1419         // current pattern is to create SequenceDummy objects - a convenience
1420         // constructor for a Sequence.
1421         // problem is that when promoted to a real sequence, all references
1422         // need to be updated somehow. We avoid that by keeping the same object.
1423         ((SequenceDummy) dummyseq).become(includedSeq);
1424         dummyseq.createDatasetSequence();
1425
1426         /*
1427          * Update mappings so they are now to the dataset sequence
1428          */
1429         for (AlignedCodonFrame mapping : align.getCodonFrames())
1430         {
1431           mapping.updateToDataset(dummyseq);
1432         }
1433
1434         /*
1435          * replace parsed sequence with the realised forward reference
1436          */
1437         includedseqs.set(p, dummyseq);
1438
1439         /*
1440          * and remove from the newseqs list
1441          */
1442         newseqs.remove(dummyseq);
1443       }
1444     }
1445
1446     /*
1447      * finally add sequences to the dataset
1448      */
1449     for (SequenceI seq : includedseqs)
1450     {
1451       // experimental: mapping-based 'alignment' to query sequence
1452       AlignmentUtils.alignSequenceAs(seq, align,
1453               String.valueOf(align.getGapCharacter()), false, true);
1454
1455       // rename sequences if GFF handler requested this
1456       // TODO a more elegant way e.g. gffHelper.postProcess(newseqs) ?
1457       List<SequenceFeature> sfs = seq.getFeatures().getPositionalFeatures();
1458       if (!sfs.isEmpty())
1459       {
1460         String newName = (String) sfs.get(0).getValue(
1461                 GffHelperI.RENAME_TOKEN);
1462         if (newName != null)
1463         {
1464           seq.setName(newName);
1465         }
1466       }
1467       align.addSequence(seq);
1468     }
1469   }
1470
1471   /**
1472    * Process a ## directive
1473    *
1474    * @param line
1475    * @param gffProps
1476    * @param align
1477    * @param newseqs
1478    * @throws IOException
1479    */
1480   protected void processGffPragma(String line, Map<String, String> gffProps,
1481           AlignmentI align, List<SequenceI> newseqs) throws IOException
1482   {
1483     line = line.trim();
1484     if ("###".equals(line))
1485     {
1486       // close off any open 'forward references'
1487       return;
1488     }
1489
1490     String[] tokens = line.substring(2).split(" ");
1491     String pragma = tokens[0];
1492     String value = tokens.length == 1 ? null : tokens[1];
1493
1494     if ("gff-version".equalsIgnoreCase(pragma))
1495     {
1496       if (value != null)
1497       {
1498         try
1499         {
1500           // value may be e.g. "3.1.2"
1501           gffVersion = Integer.parseInt(value.split("\\.")[0]);
1502         } catch (NumberFormatException e)
1503         {
1504           // ignore
1505         }
1506       }
1507     }
1508     else if ("sequence-region".equalsIgnoreCase(pragma))
1509     {
1510       // could capture <seqid start end> if wanted here
1511     }
1512     else if ("feature-ontology".equalsIgnoreCase(pragma))
1513     {
1514       // should resolve against the specified feature ontology URI
1515     }
1516     else if ("attribute-ontology".equalsIgnoreCase(pragma))
1517     {
1518       // URI of attribute ontology - not currently used in GFF3
1519     }
1520     else if ("source-ontology".equalsIgnoreCase(pragma))
1521     {
1522       // URI of source ontology - not currently used in GFF3
1523     }
1524     else if ("species-build".equalsIgnoreCase(pragma))
1525     {
1526       // save URI of specific NCBI taxon version of annotations
1527       gffProps.put("species-build", value);
1528     }
1529     else if ("fasta".equalsIgnoreCase(pragma))
1530     {
1531       // process the rest of the file as a fasta file and replace any dummy
1532       // sequence IDs
1533       processAsFasta(align, newseqs);
1534     }
1535     else
1536     {
1537       System.err.println("Ignoring unknown pragma: " + line);
1538     }
1539   }
1540 }