src/jalview/datamodel/MappedFeatures.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.datamodel;
  22
  23 import java.util.Locale;
  24
  25 import java.util.HashSet;
  26 import java.util.List;
  27 import java.util.Set;
  28
  29 import jalview.io.gff.Gff3Helper;
  30 import jalview.schemes.ResidueProperties;
  31 import jalview.util.MapList;
  32 import jalview.util.MappingUtils;
  33 import jalview.util.StringUtils;
  34
  35 /**
  36  * A data bean to hold a list of mapped sequence features (e.g. CDS features
  37  * mapped from protein), and the mapping between the sequences. It also provides
  38  * a method to derive peptide variants from codon variants.
  39  *
  40  * @author gmcarstairs
  41  */
  42 public class MappedFeatures
  43 {
  44   /*
  45    * VEP CSQ:HGVSp (if present) is a short-cut to the protein variant consequence
  46    */
  47   private static final String HGV_SP = "HGVSp";
  48
  49   private static final String CSQ = "CSQ";
  50
  51   /*
  52    * the sequence the mapped features are on
  53    */
  54   private final SequenceI featureSequence;
  55
  56   /*
  57    * the mapping between sequences;
  58    * NB this could be in either sense (from or to featureSequence)
  59    */
  60   private final Mapping mapping;
  61
  62   /*
  63    * features on featureSequence that overlap the mapped positions
  64    */
  65   public final List<SequenceFeature> features;
  66
  67   /*
  68    * the residue position in the sequence mapped to
  69    */
  70   private final int toPosition;
  71
  72   /*
  73    * the residue at toPosition
  74    */
  75   private final char toResidue;
  76
  77   /*
  78    * if the mapping is 3:1 or 1:3 (peptide to CDS), this holds the
  79    * mapped positions i.e. codon base positions in CDS; to
  80    * support calculation of peptide variants from alleles
  81    */
  82   private final int[] codonPos;
  83
  84   private final char[] baseCodon;
  85
  86   /**
  87    * Constructor
  88    *
  89    * @param theMapping
  90    *          sequence mapping (which may be either to, or from, the sequence
  91    *          holding the linked features)
  92    * @param featureSeq
  93    *          the sequence hosting the virtual features
  94    * @param pos
  95    *          the residue position in the sequence mapped to
  96    * @param res
  97    *          the residue character at position pos
  98    * @param theFeatures
  99    *          list of mapped features found in the 'featureSeq' sequence at the
 100    *          mapped position(s)
 101    */
 102   public MappedFeatures(Mapping theMapping, SequenceI featureSeq, int pos,
 103           char res, List<SequenceFeature> theFeatures)
 104   {
 105     mapping = theMapping;
 106     featureSequence = featureSeq;
 107     toPosition = pos;
 108     toResidue = res;
 109     features = theFeatures;
 110
 111     /*
 112      * determine codon positions and canonical codon
 113      * for a peptide-to-CDS mapping
 114      */
 115     int[] codonIntervals = mapping.getMap().locateInFrom(toPosition, toPosition);
 116     int[] codonPositions = codonIntervals == null ? null
 117             : MappingUtils.flattenRanges(codonIntervals);
 118     if (codonPositions != null && codonPositions.length == 3)
 119     {
 120       codonPos = codonPositions;
 121       baseCodon = new char[3];
 122       int cdsStart = featureSequence.getStart();
 123       baseCodon[0] = Character
 124               .toUpperCase(featureSequence.getCharAt(codonPos[0] - cdsStart));
 125       baseCodon[1] = Character
 126               .toUpperCase(featureSequence.getCharAt(codonPos[1] - cdsStart));
 127       baseCodon[2] = Character
 128               .toUpperCase(featureSequence.getCharAt(codonPos[2] - cdsStart));
 129     }
 130     else
 131     {
 132       codonPos = null;
 133       baseCodon = null;
 134     }
 135   }
 136
 137   /**
 138    * Computes and returns comma-delimited HGVS notation peptide variants derived
 139    * from codon allele variants. If no variants are found, answers an empty
 140    * string. The peptide variant is either simply read from the "CSQ:HGVSp"
 141    * attribute if present, else computed based on the "alleles" attribute if
 142    * present. If neither attribute is found, no variant (empty string) is
 143    * returned.
 144    *
 145    * @param sf
 146    *          a sequence feature (which must be one of those held in this
 147    *          object)
 148    * @return
 149    */
 150   public String findProteinVariants(SequenceFeature sf)
 151   {
 152     if (!features.contains(sf) || baseCodon == null)
 153     {
 154       return "";
 155     }
 156
 157     /*
 158      * VCF data may already contain the protein consequence
 159      */
 160     String hgvsp = sf.getValueAsString(CSQ, HGV_SP);
 161     if (hgvsp != null)
 162     {
 163       int colonPos = hgvsp.lastIndexOf(':');
 164       if (colonPos >= 0)
 165       {
 166         String var = hgvsp.substring(colonPos + 1);
 167         if (var.contains("p.")) // sanity check
 168         {
 169           return var;
 170         }
 171       }
 172     }
 173
 174     /*
 175      * otherwise, compute codon and peptide variant
 176      */
 177     int cdsPos = sf.getBegin();
 178     if (cdsPos != sf.getEnd())
 179     {
 180       // not handling multi-locus variant features
 181       return "";
 182     }
 183     if (cdsPos != codonPos[0] && cdsPos != codonPos[1]
 184             && cdsPos != codonPos[2])
 185     {
 186       // e.g. feature on intron within spliced codon!
 187       return "";
 188     }
 189
 190     String alls = (String) sf.getValue(Gff3Helper.ALLELES);
 191     if (alls == null)
 192     {
 193       return "";
 194     }
 195
 196     String from3 = StringUtils.toSentenceCase(
 197             ResidueProperties.aa2Triplet.get(String.valueOf(toResidue)));
 198
 199     /*
 200      * make a peptide variant for each SNP allele
 201      * e.g. C,G,T gives variants G and T for base C
 202      */
 203     Set<String> variantPeptides = new HashSet<>();
 204     String[] alleles = alls.toUpperCase(Locale.ROOT).split(",");
 205     StringBuilder vars = new StringBuilder();
 206
 207     for (String allele : alleles)
 208     {
 209       allele = allele.trim().toUpperCase(Locale.ROOT);
 210       if (allele.length() > 1 || "-".equals(allele))
 211       {
 212         continue; // multi-locus variant
 213       }
 214       char[] variantCodon = new char[3];
 215       variantCodon[0] = baseCodon[0];
 216       variantCodon[1] = baseCodon[1];
 217       variantCodon[2] = baseCodon[2];
 218
 219       /*
 220        * poke variant base into canonical codon;
 221        * ignore first 'allele' (canonical base)
 222        */
 223       final int i = cdsPos == codonPos[0] ? 0
 224               : (cdsPos == codonPos[1] ? 1 : 2);
 225       variantCodon[i] = allele.toUpperCase(Locale.ROOT).charAt(0);
 226       if (variantCodon[i] == baseCodon[i])
 227       {
 228         continue;
 229       }
 230       String codon = new String(variantCodon);
 231       String peptide = ResidueProperties.codonTranslate(codon);
 232       boolean synonymous = toResidue == peptide.charAt(0);
 233       StringBuilder var = new StringBuilder();
 234       if (synonymous)
 235       {
 236         /*
 237          * synonymous variant notation e.g. c.1062C>A(p.=)
 238          */
 239         var.append("c.").append(String.valueOf(cdsPos))
 240                 .append(String.valueOf(baseCodon[i])).append(">")
 241                 .append(String.valueOf(variantCodon[i]))
 242                 .append("(p.=)");
 243       }
 244       else
 245       {
 246         /*
 247          * missense variant notation e.g. p.Arg355Met
 248          */
 249         String to3 = ResidueProperties.STOP.equals(peptide) ? "Ter"
 250                 : StringUtils.toSentenceCase(
 251                         ResidueProperties.aa2Triplet.get(peptide));
 252         var.append("p.").append(from3).append(String.valueOf(toPosition))
 253                 .append(to3);
 254       }
 255       if (!variantPeptides.contains(peptide)) // duplicate consequence
 256       {
 257         variantPeptides.add(peptide);
 258         if (vars.length() > 0)
 259         {
 260           vars.append(",");
 261         }
 262         vars.append(var);
 263       }
 264     }
 265
 266     return vars.toString();
 267   }
 268
 269   /**
 270    * Answers the name of the linked sequence holding any mapped features
 271    *
 272    * @return
 273    */
 274   public String getLinkedSequenceName()
 275   {
 276     return featureSequence == null ? null : featureSequence.getName();
 277   }
 278
 279   /**
 280    * Answers the mapped ranges (as one or more [start, end] positions) which
 281    * correspond to the given [begin, end] range of the linked sequence.
 282    *
 283    * <pre>
 284    * Example: MappedFeatures with CDS features mapped to peptide
 285    * CDS/200-220 gtc aac TGa acGt att AAC tta
 286    * mapped to PEP/6-7 WN by mapping [206, 207, 210, 210, 215, 217] to [6, 7]
 287    * getMappedPositions(206, 206) should return [6, 6]
 288    * getMappedPositions(200, 214) should return [6, 6]
 289    * getMappedPositions(210, 215) should return [6, 7]
 290    * </pre>
 291    *
 292    * @param begin
 293    * @param end
 294    * @return
 295    */
 296   public int[] getMappedPositions(int begin, int end)
 297   {
 298     MapList map = mapping.getMap();
 299     return mapping.to == featureSequence ? map.getOverlapsInFrom(begin, end)
 300             : map.getOverlapsInTo(begin, end);
 301   }
 302
 303   /**
 304    * Answers true if the linked features are on coding sequence, false if on
 305    * peptide
 306    *
 307    * @return
 308    */
 309   public boolean isFromCds()
 310   {
 311     if (mapping.getMap().getFromRatio() == 3)
 312     {
 313       return mapping.to != featureSequence;
 314     }
 315     return mapping.to == featureSequence;
 316   }
 317 }