X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignmentUtils.java;h=bef667dfb140e63241fe91059d1d0cac4bd78cb9;hb=d2fbfb17ae43c983123e18a0ef1074c2aba907ff;hp=311a9a3afbd352cf2dde74711f790a91065ecdc9;hpb=5fba3a87f8eaa121a8139f37e6822c15283e1a2c;p=jalview.git diff --git a/src/jalview/analysis/AlignmentUtils.java b/src/jalview/analysis/AlignmentUtils.java index 311a9a3..bef667d 100644 --- a/src/jalview/analysis/AlignmentUtils.java +++ b/src/jalview/analysis/AlignmentUtils.java @@ -29,6 +29,7 @@ import jalview.datamodel.Alignment; import jalview.datamodel.AlignmentAnnotation; import jalview.datamodel.AlignmentI; import jalview.datamodel.DBRefEntry; +import jalview.datamodel.GeneLociI; import jalview.datamodel.IncompleteCodonException; import jalview.datamodel.Mapping; import jalview.datamodel.Sequence; @@ -36,6 +37,7 @@ import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceGroup; import jalview.datamodel.SequenceI; import jalview.datamodel.features.SequenceFeatures; +import jalview.io.gff.Gff3Helper; import jalview.io.gff.SequenceOntologyI; import jalview.schemes.ResidueProperties; import jalview.util.Comparison; @@ -105,6 +107,15 @@ public class AlignmentUtils { return variant == null ? null : variant.getFeatureGroup(); } + + /** + * toString for aid in the debugger only + */ + @Override + public String toString() + { + return base + ":" + (variant == null ? "" : variant.getDescription()); + } } /** @@ -384,7 +395,7 @@ public class AlignmentUtils * Answers true if the mappings include one between the given (dataset) * sequences. */ - public static boolean mappingExists(List mappings, + protected static boolean mappingExists(List mappings, SequenceI aaSeq, SequenceI cdnaSeq) { if (mappings != null) @@ -1636,8 +1647,8 @@ public class AlignmentUtils productSeqs = new HashSet(); for (SequenceI seq : products) { - productSeqs.add(seq.getDatasetSequence() == null ? seq - : seq.getDatasetSequence()); + productSeqs.add(seq.getDatasetSequence() == null ? seq : seq + .getDatasetSequence()); } } @@ -1730,9 +1741,8 @@ public class AlignmentUtils /* * add a mapping from CDS to the (unchanged) mapped to range */ - List cdsRange = Collections - .singletonList(new int[] - { 1, cdsSeq.getLength() }); + List cdsRange = Collections.singletonList(new int[] { 1, + cdsSeq.getLength() }); MapList cdsToProteinMap = new MapList(cdsRange, mapList.getToRanges(), mapList.getFromRatio(), mapList.getToRatio()); @@ -1754,7 +1764,7 @@ public class AlignmentUtils * add another mapping from original 'from' range to CDS */ AlignedCodonFrame dnaToCdsMapping = new AlignedCodonFrame(); - MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), + final MapList dnaToCdsMap = new MapList(mapList.getFromRanges(), cdsRange, 1, 1); dnaToCdsMapping.addMap(dnaSeq.getDatasetSequence(), cdsSeqDss, dnaToCdsMap); @@ -1764,6 +1774,13 @@ public class AlignmentUtils } /* + * transfer dna chromosomal loci (if known) to the CDS + * sequence (via the mapping) + */ + final MapList cdsToDnaMap = dnaToCdsMap.getInverse(); + transferGeneLoci(dnaSeq, cdsToDnaMap, cdsSeq); + + /* * 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 @@ -1782,26 +1799,30 @@ public class AlignmentUtils 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))); + /* + * create a cross-reference from CDS to the source sequence's + * primary reference and vice versa + */ + String source = primRef.getSource(); + String version = primRef.getVersion(); + DBRefEntry cdsCrossRef = new DBRefEntry(source, source + ":" + + version, primRef.getAccessionId()); + cdsCrossRef.setMap(new Mapping(dnaDss, new MapList(cdsToDnaMap))); cdsSeqDss.addDBRef(cdsCrossRef); + dnaSeq.addDBRef(new DBRefEntry(source, version, cdsSeq + .getName(), new Mapping(cdsSeqDss, dnaToCdsMap))); + // 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()); + DBRefEntry proteinToCdsRef = new DBRefEntry(source, version, + cdsSeq.getName()); // - proteinToCdsRef.setMap( - new Mapping(cdsSeqDss, cdsToProteinMap.getInverse())); + proteinToCdsRef.setMap(new Mapping(cdsSeqDss, cdsToProteinMap + .getInverse())); proteinProduct.addDBRef(proteinToCdsRef); } @@ -1814,14 +1835,46 @@ public class AlignmentUtils } } - AlignmentI cds = new Alignment( - cdsSeqs.toArray(new SequenceI[cdsSeqs.size()])); + AlignmentI cds = new Alignment(cdsSeqs.toArray(new SequenceI[cdsSeqs + .size()])); cds.setDataset(dataset); return cds; } /** + * Tries to transfer gene loci (dbref to chromosome positions) from fromSeq to + * toSeq, mediated by the given mapping between the sequences + * + * @param fromSeq + * @param targetToFrom + * Map + * @param targetSeq + */ + protected static void transferGeneLoci(SequenceI fromSeq, + MapList targetToFrom, SequenceI targetSeq) + { + if (targetSeq.getGeneLoci() != null) + { + // already have - don't override + return; + } + GeneLociI fromLoci = fromSeq.getGeneLoci(); + if (fromLoci == null) + { + return; + } + + MapList newMap = targetToFrom.traverse(fromLoci.getMap()); + + if (newMap != null) + { + targetSeq.setGeneLoci(fromLoci.getSpeciesId(), + fromLoci.getAssemblyId(), fromLoci.getChromosomeId(), newMap); + } + } + + /** * 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. @@ -1989,19 +2042,19 @@ public class AlignmentUtils } /** - * add any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to + * Adds any DBRefEntrys to cdsSeq from contig that have a Mapping congruent to * the given mapping. * * @param cdsSeq * @param contig + * @param proteinProduct * @param mapping - * @return list of DBRefEntrys added. + * @return list of DBRefEntrys added */ - public static List propagateDBRefsToCDS(SequenceI cdsSeq, + protected static List propagateDBRefsToCDS(SequenceI cdsSeq, SequenceI contig, SequenceI proteinProduct, Mapping mapping) { - - // gather direct refs from contig congrent with mapping + // gather direct refs from contig congruent with mapping List direct = new ArrayList(); HashSet directSources = new HashSet(); if (contig.getDBRefs() != null) @@ -2081,7 +2134,7 @@ public class AlignmentUtils * subtypes in the Sequence Ontology) * @param omitting */ - public static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, + protected static int transferFeatures(SequenceI fromSeq, SequenceI toSeq, MapList mapping, String select, String... omitting) { SequenceI copyTo = toSeq; @@ -2164,7 +2217,10 @@ public class AlignmentUtils /** * Returns a mapping from dna to protein by inspecting sequence features of - * type "CDS" on the dna. + * type "CDS" on the dna. A mapping is constructed if the total CDS feature + * length is 3 times the peptide length (optionally after dropping a trailing + * stop codon). This method does not check whether the CDS nucleotide sequence + * translates to the peptide sequence. * * @param dnaSeq * @param proteinSeq @@ -2176,6 +2232,15 @@ public class AlignmentUtils List ranges = findCdsPositions(dnaSeq); int mappedDnaLength = MappingUtils.getLength(ranges); + /* + * if not a whole number of codons, something is wrong, + * abort mapping + */ + if (mappedDnaLength % CODON_LENGTH > 0) + { + return null; + } + int proteinLength = proteinSeq.getLength(); int proteinStart = proteinSeq.getStart(); int proteinEnd = proteinSeq.getEnd(); @@ -2199,6 +2264,7 @@ public class AlignmentUtils if (codesForResidues == (proteinLength + 1)) { // assuming extra codon is for STOP and not in peptide + // todo: check trailing codon is indeed a STOP codon codesForResidues--; mappedDnaLength -= CODON_LENGTH; MappingUtils.removeEndPositions(CODON_LENGTH, ranges); @@ -2222,7 +2288,7 @@ public class AlignmentUtils * @param dnaSeq * @return */ - public static List findCdsPositions(SequenceI dnaSeq) + protected static List findCdsPositions(SequenceI dnaSeq) { List result = new ArrayList(); @@ -2233,7 +2299,6 @@ public class AlignmentUtils return result; } SequenceFeatures.sortFeatures(sfs, true); - int startPhase = 0; for (SequenceFeature sf : sfs) { @@ -2251,7 +2316,7 @@ public class AlignmentUtils */ int begin = sf.getBegin(); int end = sf.getEnd(); - if (result.isEmpty()) + if (result.isEmpty() && phase > 0) { begin += phase; if (begin > end) @@ -2266,16 +2331,6 @@ public class AlignmentUtils } /* - * remove 'startPhase' positions (usually 0) from the first range - * so we begin at the start of a complete codon - */ - if (!result.isEmpty()) - { - // TODO JAL-2022 correctly model start phase > 0 - result.get(0)[0] += startPhase; - } - - /* * Finally sort ranges by start position. This avoids a dependency on * keeping features in order on the sequence (if they are in order anyway, * the sort will have almost no work to do). The implicit assumption is CDS @@ -2368,7 +2423,7 @@ public class AlignmentUtils { if (var.variant != null) { - String alleles = (String) var.variant.getValue("alleles"); + String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES); if (alleles != null) { for (String base : alleles.split(",")) @@ -2390,7 +2445,7 @@ public class AlignmentUtils { if (var.variant != null) { - String alleles = (String) var.variant.getValue("alleles"); + String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES); if (alleles != null) { for (String base : alleles.split(",")) @@ -2412,7 +2467,7 @@ public class AlignmentUtils { if (var.variant != null) { - String alleles = (String) var.variant.getValue("alleles"); + String alleles = (String) var.variant.getValue(Gff3Helper.ALLELES); if (alleles != null) { for (String base : alleles.split(",")) @@ -2507,7 +2562,10 @@ public class AlignmentUtils /** * 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 + * list of the base and all variants for each corresponding codon position. + *

+ * This depends on dna variants being held as a comma-separated list as + * property "alleles" on variant features. * * @param dnaSeq * @param dnaToProtein @@ -2545,6 +2603,30 @@ public class AlignmentUtils // 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) { @@ -2563,21 +2645,6 @@ public class AlignmentUtils } /* - * 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