src/jalview/datamodel/MappedFeatures.java

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