+ /*
+ * variants in third codon base
+ */
+ for (DnaVariant var : codonVariants[2])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES);
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ if (!base3.equalsIgnoreCase(base))
+ {
+ String codon = base1.toLowerCase() + base2.toLowerCase()
+ + base.toUpperCase();
+ String canonical = base1.toLowerCase() + base2.toLowerCase()
+ + base3.toUpperCase();
+ if (addPeptideVariant(peptide, peptidePos, residue, var,
+ codon, canonical))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return count;
+ }
+
+ /**
+ * Helper method that adds a peptide variant feature. ID and
+ * clinical_significance attributes of the dna variant (if present) are copied
+ * to the new feature.
+ *
+ * @param peptide
+ * @param peptidePos
+ * @param residue
+ * @param var
+ * @param codon
+ * the variant codon e.g. aCg
+ * @param canonical
+ * the 'normal' codon e.g. aTg
+ * @return true if a feature was added, else false
+ */
+ static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
+ String residue, DnaVariant var, String codon, String canonical)
+ {
+ /*
+ * get peptide translation of codon e.g. GAT -> D
+ * note that variants which are not single alleles,
+ * e.g. multibase variants or HGMD_MUTATION etc
+ * are currently ignored here
+ */
+ String trans = codon.contains("-") ? null
+ : (codon.length() > CODON_LENGTH ? null
+ : ResidueProperties.codonTranslate(codon));
+ if (trans == null)
+ {
+ return false;
+ }
+ String desc = canonical + "/" + codon;
+ String featureType = "";
+ if (trans.equals(residue))
+ {
+ featureType = SequenceOntologyI.SYNONYMOUS_VARIANT;
+ }
+ else if (ResidueProperties.STOP.equals(trans))
+ {
+ featureType = SequenceOntologyI.STOP_GAINED;
+ }
+ else
+ {
+ String residue3Char = StringUtils
+ .toSentenceCase(ResidueProperties.aa2Triplet.get(residue));
+ String trans3Char = StringUtils
+ .toSentenceCase(ResidueProperties.aa2Triplet.get(trans));
+ desc = "p." + residue3Char + peptidePos + trans3Char;
+ featureType = SequenceOntologyI.NONSYNONYMOUS_VARIANT;
+ }
+ SequenceFeature sf = new SequenceFeature(featureType, desc, peptidePos,
+ peptidePos, var.getSource());
+
+ StringBuilder attributes = new StringBuilder(32);
+ String id = (String) var.variant.getValue(ID);
+ if (id != null)
+ {
+ if (id.startsWith(SEQUENCE_VARIANT))
+ {
+ id = id.substring(SEQUENCE_VARIANT.length());
+ }
+ sf.setValue(ID, id);
+ attributes.append(ID).append("=").append(id);
+ // TODO handle other species variants JAL-2064
+ StringBuilder link = new StringBuilder(32);
+ try
+ {
+ link.append(desc).append(" ").append(id).append(
+ "|http://www.ensembl.org/Homo_sapiens/Variation/Summary?v=")
+ .append(URLEncoder.encode(id, "UTF-8"));
+ sf.addLink(link.toString());
+ } catch (UnsupportedEncodingException e)
+ {
+ // as if
+ }
+ }
+ String clinSig = (String) var.variant.getValue(CLINICAL_SIGNIFICANCE);
+ if (clinSig != null)
+ {
+ sf.setValue(CLINICAL_SIGNIFICANCE, clinSig);
+ attributes.append(";").append(CLINICAL_SIGNIFICANCE).append("=")
+ .append(clinSig);
+ }
+ peptide.addSequenceFeature(sf);
+ if (attributes.length() > 0)
+ {
+ sf.setAttributes(attributes.toString());
+ }
+ return true;
+ }
+
+ /**
+ * Builds a map whose key is position in the protein sequence, and value is a
+ * list of the base and all variants for each corresponding codon position.
+ * <p>
+ * This depends on dna variants being held as a comma-separated list as
+ * property "alleles" on variant features.
+ *
+ * @param dnaSeq
+ * @param dnaToProtein
+ * @return
+ */
+ @SuppressWarnings("unchecked")
+ static LinkedHashMap<Integer, List<DnaVariant>[]> buildDnaVariantsMap(
+ SequenceI dnaSeq, MapList dnaToProtein)
+ {
+ /*
+ * map from peptide position to all variants of the codon which codes for it
+ * LinkedHashMap ensures we keep the peptide features in sequence order
+ */
+ LinkedHashMap<Integer, List<DnaVariant>[]> variants = new LinkedHashMap<>();
+
+ List<SequenceFeature> dnaFeatures = dnaSeq.getFeatures()
+ .getFeaturesByOntology(SequenceOntologyI.SEQUENCE_VARIANT);
+ if (dnaFeatures.isEmpty())
+ {
+ return variants;
+ }
+
+ int dnaStart = dnaSeq.getStart();
+ int[] lastCodon = null;
+ int lastPeptidePostion = 0;
+
+ /*
+ * build a map of codon variations for peptides
+ */
+ for (SequenceFeature sf : dnaFeatures)
+ {
+ int dnaCol = sf.getBegin();
+ if (dnaCol != sf.getEnd())
+ {
+ // not handling multi-locus variant features
+ continue;
+ }
+
+ /*
+ * ignore variant if not a SNP
+ */
+ String alls = (String) sf.getValue(Gff3Helper.ALLELES);
+ if (alls == null)
+ {
+ continue; // non-SNP VCF variant perhaps - can't process this
+ }
+
+ String[] alleles = alls.toUpperCase().split(",");
+ boolean isSnp = true;
+ for (String allele : alleles)
+ {
+ if (allele.trim().length() > 1)
+ {
+ isSnp = false;
+ }
+ }
+ if (!isSnp)
+ {
+ continue;
+ }
+
+ int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol);
+ if (mapsTo == null)
+ {
+ // feature doesn't lie within coding region
+ continue;
+ }
+ int peptidePosition = mapsTo[0];
+ List<DnaVariant>[] codonVariants = variants.get(peptidePosition);
+ if (codonVariants == null)
+ {
+ codonVariants = new ArrayList[CODON_LENGTH];
+ codonVariants[0] = new ArrayList<>();
+ codonVariants[1] = new ArrayList<>();
+ codonVariants[2] = new ArrayList<>();
+ variants.put(peptidePosition, codonVariants);
+ }
+
+ /*
+ * get this peptide's codon positions e.g. [3, 4, 5] or [4, 7, 10]
+ */
+ int[] codon = peptidePosition == lastPeptidePostion ? lastCodon
+ : MappingUtils.flattenRanges(dnaToProtein.locateInFrom(
+ peptidePosition, peptidePosition));
+ lastPeptidePostion = peptidePosition;
+ lastCodon = codon;
+
+ /*
+ * save nucleotide (and any variant) for each codon position
+ */
+ for (int codonPos = 0; codonPos < CODON_LENGTH; codonPos++)
+ {
+ String nucleotide = String.valueOf(
+ dnaSeq.getCharAt(codon[codonPos] - dnaStart)).toUpperCase();
+ List<DnaVariant> codonVariant = codonVariants[codonPos];
+ if (codon[codonPos] == dnaCol)