src/jalview/io/gff/ExonerateHelper.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.gff;
  22
  23 import jalview.datamodel.AlignedCodonFrame;
  24 import jalview.datamodel.AlignmentI;
  25 import jalview.datamodel.MappingType;
  26 import jalview.datamodel.SequenceFeature;
  27 import jalview.datamodel.SequenceI;
  28 import jalview.util.MapList;
  29
  30 import java.io.IOException;
  31 import java.util.List;
  32 import java.util.Map;
  33
  34 /**
  35  * A handler to parse GFF in the format generated by the exonerate tool
  36  */
  37 public class ExonerateHelper extends Gff2Helper
  38 {
  39   private static final String SIMILARITY = "similarity";
  40
  41   private static final String GENOME2GENOME = "genome2genome";
  42
  43   private static final String CDNA2GENOME = "cdna2genome";
  44
  45   private static final String CODING2GENOME = "coding2genome";
  46
  47   private static final String CODING2CODING = "coding2coding";
  48
  49   private static final String PROTEIN2GENOME = "protein2genome";
  50
  51   private static final String PROTEIN2DNA = "protein2dna";
  52
  53   private static final String ALIGN = "Align";
  54
  55   private static final String QUERY = "Query";
  56
  57   private static final String TARGET = "Target";
  58
  59   /**
  60    * Process one GFF feature line (as modelled by SequenceFeature)
  61    *
  62    * @param seq
  63    *          the sequence with which this feature is associated
  64    * @param gffColumns
  65    *          the sequence feature with ATTRIBUTES property containing any
  66    *          additional attributes
  67    * @param align
  68    *          the alignment we are adding GFF to
  69    * @param newseqs
  70    *          any new sequences referenced by the GFF
  71    * @param relaxedIdMatching
  72    *          if true, match word tokens in sequence names
  73    * @return true if the sequence feature should be added to the sequence, else
  74    *         false (i.e. it has been processed in another way e.g. to generate a
  75    *         mapping)
  76    */
  77   @Override
  78   public SequenceFeature processGff(SequenceI seq, String[] gffColumns,
  79           AlignmentI align, List<SequenceI> newseqs,
  80           boolean relaxedIdMatching)
  81   {
  82     String attr = gffColumns[ATTRIBUTES_COL];
  83     Map<String, List<String>> set = parseNameValuePairs(attr);
  84
  85     try
  86     {
  87       processGffSimilarity(set, seq, gffColumns, align, newseqs,
  88               relaxedIdMatching);
  89     } catch (IOException ivfe)
  90     {
  91       System.err.println(ivfe);
  92     }
  93
  94     /*
  95      * return null to indicate we don't want to add a sequence feature for
  96      * similarity (only process it to create mappings)
  97      */
  98     return null;
  99   }
 100
 101   /**
 102    * Processes the 'Query' (or 'Target') and 'Align' properties associated with
 103    * an exonerate GFF similarity feature; these properties define the mapping of
 104    * the annotated range to a related sequence.
 105    *
 106    * @param set
 107    *          parsed GFF column 9 key/value(s)
 108    * @param seq
 109    *          the sequence the GFF feature is on
 110    * @param gff
 111    *          the GFF column data
 112    * @param align
 113    *          the alignment the sequence belongs to, where any new mappings
 114    *          should be added
 115    * @param newseqs
 116    *          a list of new 'virtual sequences' generated while parsing GFF
 117    * @param relaxedIdMatching
 118    *          if true allow fuzzy search for a matching target sequence
 119    * @throws IOException
 120    */
 121   protected void processGffSimilarity(Map<String, List<String>> set,
 122           SequenceI seq, String[] gff, AlignmentI align,
 123           List<SequenceI> newseqs, boolean relaxedIdMatching)
 124           throws IOException
 125   {
 126     /*
 127      * exonerate may be run with
 128      * --showquerygff - outputs 'features on the query' e.g. (protein2genome)
 129      *     Target <dnaseqid> ; Align proteinStartPos dnaStartPos proteinCount
 130      * --showtargetgff - outputs 'features on the target' e.g. (protein2genome)
 131      *     Query <proteinseqid> ; Align dnaStartPos proteinStartPos nucleotideCount
 132      * where the Align spec may repeat
 133      */
 134     // TODO handle coding2coding and similar as well
 135     boolean featureIsOnTarget = true;
 136     List<String> mapTo = set.get(QUERY);
 137     if (mapTo == null)
 138     {
 139       mapTo = set.get(TARGET);
 140       featureIsOnTarget = false;
 141     }
 142     MappingType type = getMappingType(gff[SOURCE_COL]);
 143
 144     if (type == null)
 145     {
 146       throw new IOException("Sorry, I don't handle " + gff[SOURCE_COL]);
 147     }
 148
 149     if (mapTo == null || mapTo.size() != 1)
 150     {
 151       throw new IOException(
 152               "Expecting exactly one sequence in Query or Target field (got "
 153                       + mapTo + ")");
 154     }
 155
 156     /*
 157      * locate the mapped sequence in the alignment or 'new' (GFF file) sequences;
 158      */
 159     SequenceI mappedSequence = findSequence(mapTo.get(0), align, newseqs,
 160             relaxedIdMatching);
 161
 162     /*
 163      * If mapping is from protein to dna, we store it as dna to protein instead
 164      */
 165     SequenceI mapFromSequence = seq;
 166     SequenceI mapToSequence = mappedSequence;
 167     if ((MappingType.isDnaToPeptide(type) && featureIsOnTarget)
 168             || (MappingType.isPeptideToDna(type) && !featureIsOnTarget))
 169     {
 170       mapFromSequence = mappedSequence;
 171       mapToSequence = seq;
 172     }
 173
 174     /*
 175      * Process the Align maps and create mappings.
 176      * These may be cdna-genome, cdna-protein, genome-protein.
 177      * The mapped sequences may or may not be in the alignment
 178      * (they may be included later in the GFF file).
 179      */
 180
 181     /*
 182      * get any existing mapping for these sequences (or start one),
 183      * and add this mapped range
 184      */
 185     AlignedCodonFrame acf = getMapping(align, mapFromSequence,
 186             mapToSequence);
 187
 188     /*
 189      * exonerate GFF has the strand of the target in column 7
 190      * (differs from GFF3 which has it in the Target descriptor)
 191      */
 192     String strand = gff[STRAND_COL];
 193     boolean forwardStrand = true;
 194     if ("-".equals(strand))
 195     {
 196       forwardStrand = false;
 197     }
 198     else if (!"+".equals(strand))
 199     {
 200       System.err.println("Strand must be specified for alignment");
 201       return;
 202     }
 203
 204     List<String> alignedRegions = set.get(ALIGN);
 205     for (String region : alignedRegions)
 206     {
 207       MapList mapping = buildMapping(region, type, forwardStrand,
 208               featureIsOnTarget, gff);
 209
 210       if (mapping == null)
 211       {
 212         continue;
 213       }
 214
 215       acf.addMap(mapFromSequence, mapToSequence, mapping);
 216     }
 217     align.addCodonFrame(acf);
 218   }
 219
 220   /**
 221    * Construct the mapping
 222    *
 223    * @param region
 224    * @param type
 225    * @param forwardStrand
 226    * @param featureIsOnTarget
 227    * @param gff
 228    * @return
 229    */
 230   protected MapList buildMapping(String region, MappingType type,
 231           boolean forwardStrand, boolean featureIsOnTarget, String[] gff)
 232   {
 233     /*
 234      * process one "fromStart toStart fromCount" descriptor
 235      */
 236     String[] tokens = region.split(" ");
 237     if (tokens.length != 3)
 238     {
 239       System.err.println("Malformed Align descriptor: " + region);
 240       return null;
 241     }
 242
 243     /*
 244      * get start/end of from/to mappings
 245      * if feature is on the target sequence we have to invert the sense
 246      */
 247     int alignFromStart;
 248     int alignToStart;
 249     int alignCount;
 250     try
 251     {
 252       alignFromStart = Integer.parseInt(tokens[0]);
 253       alignToStart = Integer.parseInt(tokens[1]);
 254       alignCount = Integer.parseInt(tokens[2]);
 255     } catch (NumberFormatException nfe)
 256     {
 257       System.err.println(nfe.toString());
 258       return null;
 259     }
 260
 261     int fromStart;
 262     int fromEnd;
 263     int toStart;
 264     int toEnd;
 265
 266     if (featureIsOnTarget)
 267     {
 268       fromStart = alignToStart;
 269       toStart = alignFromStart;
 270       toEnd = forwardStrand ? toStart + alignCount - 1
 271               : toStart - (alignCount - 1);
 272       int toLength = Math.abs(toEnd - toStart) + 1;
 273       int fromLength = toLength * type.getFromRatio() / type.getToRatio();
 274       fromEnd = fromStart + fromLength - 1;
 275     }
 276     else
 277     {
 278       // we use the 'Align' values here not the feature start/end
 279       // not clear why they may differ but it seems they can
 280       fromStart = alignFromStart;
 281       fromEnd = alignFromStart + alignCount - 1;
 282       int fromLength = fromEnd - fromStart + 1;
 283       int toLength = fromLength * type.getToRatio() / type.getFromRatio();
 284       toStart = alignToStart;
 285       if (forwardStrand)
 286       {
 287         toEnd = toStart + toLength - 1;
 288       }
 289       else
 290       {
 291         toEnd = toStart - (toLength - 1);
 292       }
 293     }
 294
 295     MapList codonmapping = constructMappingFromAlign(fromStart, fromEnd,
 296             toStart, toEnd, type);
 297     return codonmapping;
 298   }
 299
 300   /**
 301    * Returns a MappingType depending on the exonerate 'model' value
 302    *
 303    * @param model
 304    * @return
 305    */
 306   protected static MappingType getMappingType(String model)
 307   {
 308     MappingType result = null;
 309
 310     if (model.contains(PROTEIN2GENOME))
 311     {
 312       result = MappingType.PeptideToGenome;
 313     }
 314     else if (model.contains(PROTEIN2DNA))
 315     {
 316       result = MappingType.PeptideToCds;
 317     }
 318     else if (model.contains(CODING2GENOME))
 319     {
 320       result = MappingType.CdsToGenome;
 321     }
 322     else if (model.contains(CDNA2GENOME))
 323     {
 324       result = MappingType.CdnaToGenome;
 325     }
 326     else if (model.contains(CODING2CODING) || model.contains(GENOME2GENOME))
 327     {
 328       result = MappingType.GenomeToCdna;
 329     }
 330     return result;
 331   }
 332
 333   /**
 334    * Tests whether the GFF data looks like it was generated by exonerate, and is
 335    * a format we are willing to handle
 336    *
 337    * @param columns
 338    * @return
 339    */
 340   public static boolean recognises(String[] columns)
 341   {
 342     if (!SIMILARITY.equalsIgnoreCase(columns[TYPE_COL]))
 343     {
 344       return false;
 345     }
 346
 347     /*
 348      * inspect alignment model
 349      */
 350     String model = columns[SOURCE_COL];
 351     // e.g. exonerate:protein2genome:local
 352     if (model != null)
 353     {
 354       String mdl = model.toLowerCase();
 355       if (mdl.contains(PROTEIN2DNA) || mdl.contains(PROTEIN2GENOME)
 356               || mdl.contains(CODING2CODING) || mdl.contains(CODING2GENOME)
 357               || mdl.contains(CDNA2GENOME) || mdl.contains(GENOME2GENOME))
 358       {
 359         return true;
 360       }
 361     }
 362     System.err.println("Sorry, I don't handle exonerate model " + model);
 363     return false;
 364   }
 365
 366   /**
 367    * An override to set feature group to "exonerate" instead of the default GFF
 368    * source value (column 2)
 369    */
 370   @Override
 371   protected SequenceFeature buildSequenceFeature(String[] gff,
 372           Map<String, List<String>> set)
 373   {
 374     SequenceFeature sf = super.buildSequenceFeature(gff, TYPE_COL,
 375             "exonerate", set);
 376
 377     return sf;
 378   }
 379
 380 }