+ /*
+ * skip if mapping is not to one of the target set of proteins
+ */
+ SequenceI proteinProduct = aMapping.getTo();
+ if (productSeqs != null && !productSeqs.contains(proteinProduct))
+ {
+ continue;
+ }
+
+ /*
+ * try to locate the CDS from the dataset mappings;
+ * guard against duplicate results (for the case that protein has
+ * dbrefs to both dna and cds sequences)
+ */
+ SequenceI cdsSeq = findCdsForProtein(mappings, dnaSeq,
+ seqMappings, aMapping);
+ if (cdsSeq != null)
+ {
+ if (!foundSeqs.contains(cdsSeq))
+ {
+ foundSeqs.add(cdsSeq);
+ SequenceI derivedSequence = cdsSeq.deriveSequence();
+ cdsSeqs.add(derivedSequence);
+ if (!dataset.getSequences().contains(cdsSeq))
+ {
+ dataset.addSequence(cdsSeq);
+ }
+ }
+ continue;
+ }
+
+ /*
+ * didn't find mapped CDS sequence - construct it and add
+ * its dataset sequence to the dataset
+ */
+ cdsSeq = makeCdsSequence(dnaSeq.getDatasetSequence(), aMapping,
+ dataset).deriveSequence();
+ // cdsSeq has a name constructed as CDS|<dbref>
+ // <dbref> will be either the accession for the coding sequence,
+ // marked in the /via/ dbref to the protein product accession
+ // or it will be the original nucleotide accession.
+ SequenceI cdsSeqDss = cdsSeq.getDatasetSequence();
+
+ cdsSeqs.add(cdsSeq);
+
+ if (!dataset.getSequences().contains(cdsSeqDss))
+ {
+ // check if this sequence is a newly created one
+ // so needs adding to the dataset
+ dataset.addSequence(cdsSeqDss);
+ }
+
+ /*
+ * add a mapping from CDS to the (unchanged) mapped to range
+ */
+ List<int[]> cdsRange = Collections
+ .singletonList(new int[]
+ { 1, cdsSeq.getLength() });
+ MapList cdsToProteinMap = new MapList(cdsRange,
+ mapList.getToRanges(), mapList.getFromRatio(),
+ mapList.getToRatio());
+ AlignedCodonFrame cdsToProteinMapping = new AlignedCodonFrame();
+ cdsToProteinMapping.addMap(cdsSeqDss, proteinProduct,
+ cdsToProteinMap);
+
+ /*
+ * guard against duplicating the mapping if repeating this action
+ */
+ if (!mappings.contains(cdsToProteinMapping))
+ {
+ mappings.add(cdsToProteinMapping);
+ }
+
+ propagateDBRefsToCDS(cdsSeqDss, dnaSeq.getDatasetSequence(),
+ proteinProduct, aMapping);
+ /*
+ * add another mapping from original 'from' range to CDS
+ */
+ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame();
+ MapList dnaToCdsMap = new MapList(mapList.getFromRanges(),
+ cdsRange, 1, 1);
+ dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss,
+ dnaToCdsMap);
+ if (!mappings.contains(dnaToCdsMapping))
+ {
+ mappings.add(dnaToCdsMapping);
+ }
+
+ /*
+ * add DBRef with mapping from protein to CDS
+ * (this enables Get Cross-References from protein alignment)
+ * This is tricky because we can't have two DBRefs with the
+ * same source and accession, so need a different accession for
+ * the CDS from the dna sequence
+ */
+
+ // specific use case:
+ // Genomic contig ENSCHR:1, contains coding regions for ENSG01,
+ // ENSG02, ENSG03, with transcripts and products similarly named.
+ // cannot add distinct dbrefs mapping location on ENSCHR:1 to ENSG01
+
+ // JBPNote: ?? can't actually create an example that demonstrates we
+ // need to
+ // synthesize an xref.
+
+ for (DBRefEntry primRef : dnaDss.getPrimaryDBRefs())
+ {
+ // creates a complementary cross-reference to the source sequence's
+ // primary reference.
+
+ DBRefEntry cdsCrossRef = new DBRefEntry(primRef.getSource(),
+ primRef.getSource() + ":" + primRef.getVersion(),
+ primRef.getAccessionId());
+ cdsCrossRef
+ .setMap(new Mapping(dnaDss, new MapList(dnaToCdsMap)));
+ cdsSeqDss.addDBRef(cdsCrossRef);
+
+ // problem here is that the cross-reference is synthesized -
+ // cdsSeq.getName() may be like 'CDS|dnaaccession' or
+ // 'CDS|emblcdsacc'
+ // assuming cds version same as dna ?!?
+
+ DBRefEntry proteinToCdsRef = new DBRefEntry(primRef.getSource(),
+ primRef.getVersion(), cdsSeq.getName());
+ //
+ proteinToCdsRef.setMap(
+ new Mapping(cdsSeqDss, cdsToProteinMap.getInverse()));
+ proteinProduct.addDBRef(proteinToCdsRef);
+ }
+
+ /*
+ * transfer any features on dna that overlap the CDS
+ */
+ transferFeatures(dnaSeq, cdsSeq, dnaToCdsMap, null,
+ SequenceOntologyI.CDS);
+ }
+ }
+ }
+
+ AlignmentI cds = new Alignment(
+ cdsSeqs.toArray(new SequenceI[cdsSeqs.size()]));
+ cds.setDataset(dataset);
+
+ return cds;
+ }
+
+ /**
+ * A helper method that finds a CDS sequence in the alignment dataset that is
+ * mapped to the given protein sequence, and either is, or has a mapping from,
+ * the given dna sequence.
+ *
+ * @param mappings
+ * set of all mappings on the dataset
+ * @param dnaSeq
+ * a dna (or cds) sequence we are searching from
+ * @param seqMappings
+ * the set of mappings involving dnaSeq
+ * @param aMapping
+ * an initial candidate from seqMappings
+ * @return
+ */
+ static SequenceI findCdsForProtein(List<AlignedCodonFrame> mappings,
+ SequenceI dnaSeq, List<AlignedCodonFrame> seqMappings,
+ Mapping aMapping)
+ {
+ /*
+ * TODO a better dna-cds-protein mapping data representation to allow easy
+ * navigation; until then this clunky looping around lists of mappings
+ */
+ SequenceI seqDss = dnaSeq.getDatasetSequence() == null ? dnaSeq
+ : dnaSeq.getDatasetSequence();
+ SequenceI proteinProduct = aMapping.getTo();
+
+ /*
+ * is this mapping from the whole dna sequence (i.e. CDS)?
+ * allowing for possible stop codon on dna but not peptide
+ */
+ int mappedFromLength = MappingUtils
+ .getLength(aMapping.getMap().getFromRanges());
+ int dnaLength = seqDss.getLength();
+ if (mappedFromLength == dnaLength
+ || mappedFromLength == dnaLength - CODON_LENGTH)
+ {
+ return seqDss;
+ }
+
+ /*
+ * looks like we found the dna-to-protein mapping; search for the
+ * corresponding cds-to-protein mapping
+ */
+ List<AlignedCodonFrame> mappingsToPeptide = MappingUtils
+ .findMappingsForSequence(proteinProduct, mappings);
+ for (AlignedCodonFrame acf : mappingsToPeptide)
+ {
+ for (SequenceToSequenceMapping map : acf.getMappings())
+ {
+ Mapping mapping = map.getMapping();
+ if (mapping != aMapping
+ && mapping.getMap().getFromRatio() == CODON_LENGTH
+ && proteinProduct == mapping.getTo()
+ && seqDss != map.getFromSeq())
+ {
+ mappedFromLength = MappingUtils
+ .getLength(mapping.getMap().getFromRanges());
+ if (mappedFromLength == map.getFromSeq().getLength())
+ {
+ /*
+ * found a 3:1 mapping to the protein product which covers
+ * the whole dna sequence i.e. is from CDS; finally check it
+ * is from the dna start sequence
+ */
+ SequenceI cdsSeq = map.getFromSeq();
+ List<AlignedCodonFrame> dnaToCdsMaps = MappingUtils
+ .findMappingsForSequence(cdsSeq, seqMappings);
+ if (!dnaToCdsMaps.isEmpty())
+ {
+ return cdsSeq;
+ }
+ }
+ }
+ }
+ }
+ return null;
+ }
+
+ /**
+ * Helper method that makes a CDS sequence as defined by the mappings from the
+ * given sequence i.e. extracts the 'mapped from' ranges (which may be on
+ * forward or reverse strand).
+ *
+ * @param seq
+ * @param mapping
+ * @param dataset
+ * - existing dataset. We check for sequences that look like the CDS
+ * we are about to construct, if one exists already, then we will
+ * just return that one.
+ * @return CDS sequence (as a dataset sequence)
+ */
+ static SequenceI makeCdsSequence(SequenceI seq, Mapping mapping,
+ AlignmentI dataset)
+ {
+ char[] seqChars = seq.getSequence();
+ List<int[]> fromRanges = mapping.getMap().getFromRanges();
+ int cdsWidth = MappingUtils.getLength(fromRanges);
+ char[] newSeqChars = new char[cdsWidth];
+
+ int newPos = 0;
+ for (int[] range : fromRanges)
+ {
+ if (range[0] <= range[1])
+ {
+ // forward strand mapping - just copy the range
+ int length = range[1] - range[0] + 1;
+ System.arraycopy(seqChars, range[0] - 1, newSeqChars, newPos,
+ length);
+ newPos += length;
+ }
+ else
+ {
+ // reverse strand mapping - copy and complement one by one
+ for (int i = range[0]; i >= range[1]; i--)
+ {
+ newSeqChars[newPos++] = Dna.getComplement(seqChars[i - 1]);
+ }
+ }
+ }
+
+ /*
+ * assign 'from id' held in the mapping if set (e.g. EMBL protein_id),
+ * else generate a sequence name
+ */
+ String mapFromId = mapping.getMappedFromId();
+ String seqId = "CDS|" + (mapFromId != null ? mapFromId : seq.getName());
+ SequenceI newSeq = new Sequence(seqId, newSeqChars, 1, newPos);
+ if (dataset != null)
+ {
+ SequenceI[] matches = dataset.findSequenceMatch(newSeq.getName());
+ if (matches != null)
+ {
+ boolean matched = false;
+ for (SequenceI mtch : matches)
+ {
+ if (mtch.getStart() != newSeq.getStart())
+ {
+ continue;
+ }
+ if (mtch.getEnd() != newSeq.getEnd())
+ {
+ continue;
+ }
+ if (!Arrays.equals(mtch.getSequence(), newSeq.getSequence()))
+ {
+ continue;
+ }
+ if (!matched)
+ {
+ matched = true;
+ newSeq = mtch;
+ }
+ else
+ {
+ System.err.println(
+ "JAL-2154 regression: warning - found (and ignnored a duplicate CDS sequence):"
+ + mtch.toString());
+ }
+ }
+ }
+ }
+ // newSeq.setDescription(mapFromId);
+
+ return newSeq;
+ }
+
+ /**
+ * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to
+ * the given mapping.
+ *
+ * @param cdsSeq
+ * @param contig
+ * @param mapping
+ * @return list of DBRefEntrys added.
+ */
+ public static List<DBRefEntry> propagateDBRefsToCDS(SequenceI cdsSeq,
+ SequenceI contig, SequenceI proteinProduct, Mapping mapping)
+ {
+
+ // gather direct refs from contig congrent with mapping
+ List<DBRefEntry> direct = new ArrayList<DBRefEntry>();
+ HashSet<String> directSources = new HashSet<String>();
+ if (contig.getDBRefs() != null)
+ {
+ for (DBRefEntry dbr : contig.getDBRefs())
+ {
+ if (dbr.hasMap() && dbr.getMap().getMap().isTripletMap())
+ {
+ MapList map = dbr.getMap().getMap();
+ // check if map is the CDS mapping
+ if (mapping.getMap().equals(map))
+ {
+ direct.add(dbr);
+ directSources.add(dbr.getSource());
+ }
+ }
+ }
+ }
+ DBRefEntry[] onSource = DBRefUtils.selectRefs(
+ proteinProduct.getDBRefs(),
+ directSources.toArray(new String[0]));
+ List<DBRefEntry> propagated = new ArrayList<DBRefEntry>();
+
+ // and generate appropriate mappings
+ for (DBRefEntry cdsref : direct)
+ {
+ // clone maplist and mapping
+ MapList cdsposmap = new MapList(
+ Arrays.asList(new int[][]
+ { new int[] { cdsSeq.getStart(), cdsSeq.getEnd() } }),
+ cdsref.getMap().getMap().getToRanges(), 3, 1);
+ Mapping cdsmap = new Mapping(cdsref.getMap().getTo(),
+ cdsref.getMap().getMap());
+
+ // create dbref
+ DBRefEntry newref = new DBRefEntry(cdsref.getSource(),
+ cdsref.getVersion(), cdsref.getAccessionId(),
+ new Mapping(cdsmap.getTo(), cdsposmap));
+
+ // and see if we can map to the protein product for this mapping.
+ // onSource is the filtered set of accessions on protein that we are
+ // tranferring, so we assume accession is the same.
+ if (cdsmap.getTo() == null && onSource != null)
+ {
+ List<DBRefEntry> sourceRefs = DBRefUtils.searchRefs(onSource,
+ cdsref.getAccessionId());
+ if (sourceRefs != null)
+ {
+ for (DBRefEntry srcref : sourceRefs)
+ {
+ if (srcref.getSource().equalsIgnoreCase(cdsref.getSource()))
+ {
+ // we have found a complementary dbref on the protein product, so
+ // update mapping's getTo
+ newref.getMap().setTo(proteinProduct);
+ }
+ }
+ }
+ }
+ cdsSeq.addDBRef(newref);
+ propagated.add(newref);
+ }
+ return propagated;
+ }
+
+ /**
+ * Transfers co-located features on 'fromSeq' to 'toSeq', adjusting the
+ * feature start/end ranges, optionally omitting specified feature types.
+ * Returns the number of features copied.
+ *
+ * @param fromSeq
+ * @param toSeq
+ * @param mapping
+ * the mapping from 'fromSeq' to 'toSeq'
+ * @param select
+ * if not null, only features of this type are copied (including
+ * subtypes in the Sequence Ontology)
+ * @param omitting
+ */
+ public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq,
+ MapList mapping, String select, String... omitting)
+ {
+ SequenceI copyTo = toSeq;
+ while (copyTo.getDatasetSequence() != null)
+ {
+ copyTo = copyTo.getDatasetSequence();
+ }
+
+ /*
+ * get features, optionally restricted by an ontology term
+ */
+ List<SequenceFeature> sfs = select == null ? fromSeq.getFeatures()
+ .getPositionalFeatures() : fromSeq.getFeatures()
+ .getFeaturesByOntology(select);
+
+ int count = 0;
+ for (SequenceFeature sf : sfs)
+ {
+ String type = sf.getType();
+ boolean omit = false;
+ for (String toOmit : omitting)
+ {
+ if (type.equals(toOmit))
+ {
+ omit = true;
+ }
+ }
+ if (omit)
+ {
+ continue;
+ }
+
+ /*
+ * locate the mapped range - null if either start or end is
+ * not mapped (no partial overlaps are calculated)
+ */
+ int start = sf.getBegin();
+ int end = sf.getEnd();
+ int[] mappedTo = mapping.locateInTo(start, end);
+ /*
+ * if whole exon range doesn't map, try interpreting it
+ * as 5' or 3' exon overlapping the CDS range
+ */
+ if (mappedTo == null)
+ {
+ mappedTo = mapping.locateInTo(end, end);