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