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
  15  *
  16  * @author gmcarstairs
  17  */
  18 public class MappedFeatures
  19 {
  20   private static final String HGV_SP = "HGVSp";
  21
  22   private static final String CSQ = "CSQ";
  23
  24   /*
  25    * the mapping from one sequence to another
  26    */
  27   public final Mapping mapping;
  28
  29   /**
  30    * the sequence mapped to
  31    */
  32   public final SequenceI fromSeq;
  33
  34   /*
  35    * the residue position in the sequence mapped from
  36    */
  37   public final int fromPosition;
  38
  39   /*
  40    * the residue at fromPosition
  41    */
  42   public final char fromResidue;
  43
  44   /*
  45    * features on the sequence mapped to that overlap the mapped positions
  46    */
  47   public final List<SequenceFeature> features;
  48
  49   /*
  50    * if the mapping is 1:3 (peptide to CDS), this holds the
  51    * mapped positions i.e. codon base positions in CDS; to
  52    * support calculation of peptide variants from alleles
  53    */
  54   public final int[] codonPos;
  55
  56   private final char[] baseCodon;
  57
  58   /**
  59    * Constructor
  60    *
  61    * @param theMapping
  62    * @param pos
  63    * @param res
  64    * @param theFeatures
  65    */
  66   public MappedFeatures(Mapping theMapping, SequenceI from, int pos,
  67           char res,
  68           List<SequenceFeature> theFeatures)
  69   {
  70     mapping = theMapping;
  71     fromSeq = from;
  72     fromPosition = pos;
  73     fromResidue = res;
  74     features = theFeatures;
  75
  76     /*
  77      * determine codon positions and canonical codon
  78      * for a peptide-to-CDS mapping
  79      */
  80     codonPos = MappingUtils.flattenRanges(
  81             mapping.getMap().locateInFrom(fromPosition, fromPosition));
  82     if (codonPos.length == 3)
  83     {
  84       baseCodon = new char[3];
  85       int cdsStart = fromSeq.getStart();
  86       baseCodon[0] = fromSeq.getCharAt(codonPos[0] - cdsStart);
  87       baseCodon[1] = fromSeq.getCharAt(codonPos[1] - cdsStart);
  88       baseCodon[2] = fromSeq.getCharAt(codonPos[2] - cdsStart);
  89     }
  90     else
  91     {
  92       baseCodon = null;
  93     }
  94   }
  95
  96   /**
  97    * Computes and returns comma-delimited HGVS notation peptide variants derived
  98    * from codon allele variants. If no variants are found, answers an empty
  99    * string.
 100    *
 101    * @return
 102    */
 103   public String findProteinVariants(SequenceFeature sf)
 104   {
 105     if (!features.contains(sf) || baseCodon == null)
 106     {
 107       return "";
 108     }
 109
 110     StringBuilder vars = new StringBuilder();
 111
 112     /*
 113      * VCF data may already contain the protein consequence
 114      */
 115     String hgvsp = sf.getValueAsString(CSQ, HGV_SP);
 116     if (hgvsp != null)
 117     {
 118       int colonPos = hgvsp.indexOf(':');
 119       if (colonPos >= 0)
 120       {
 121         String var = hgvsp.substring(colonPos + 1);
 122         return var;
 123       }
 124     }
 125
 126     /*
 127      * otherwise, compute codon and peptide variant
 128      */
 129     // todo avoid duplication of code in AlignmentUtils.buildDnaVariantsMap
 130     int cdsPos = sf.getBegin();
 131     if (cdsPos != sf.getEnd())
 132     {
 133       // not handling multi-locus variant features
 134       return "";
 135     }
 136     if (cdsPos != codonPos[0] && cdsPos != codonPos[1]
 137             && cdsPos != codonPos[2])
 138     {
 139       // e.g. feature on intron within spliced codon!
 140       return "";
 141     }
 142
 143     String alls = (String) sf.getValue(Gff3Helper.ALLELES);
 144     if (alls == null)
 145     {
 146       return "";
 147     }
 148
 149     String from3 = StringUtils.toSentenceCase(
 150             ResidueProperties.aa2Triplet.get(String.valueOf(fromResidue)));
 151
 152     /*
 153      * make a peptide variant for each SNP allele
 154      * e.g. C,G,T gives variants G and T for base C
 155      */
 156     Set<String> variantPeptides = new HashSet<>();
 157     String[] alleles = alls.toUpperCase().split(",");
 158     for (String allele : alleles)
 159     {
 160       allele = allele.trim().toUpperCase();
 161       if (allele.length() > 1 || "-".equals(allele))
 162       {
 163         continue; // multi-locus variant
 164       }
 165       char[] variantCodon = new char[3];
 166       variantCodon[0] = baseCodon[0];
 167       variantCodon[1] = baseCodon[1];
 168       variantCodon[2] = baseCodon[2];
 169
 170       /*
 171        * poke variant base into canonical codon
 172        */
 173       int i = cdsPos == codonPos[0] ? 0 : (cdsPos == codonPos[1] ? 1 : 2);
 174       variantCodon[i] = allele.toUpperCase().charAt(0);
 175       String codon = new String(variantCodon);
 176       String peptide = ResidueProperties.codonTranslate(codon);
 177       if (fromResidue != peptide.charAt(0))
 178       {
 179         String to3 = ResidueProperties.STOP.equals(peptide) ? "STOP"
 180                 : StringUtils.toSentenceCase(
 181                         ResidueProperties.aa2Triplet.get(peptide));
 182         String var = "p." + from3 + fromPosition + to3;
 183         if (!variantPeptides.contains(peptide)) // duplicate consequence
 184         {
 185           variantPeptides.add(peptide);
 186           if (vars.length() > 0)
 187           {
 188             vars.append(",");
 189           }
 190           vars.append(var);
 191         }
 192       }
 193     }
 194
 195     return vars.toString();
 196   }
 197 }