src/jalview/datamodel/MappedFeatures.java

   1 package jalview.datamodel;
   2
   3 import jalview.io.gff.Gff3Helper;
   4 import jalview.schemes.ResidueProperties;
   5 import jalview.util.MappingUtils;
   6 import jalview.util.StringUtils;
   7
   8 import java.util.HashSet;
   9 import java.util.List;
  10 import java.util.Set;
  11
  12 /**
  13  * A data bean to hold a list of mapped sequence features (e.g. CDS features
  14  * mapped from protein), and the mapping between the sequences. It also provides
  15  * a method to derive peptide variants from codon variants.
  16  *
  17  * @author gmcarstairs
  18  */
  19 public class MappedFeatures
  20 {
  21   private static final String HGV_SP = "HGVSp";
  22
  23   private static final String CSQ = "CSQ";
  24
  25   /*
  26    * the mapping from one sequence to another
  27    */
  28   public final Mapping mapping;
  29
  30   /**
  31    * the sequence mapped from
  32    */
  33   public final SequenceI fromSeq;
  34
  35   /*
  36    * features on the sequence mapped to that overlap the mapped positions
  37    */
  38   public final List<SequenceFeature> features;
  39
  40   /*
  41    * the residue position in the sequence mapped to
  42    */
  43   private final int toPosition;
  44
  45   /*
  46    * the residue at toPosition
  47    */
  48   private final char toResidue;
  49
  50   /*
  51    * if the mapping is 3:1 or 1:3 (peptide to CDS), this holds the
  52    * mapped positions i.e. codon base positions in CDS; to
  53    * support calculation of peptide variants from alleles
  54    */
  55   private final int[] codonPos;
  56
  57   private final char[] baseCodon;
  58
  59   /**
  60    * Constructor
  61    *
  62    * @param theMapping
  63    * @param from
  64    *                      the sequence mapped from (e.g. CDS)
  65    * @param pos
  66    *                      the residue position in the sequence mapped to
  67    * @param res
  68    *                      the residue character at position pos
  69    * @param theFeatures
  70    *                      list of mapped features found in the 'from' sequence at
  71    *                      the mapped position(s)
  72    */
  73   public MappedFeatures(Mapping theMapping, SequenceI from, int pos,
  74           char res, List<SequenceFeature> theFeatures)
  75   {
  76     mapping = theMapping;
  77     fromSeq = from;
  78     toPosition = pos;
  79     toResidue = res;
  80     features = theFeatures;
  81
  82     /*
  83      * determine codon positions and canonical codon
  84      * for a peptide-to-CDS mapping
  85      */
  86     int[] codonIntervals = mapping.getMap().locateInFrom(toPosition, toPosition);
  87     if (codonIntervals != null)
  88     {
  89       codonPos = MappingUtils.flattenRanges(codonIntervals);
  90       if (codonPos.length == 3)
  91       {
  92         baseCodon = new char[3];
  93         int cdsStart = fromSeq.getStart();
  94         baseCodon[0] = fromSeq.getCharAt(codonPos[0] - cdsStart);
  95         baseCodon[1] = fromSeq.getCharAt(codonPos[1] - cdsStart);
  96         baseCodon[2] = fromSeq.getCharAt(codonPos[2] - cdsStart);
  97       }
  98       else
  99       {
 100         baseCodon = null;
 101       }
 102     }
 103     else
 104     {
 105       codonPos = null;
 106       baseCodon = null; // todo tidy!
 107     }
 108   }
 109
 110   /**
 111    * Computes and returns comma-delimited HGVS notation peptide variants derived
 112    * from codon allele variants. If no variants are found, answers an empty
 113    * string.
 114    *
 115    * @param sf
 116    *             a sequence feature (which must be one of those held in this
 117    *             object)
 118    * @return
 119    */
 120   public String findProteinVariants(SequenceFeature sf)
 121   {
 122     if (!features.contains(sf) || baseCodon == null)
 123     {
 124       return "";
 125     }
 126
 127     /*
 128      * VCF data may already contain the protein consequence
 129      */
 130     String hgvsp = sf.getValueAsString(CSQ, HGV_SP);
 131     if (hgvsp != null)
 132     {
 133       int colonPos = hgvsp.lastIndexOf(':');
 134       if (colonPos >= 0)
 135       {
 136         String var = hgvsp.substring(colonPos + 1);
 137         return var;
 138       }
 139     }
 140
 141     /*
 142      * otherwise, compute codon and peptide variant
 143      */
 144     // todo avoid duplication of code in AlignmentUtils.buildDnaVariantsMap
 145     int cdsPos = sf.getBegin();
 146     if (cdsPos != sf.getEnd())
 147     {
 148       // not handling multi-locus variant features
 149       return "";
 150     }
 151     if (cdsPos != codonPos[0] && cdsPos != codonPos[1]
 152             && cdsPos != codonPos[2])
 153     {
 154       // e.g. feature on intron within spliced codon!
 155       return "";
 156     }
 157
 158     String alls = (String) sf.getValue(Gff3Helper.ALLELES);
 159     if (alls == null)
 160     {
 161       return "";
 162     }
 163
 164     String from3 = StringUtils.toSentenceCase(
 165             ResidueProperties.aa2Triplet.get(String.valueOf(toResidue)));
 166
 167     /*
 168      * make a peptide variant for each SNP allele
 169      * e.g. C,G,T gives variants G and T for base C
 170      */
 171     Set<String> variantPeptides = new HashSet<>();
 172     String[] alleles = alls.toUpperCase().split(",");
 173     StringBuilder vars = new StringBuilder();
 174
 175     for (String allele : alleles)
 176     {
 177       allele = allele.trim().toUpperCase();
 178       if (allele.length() > 1 || "-".equals(allele))
 179       {
 180         continue; // multi-locus variant
 181       }
 182       char[] variantCodon = new char[3];
 183       variantCodon[0] = baseCodon[0];
 184       variantCodon[1] = baseCodon[1];
 185       variantCodon[2] = baseCodon[2];
 186
 187       /*
 188        * poke variant base into canonical codon;
 189        * ignore first 'allele' (canonical base)
 190        */
 191       final int i = cdsPos == codonPos[0] ? 0
 192               : (cdsPos == codonPos[1] ? 1 : 2);
 193       variantCodon[i] = allele.toUpperCase().charAt(0);
 194       if (variantCodon[i] == baseCodon[i])
 195       {
 196         continue;
 197       }
 198       String codon = new String(variantCodon);
 199       String peptide = ResidueProperties.codonTranslate(codon);
 200       boolean synonymous = toResidue == peptide.charAt(0);
 201       StringBuilder var = new StringBuilder();
 202       if (synonymous)
 203       {
 204         /*
 205          * synonymous variant notation e.g. c.1062C>A(p.=)
 206          */
 207         var.append("c.").append(String.valueOf(cdsPos))
 208                 .append(String.valueOf(baseCodon[i])).append(">")
 209                 .append(String.valueOf(variantCodon[i]))
 210                 .append("(p.=)");
 211       }
 212       else
 213       {
 214         /*
 215          * missense variant notation e.g. p.Arg355Met
 216          */
 217         String to3 = ResidueProperties.STOP.equals(peptide) ? "Ter"
 218                 : StringUtils.toSentenceCase(
 219                         ResidueProperties.aa2Triplet.get(peptide));
 220         var.append("p.").append(from3).append(String.valueOf(toPosition))
 221                 .append(to3);
 222       }
 223       if (!variantPeptides.contains(peptide)) // duplicate consequence
 224       {
 225         variantPeptides.add(peptide);
 226         if (vars.length() > 0)
 227         {
 228           vars.append(",");
 229         }
 230         vars.append(var);
 231       }
 232     }
 233
 234     return vars.toString();
 235   }
 236 }