+ /*
+ * build a map with codon variations for each potentially varying peptide
+ */
+ LinkedHashMap<Integer, List<DnaVariant>[]> variants = buildDnaVariantsMap(
+ dnaSeq, dnaToProtein);
+
+ /*
+ * scan codon variations, compute peptide variants and add to peptide sequence
+ */
+ int count = 0;
+ for (Entry<Integer, List<DnaVariant>[]> variant : variants.entrySet())
+ {
+ int peptidePos = variant.getKey();
+ List<DnaVariant>[] codonVariants = variant.getValue();
+ count += computePeptideVariants(peptide, peptidePos, codonVariants);
+ }
+
+ /*
+ * sort to get sequence features in start position order
+ * - would be better to store in Sequence as a TreeSet or NCList?
+ */
+ Arrays.sort(peptide.getSequenceFeatures(),
+ new Comparator<SequenceFeature>()
+ {
+ @Override
+ public int compare(SequenceFeature o1, SequenceFeature o2)
+ {
+ int c = Integer.compare(o1.getBegin(), o2.getBegin());
+ return c == 0 ? Integer.compare(o1.getEnd(), o2.getEnd())
+ : c;
+ }
+ });
+ return count;
+ }
+
+ /**
+ * Computes non-synonymous peptide variants from codon variants and adds them
+ * as sequence_variant features on the protein sequence (one feature per
+ * allele variant). Selected attributes (variant id, clinical significance)
+ * are copied over to the new features.
+ *
+ * @param peptide
+ * the protein sequence
+ * @param peptidePos
+ * the position to compute peptide variants for
+ * @param codonVariants
+ * a list of dna variants per codon position
+ * @return the number of features added
+ */
+ static int computePeptideVariants(SequenceI peptide, int peptidePos,
+ List<DnaVariant>[] codonVariants)
+ {
+ String residue = String.valueOf(peptide.getCharAt(peptidePos - 1));
+ int count = 0;
+ String base1 = codonVariants[0].get(0).base;
+ String base2 = codonVariants[1].get(0).base;
+ String base3 = codonVariants[2].get(0).base;
+
+ /*
+ * variants in first codon base
+ */
+ for (DnaVariant var : codonVariants[0])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base + base2 + base3;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
+
+ /*
+ * variants in second codon base
+ */
+ for (DnaVariant var : codonVariants[1])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base1 + base + base3;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
+
+ /*
+ * variants in third codon base
+ */
+ for (DnaVariant var : codonVariants[2])
+ {
+ if (var.variant != null)
+ {
+ String alleles = (String) var.variant.getValue("alleles");
+ if (alleles != null)
+ {
+ for (String base : alleles.split(","))
+ {
+ String codon = base1 + base2 + base;
+ if (addPeptideVariant(peptide, peptidePos, residue, var, codon))
+ {
+ count++;
+ }
+ }
+ }
+ }
+ }
+
+ return count;
+ }
+
+ /**
+ * Helper method that adds a peptide variant feature, provided the given codon
+ * translates to a value different to the current residue (is a non-synonymous
+ * variant). 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
+ * @return true if a feature was added, else false
+ */
+ static boolean addPeptideVariant(SequenceI peptide, int peptidePos,
+ String residue, DnaVariant var, String codon)
+ {
+ /*
+ * 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("-") ? "-"
+ : (codon.length() > 3 ? null : ResidueProperties
+ .codonTranslate(codon));
+ if (trans != null && !trans.equals(residue))
+ {
+ String desc = residue + "->" + trans;
+ // set score to 0f so 'graduated colour' option is offered!
+ SequenceFeature sf = new SequenceFeature(
+ SequenceOntologyI.SEQUENCE_VARIANT, desc, peptidePos,
+ peptidePos, 0f, null);
+ 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);
+ // TODO handle other species variants
+ 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);
+ }
+ peptide.addSequenceFeature(sf);
+ return true;
+ }
+ return false;
+ }
+
+ /**
+ * 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
+ *
+ * @param dnaSeq
+ * @param dnaToProtein
+ * @return
+ */
+ 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<Integer, List<DnaVariant>[]>();
+ SequenceOntologyI so = SequenceOntologyFactory.getInstance();
+
+ SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures();
+ if (dnaFeatures == null)
+ {
+ 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;
+ }
+ if (so.isA(sf.getType(), SequenceOntologyI.SEQUENCE_VARIANT))
+ {
+ 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[3];
+ codonVariants[0] = new ArrayList<DnaVariant>();
+ codonVariants[1] = new ArrayList<DnaVariant>();
+ codonVariants[2] = new ArrayList<DnaVariant>();
+ variants.put(peptidePosition, codonVariants);
+ }
+
+ /*
+ * extract dna variants to a string array
+ */
+ String alls = (String) sf.getValue("alleles");
+ if (alls == null)
+ {
+ continue;
+ }
+ String[] alleles = alls.toUpperCase().split(",");
+ int i = 0;
+ for (String allele : alleles)
+ {
+ alleles[i++] = allele.trim(); // lose any space characters "A, G"
+ }
+
+ /*
+ * 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 < 3; codonPos++)
+ {
+ String nucleotide = String.valueOf(
+ dnaSeq.getCharAt(codon[codonPos] - dnaStart))
+ .toUpperCase();
+ List<DnaVariant> codonVariant = codonVariants[codonPos];
+ if (codon[codonPos] == dnaCol)
+ {
+ if (!codonVariant.isEmpty()
+ && codonVariant.get(0).variant == null)
+ {
+ /*
+ * already recorded base value, add this variant
+ */
+ codonVariant.get(0).variant = sf;
+ }
+ else
+ {
+ /*
+ * add variant with base value
+ */
+ codonVariant.add(new DnaVariant(nucleotide, sf));
+ }
+ }
+ else if (codonVariant.isEmpty())
+ {
+ /*
+ * record (possibly non-varying) base value
+ */
+ codonVariant.add(new DnaVariant(nucleotide));
+ }
+ }
+ }
+ }
+ return variants;
+ }
+
+ /**
+ * Makes an alignment with a copy of the given sequences, adding in any
+ * non-redundant sequences which are mapped to by the cross-referenced
+ * sequences.
+ *
+ * @param seqs
+ * @param xrefs
+ * @return
+ */
+ public static AlignmentI makeCopyAlignment(SequenceI[] seqs,
+ SequenceI[] xrefs)
+ {
+ AlignmentI copy = new Alignment(new Alignment(seqs));
+
+ SequenceIdMatcher matcher = new SequenceIdMatcher(seqs);
+ if (xrefs != null)
+ {
+ for (SequenceI xref : xrefs)
+ {
+ DBRefEntry[] dbrefs = xref.getDBRefs();
+ if (dbrefs != null)
+ {
+ for (DBRefEntry dbref : dbrefs)
+ {
+ if (dbref.getMap() == null || dbref.getMap().getTo() == null)
+ {
+ continue;
+ }
+ SequenceI mappedTo = dbref.getMap().getTo();
+ SequenceI match = matcher.findIdMatch(mappedTo);
+ if (match == null)
+ {
+ matcher.add(mappedTo);
+ copy.addSequence(mappedTo);
+ }
+ }
+ }
+ }
+ }
+ return copy;
+ }
+
+ /**
+ * Try to align sequences in 'unaligned' to match the alignment of their
+ * mapped regions in 'aligned'. For example, could use this to align CDS
+ * sequences which are mapped to their parent cDNA sequences.
+ *
+ * This method handles 1:1 mappings (dna-to-dna or protein-to-protein). For
+ * dna-to-protein or protein-to-dna use alternative methods.
+ *
+ * @param unaligned
+ * sequences to be aligned
+ * @param aligned
+ * holds aligned sequences and their mappings
+ * @return
+ */
+ public static int alignAs(AlignmentI unaligned, AlignmentI aligned)
+ {
+ List<SequenceI> unmapped = new ArrayList<SequenceI>();
+ Map<Integer, Map<SequenceI, Character>> columnMap = buildMappedColumnsMap(
+ unaligned, aligned, unmapped);
+ int width = columnMap.size();
+ char gap = unaligned.getGapCharacter();
+ int realignedCount = 0;
+
+ for (SequenceI seq : unaligned.getSequences())
+ {
+ if (!unmapped.contains(seq))
+ {
+ char[] newSeq = new char[width];
+ Arrays.fill(newSeq, gap);
+ int newCol = 0;
+ int lastCol = 0;
+
+ /*
+ * traverse the map to find columns populated
+ * by our sequence
+ */
+ for (Integer column : columnMap.keySet())
+ {
+ Character c = columnMap.get(column).get(seq);
+ if (c != null)
+ {
+ /*
+ * sequence has a character at this position
+ *
+ */
+ newSeq[newCol] = c;
+ lastCol = newCol;
+ }
+ newCol++;
+ }
+
+ /*
+ * trim trailing gaps
+ */
+ if (lastCol < width)
+ {
+ char[] tmp = new char[lastCol + 1];
+ System.arraycopy(newSeq, 0, tmp, 0, lastCol + 1);
+ newSeq = tmp;
+ }
+ seq.setSequence(String.valueOf(newSeq));
+ realignedCount++;
+ }
+ }
+ return realignedCount;