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     int[] codonPositions = codonIntervals == null ? null
  88             : MappingUtils.flattenRanges(codonIntervals);
  89     if (codonPositions != null && codonPositions.length == 3)
  90     {
  91       codonPos = codonPositions;
  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       codonPos = null;
 101       baseCodon = null;
 102     }
 103   }
 104
 105   /**
 106    * Computes and returns comma-delimited HGVS notation peptide variants derived
 107    * from codon allele variants. If no variants are found, answers an empty
 108    * string.
 109    *
 110    * @param sf
 111    *             a sequence feature (which must be one of those held in this
 112    *             object)
 113    * @return
 114    */
 115   public String findProteinVariants(SequenceFeature sf)
 116   {
 117     if (!features.contains(sf) || baseCodon == null)
 118     {
 119       return "";
 120     }
 121
 122     /*
 123      * VCF data may already contain the protein consequence
 124      */
 125     String hgvsp = sf.getValueAsString(CSQ, HGV_SP);
 126     if (hgvsp != null)
 127     {
 128       int colonPos = hgvsp.lastIndexOf(':');
 129       if (colonPos >= 0)
 130       {
 131         String var = hgvsp.substring(colonPos + 1);
 132         if (var.contains("p.")) // sanity check
 133         {
 134           return var;
 135         }
 136       }
 137     }
 138
 139     /*
 140      * otherwise, compute codon and peptide variant
 141      */
 142     int cdsPos = sf.getBegin();
 143     if (cdsPos != sf.getEnd())
 144     {
 145       // not handling multi-locus variant features
 146       return "";
 147     }
 148     if (cdsPos != codonPos[0] && cdsPos != codonPos[1]
 149             && cdsPos != codonPos[2])
 150     {
 151       // e.g. feature on intron within spliced codon!
 152       return "";
 153     }
 154
 155     String alls = (String) sf.getValue(Gff3Helper.ALLELES);
 156     if (alls == null)
 157     {
 158       return "";
 159     }
 160
 161     String from3 = StringUtils.toSentenceCase(
 162             ResidueProperties.aa2Triplet.get(String.valueOf(toResidue)));
 163
 164     /*
 165      * make a peptide variant for each SNP allele
 166      * e.g. C,G,T gives variants G and T for base C
 167      */
 168     Set<String> variantPeptides = new HashSet<>();
 169     String[] alleles = alls.toUpperCase().split(",");
 170     StringBuilder vars = new StringBuilder();
 171
 172     for (String allele : alleles)
 173     {
 174       allele = allele.trim().toUpperCase();
 175       if (allele.length() > 1 || "-".equals(allele))
 176       {
 177         continue; // multi-locus variant
 178       }
 179       char[] variantCodon = new char[3];
 180       variantCodon[0] = baseCodon[0];
 181       variantCodon[1] = baseCodon[1];
 182       variantCodon[2] = baseCodon[2];
 183
 184       /*
 185        * poke variant base into canonical codon;
 186        * ignore first 'allele' (canonical base)
 187        */
 188       final int i = cdsPos == codonPos[0] ? 0
 189               : (cdsPos == codonPos[1] ? 1 : 2);
 190       variantCodon[i] = allele.toUpperCase().charAt(0);
 191       if (variantCodon[i] == baseCodon[i])
 192       {
 193         continue;
 194       }
 195       String codon = new String(variantCodon);
 196       String peptide = ResidueProperties.codonTranslate(codon);
 197       boolean synonymous = toResidue == peptide.charAt(0);
 198       StringBuilder var = new StringBuilder();
 199       if (synonymous)
 200       {
 201         /*
 202          * synonymous variant notation e.g. c.1062C>A(p.=)
 203          */
 204         var.append("c.").append(String.valueOf(cdsPos))
 205                 .append(String.valueOf(baseCodon[i])).append(">")
 206                 .append(String.valueOf(variantCodon[i]))
 207                 .append("(p.=)");
 208       }
 209       else
 210       {
 211         /*
 212          * missense variant notation e.g. p.Arg355Met
 213          */
 214         String to3 = ResidueProperties.STOP.equals(peptide) ? "Ter"
 215                 : StringUtils.toSentenceCase(
 216                         ResidueProperties.aa2Triplet.get(peptide));
 217         var.append("p.").append(from3).append(String.valueOf(toPosition))
 218                 .append(to3);
 219       }
 220       if (!variantPeptides.contains(peptide)) // duplicate consequence
 221       {
 222         variantPeptides.add(peptide);
 223         if (vars.length() > 0)
 224         {
 225           vars.append(",");
 226         }
 227         vars.append(var);
 228       }
 229     }
 230
 231     return vars.toString();
 232   }
 233 }