import java.util.HashSet;
import java.util.Iterator;
import java.util.LinkedHashMap;
+import java.util.LinkedHashSet;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
import jalview.datamodel.AlignedCodon;
import jalview.datamodel.AlignedCodonFrame;
+import jalview.datamodel.Alignment;
import jalview.datamodel.AlignmentAnnotation;
import jalview.datamodel.AlignmentI;
import jalview.datamodel.DBRefEntry;
+import jalview.datamodel.DBRefSource;
+import jalview.datamodel.FeatureProperties;
import jalview.datamodel.Mapping;
import jalview.datamodel.SearchResults;
import jalview.datamodel.Sequence;
import jalview.datamodel.SequenceGroup;
import jalview.datamodel.SequenceI;
import jalview.schemes.ResidueProperties;
+import jalview.util.DBRefUtils;
import jalview.util.MapList;
+import jalview.util.MappingUtils;
/**
* grab bag of useful alignment manipulation operations Expect these to be
}
return false;
}
+
+ /**
+ * Constructs an alignment consisting of the mapped exon regions in the given
+ * nucleotide sequences, and updates mappings to match.
+ *
+ * @param dna
+ * aligned dna sequences
+ * @param mappings
+ * from dna to protein; these are replaced with new mappings
+ * @return an alignment whose sequences are the exon-only parts of the dna
+ * sequences (or null if no exons are found)
+ */
+ public static AlignmentI makeExonAlignment(SequenceI[] dna,
+ Set<AlignedCodonFrame> mappings)
+ {
+ Set<AlignedCodonFrame> newMappings = new LinkedHashSet<AlignedCodonFrame>();
+ List<SequenceI> exonSequences = new ArrayList<SequenceI>();
+
+ for (SequenceI dnaSeq : dna)
+ {
+ final SequenceI ds = dnaSeq.getDatasetSequence();
+ List<AlignedCodonFrame> seqMappings = MappingUtils
+ .findMappingsForSequence(ds, mappings);
+ for (AlignedCodonFrame acf : seqMappings)
+ {
+ AlignedCodonFrame newMapping = new AlignedCodonFrame();
+ final List<SequenceI> mappedExons = makeExonSequences(ds, acf,
+ newMapping);
+ if (!mappedExons.isEmpty())
+ {
+ exonSequences.addAll(mappedExons);
+ newMappings.add(newMapping);
+ }
+ }
+ }
+ AlignmentI al = new Alignment(
+ exonSequences.toArray(new SequenceI[exonSequences.size()]));
+ al.setDataset(null);
+
+ /*
+ * Replace the old mappings with the new ones
+ */
+ mappings.clear();
+ mappings.addAll(newMappings);
+
+ return al;
+ }
+
+ /**
+ * Helper method to make exon-only sequences and populate their mappings to
+ * protein products
+ * <p>
+ * For example, if ggCCaTTcGAg has mappings [3, 4, 6, 7, 9, 10] to protein
+ * then generate a sequence CCTTGA with mapping [1, 6] to the same protein
+ * residues
+ * <p>
+ * Typically eukaryotic dna will include exons encoding for a single peptide
+ * sequence i.e. return a single result. Bacterial dna may have overlapping
+ * exon mappings coding for multiple peptides so return multiple results
+ * (example EMBL KF591215).
+ *
+ * @param dnaSeq
+ * a dna dataset sequence
+ * @param mapping
+ * containing one or more mappings of the sequence to protein
+ * @param newMapping
+ * the new mapping to populate, from the exon-only sequences to their
+ * mapped protein sequences
+ * @return
+ */
+ protected static List<SequenceI> makeExonSequences(SequenceI dnaSeq,
+ AlignedCodonFrame mapping, AlignedCodonFrame newMapping)
+ {
+ List<SequenceI> exonSequences = new ArrayList<SequenceI>();
+ List<Mapping> seqMappings = mapping.getMappingsForSequence(dnaSeq);
+ final char[] dna = dnaSeq.getSequence();
+ for (Mapping seqMapping : seqMappings)
+ {
+ StringBuilder newSequence = new StringBuilder(dnaSeq.getLength());
+
+ /*
+ * Get the codon regions as { [2, 5], [7, 12], [14, 14] etc }
+ */
+ final List<int[]> dnaExonRanges = seqMapping.getMap().getFromRanges();
+ for (int[] range : dnaExonRanges)
+ {
+ for (int pos = range[0]; pos <= range[1]; pos++)
+ {
+ newSequence.append(dna[pos - 1]);
+ }
+ }
+
+ SequenceI exon = new Sequence(dnaSeq.getName(),
+ newSequence.toString());
+
+ /*
+ * Locate any xrefs to CDS database on the protein product and attach to
+ * the CDS sequence. Also add as a sub-token of the sequence name.
+ */
+ // default to "CDS" if we can't locate an actual gene id
+ String cdsAccId = FeatureProperties
+ .getCodingFeature(DBRefSource.EMBL);
+ DBRefEntry[] cdsRefs = DBRefUtils.selectRefs(seqMapping.getTo()
+ .getDBRef(), DBRefSource.CODINGDBS);
+ if (cdsRefs != null)
+ {
+ for (DBRefEntry cdsRef : cdsRefs)
+ {
+ exon.addDBRef(new DBRefEntry(cdsRef));
+ cdsAccId = cdsRef.getAccessionId();
+ }
+ }
+ exon.setName(exon.getName() + "|" + cdsAccId);
+ exon.createDatasetSequence();
+
+ /*
+ * Build new mappings - from the same protein regions, but now to
+ * contiguous exons
+ */
+ List<int[]> exonRange = new ArrayList<int[]>();
+ exonRange.add(new int[]
+ { 1, newSequence.length() });
+ MapList map = new MapList(exonRange, seqMapping.getMap()
+ .getToRanges(),
+ 3, 1);
+ newMapping.addMap(exon.getDatasetSequence(), seqMapping.getTo(), map);
+ MapList cdsToDnaMap = new MapList(dnaExonRanges, exonRange, 1, 1);
+ newMapping.addMap(dnaSeq, exon.getDatasetSequence(), cdsToDnaMap);
+
+ exonSequences.add(exon);
+ }
+ return exonSequences;
+ }
}
git://source.jalview.org / jalview.git / blobdiff