From f28d892d6d2584e7eb44ff7333d49d60d787f706 Mon Sep 17 00:00:00 2001 From: gmungoc Date: Thu, 30 Jun 2016 15:49:51 +0100 Subject: [PATCH] JAL-2110 fixes and tests for synthesized EMBLCDSPROTEIN dbrefs --- src/jalview/datamodel/xdb/embl/EmblEntry.java | 222 +++++++++++--------- test/jalview/datamodel/xdb/embl/EmblEntryTest.java | 58 ++++- test/jalview/datamodel/xdb/embl/EmblFileTest.java | 21 +- .../jalview/datamodel/xdb/embl/EmblTestHelper.java | 8 + 4 files changed, 199 insertions(+), 110 deletions(-) diff --git a/src/jalview/datamodel/xdb/embl/EmblEntry.java b/src/jalview/datamodel/xdb/embl/EmblEntry.java index 5409d5b..a2354ed 100644 --- a/src/jalview/datamodel/xdb/embl/EmblEntry.java +++ b/src/jalview/datamodel/xdb/embl/EmblEntry.java @@ -254,11 +254,11 @@ public class EmblEntry { boolean isEmblCdna = sourceDb.equals(DBRefSource.EMBLCDS); - int[] exon = getCdsRanges(feature); + int[] exons = getCdsRanges(feature); - String prseq = null; - String prname = ""; - String prid = null; + String translation = null; + String proteinName = ""; + String proteinId = null; Map vals = new Hashtable(); /* @@ -279,11 +279,11 @@ public class EmblEntry if (qname.equals("translation")) { // remove all spaces (precompiled String.replaceAll(" ", "")) - prseq = SPACE_PATTERN.matcher(q.getValues()[0]).replaceAll(""); + translation = SPACE_PATTERN.matcher(q.getValues()[0]).replaceAll(""); } else if (qname.equals("protein_id")) { - prid = q.getValues()[0].trim(); + proteinId = q.getValues()[0].trim(); } else if (qname.equals("codon_start")) { @@ -299,7 +299,7 @@ public class EmblEntry else if (qname.equals("product")) { // sometimes name is returned e.g. for V00488 - prname = q.getValues()[0].trim(); + proteinName = q.getValues()[0].trim(); } else { @@ -315,53 +315,54 @@ public class EmblEntry } } - DBRefEntry protEMBLCDS = null; - exon = MappingUtils.removeStartPositions(codonStart - 1, exon); - boolean noProteinDbref = true; + DBRefEntry proteinToEmblProteinRef = null; + exons = MappingUtils.removeStartPositions(codonStart - 1, exons); SequenceI product = null; - Mapping map = null; - if (prseq != null && prname != null && prid != null) + Mapping dnaToProteinMapping = null; + if (translation != null && proteinName != null && proteinId != null) { /* * look for product in peptides list, if not found, add it */ - product = matcher.findIdMatch(prid); + product = matcher.findIdMatch(proteinId); if (product == null) { - product = new Sequence(prid, prseq, 1, prseq.length()); - product.setDescription(((prname.length() == 0) ? "Protein Product from " + product = new Sequence(proteinId, translation, 1, translation.length()); + product.setDescription(((proteinName.length() == 0) ? "Protein Product from " + sourceDb - : prname)); + : proteinName)); peptides.add(product); matcher.add(product); } // we have everything - create the mapping and perhaps the protein // sequence - if (exon == null || exon.length == 0) + if (exons == null || exons.length == 0) { System.err .println("Implementation Notice: EMBLCDS records not properly supported yet - Making up the CDNA region of this sequence... may be incorrect (" + sourceDb + ":" + getAccession() + ")"); - if (prseq.length() * 3 == (1 - codonStart + dna.getSequence().length)) + if (translation.length() * 3 == (1 - codonStart + dna.getSequence().length)) { System.err .println("Not allowing for additional stop codon at end of cDNA fragment... !"); // this might occur for CDS sequences where no features are // marked. - exon = new int[] { dna.getStart() + (codonStart - 1), + exons = new int[] { dna.getStart() + (codonStart - 1), dna.getEnd() }; - map = new Mapping(product, exon, new int[] { 1, prseq.length() }, + dnaToProteinMapping = new Mapping(product, exons, new int[] { 1, + translation.length() }, 3, 1); } - if ((prseq.length() + 1) * 3 == (1 - codonStart + dna.getSequence().length)) + if ((translation.length() + 1) * 3 == (1 - codonStart + dna.getSequence().length)) { System.err .println("Allowing for additional stop codon at end of cDNA fragment... will probably cause an error in VAMSAs!"); - exon = new int[] { dna.getStart() + (codonStart - 1), + exons = new int[] { dna.getStart() + (codonStart - 1), dna.getEnd() - 3 }; - map = new Mapping(product, exon, new int[] { 1, prseq.length() }, + dnaToProteinMapping = new Mapping(product, exons, new int[] { 1, + translation.length() }, 3, 1); } } @@ -381,37 +382,42 @@ public class EmblEntry else { // final product length truncation check - // TODO should from range include stop codon even if not in protein - // in order to include stop codon in CDS sequence (as done for - // Ensembl)? - int[] cdsRanges = adjustForProteinLength(prseq.length(), exon); - map = new Mapping(product, cdsRanges, new int[] { 1, - prseq.length() }, 3, 1); - // reconstruct the EMBLCDS entry - // TODO: this is only necessary when there codon annotation is - // complete (I think JBPNote) - DBRefEntry pcdnaref = new DBRefEntry(); - pcdnaref.setAccessionId(prid); - pcdnaref.setSource(DBRefSource.EMBLCDS); - pcdnaref.setVersion(getSequenceVersion()); // same as parent EMBL - // version. - MapList mp = new MapList(new int[] { 1, prseq.length() }, - new int[] { 1 + (codonStart - 1), - (codonStart - 1) + 3 * prseq.length() }, 1, 3); - pcdnaref.setMap(new Mapping(mp)); + int[] cdsRanges = adjustForProteinLength(translation.length(), exons); + dnaToProteinMapping = new Mapping(product, cdsRanges, new int[] { 1, + translation.length() }, 3, 1); if (product != null) { - product.addDBRef(pcdnaref); - protEMBLCDS = new DBRefEntry(pcdnaref); - protEMBLCDS.setSource(DBRefSource.EMBLCDSProduct); - product.addDBRef(protEMBLCDS); + /* + * make xrefs from protein to EMBLCDS and EMBLCDSPROTEIN + */ + DBRefEntry proteinToEmblCdsRef = new DBRefEntry(); + proteinToEmblCdsRef.setAccessionId(proteinId); + proteinToEmblCdsRef.setSource(DBRefSource.EMBLCDS); + proteinToEmblCdsRef.setVersion(getSequenceVersion()); // same as + // parent EMBL + // version. + MapList mp = new MapList(new int[] { 1, translation.length() }, + new int[] { 1 + (codonStart - 1), + (codonStart - 1) + 3 * translation.length() }, 1, 3); + proteinToEmblCdsRef.setMap(new Mapping(mp)); + product.addDBRef(proteinToEmblCdsRef); + proteinToEmblProteinRef = new DBRefEntry(proteinToEmblCdsRef); + MapList mp2 = new MapList( + new int[] { 1, translation.length() }, new int[] { 1, + translation.length() }, 1, 1); + proteinToEmblProteinRef.setMap(new Mapping(mp2)); + proteinToEmblProteinRef.setSource(DBRefSource.EMBLCDSProduct); + product.addDBRef(proteinToEmblProteinRef); } } } - // add cds feature to dna seq - this may include the stop codon - for (int xint = 0; exon != null && xint < exon.length; xint += 2) + + /* + * add cds features to dna sequence + */ + for (int xint = 0; exons != null && xint < exons.length; xint += 2) { - SequenceFeature sf = makeCdsFeature(exon, xint, prname, prid, vals, + SequenceFeature sf = makeCdsFeature(exons, xint, proteinName, proteinId, vals, codonStart); sf.setType(feature.getName()); // "CDS" sf.setEnaLocation(feature.getLocation()); @@ -421,8 +427,9 @@ public class EmblEntry } /* - * add dbRefs to sequence, and mappings for Uniprot xrefs + * add feature dbRefs to sequence, and mappings for Uniprot xrefs */ + boolean hasUniprotDbref = false; if (feature.dbRefs != null) { boolean mappingUsed = false; @@ -431,12 +438,15 @@ public class EmblEntry /* * ensure UniProtKB/Swiss-Prot converted to UNIPROT */ - ref.setSource(DBRefUtils.getCanonicalName(ref.getSource())); - if (ref.getSource().equals(DBRefSource.UNIPROT)) + String source = DBRefUtils.getCanonicalName(ref.getSource()); + ref.setSource(source); + DBRefEntry proteinToDnaRef = new DBRefEntry(ref.getSource(), ref.getVersion(), ref + .getAccessionId()); + if (source.equals(DBRefSource.UNIPROT)) { String proteinSeqName = DBRefSource.UNIPROT + "|" + ref.getAccessionId(); - if (map != null && map.getTo() != null) + if (dnaToProteinMapping != null && dnaToProteinMapping.getTo() != null) { if (mappingUsed) { @@ -444,13 +454,14 @@ public class EmblEntry * two or more Uniprot xrefs for the same CDS - * each needs a distinct Mapping (as to a different sequence) */ - map = new Mapping(map); + dnaToProteinMapping = new Mapping(dnaToProteinMapping); } mappingUsed = true; /* * try to locate the protein mapped to (possibly by a - * previous CDS feature) + * previous CDS feature); if not found, construct it from + * the EMBL translation */ SequenceI proteinSeq = matcher.findIdMatch(proteinSeqName); if (proteinSeq == null) @@ -460,61 +471,62 @@ public class EmblEntry matcher.add(proteinSeq); peptides.add(proteinSeq); } - map.setTo(proteinSeq); - map.getTo().addDBRef( - new DBRefEntry(ref.getSource(), ref.getVersion(), ref - .getAccessionId())); - ref.setMap(map); + dnaToProteinMapping.setTo(proteinSeq); + proteinSeq.addDBRef(proteinToDnaRef); + ref.setMap(dnaToProteinMapping); } - noProteinDbref = false; + hasUniprotDbref = true; } if (product != null) { - DBRefEntry pref = new DBRefEntry(ref.getSource(), - ref.getVersion(), ref.getAccessionId()); + /* + * copy feature dbref to our protein product + */ + DBRefEntry pref = proteinToDnaRef; pref.setMap(null); // reference is direct product.addDBRef(pref); // Add converse mapping reference - if (map != null) + if (dnaToProteinMapping != null) { - Mapping pmap = new Mapping(dna, map.getMap().getInverse()); + Mapping pmap = new Mapping(dna, dnaToProteinMapping.getMap() + .getInverse()); pref = new DBRefEntry(sourceDb, getSequenceVersion(), this.getAccession()); pref.setMap(pmap); - if (map.getTo() != null) + if (dnaToProteinMapping.getTo() != null) { - map.getTo().addDBRef(pref); + dnaToProteinMapping.getTo().addDBRef(pref); } } } dna.addDBRef(ref); } - if (noProteinDbref && product != null) + } + /* + * if we have a product (translation) but no explicit Uniprot dbref + * (example: EMBL AAFI02000057 protein_id EAL65544.1 + * construct mappings to an assumed EMBLCDSPROTEIN accession + */ + if (!hasUniprotDbref && product != null) + { + if (proteinToEmblProteinRef == null) { - // add protein coding reference to dna sequence so xref matches - if (protEMBLCDS == null) - { - protEMBLCDS = new DBRefEntry(); - protEMBLCDS.setAccessionId(prid); - protEMBLCDS.setSource(DBRefSource.EMBLCDSProduct); - protEMBLCDS.setVersion(getSequenceVersion()); - protEMBLCDS - .setMap(new Mapping(product, map.getMap().getInverse())); - } - product.addDBRef(protEMBLCDS); + proteinToEmblProteinRef = new DBRefEntry(); + proteinToEmblProteinRef.setAccessionId(proteinId); + proteinToEmblProteinRef.setSource(DBRefSource.EMBLCDSProduct); + proteinToEmblProteinRef.setVersion(getSequenceVersion()); + proteinToEmblProteinRef.setMap(new Mapping(product, + dnaToProteinMapping.getMap().getInverse())); + } + product.addDBRef(proteinToEmblProteinRef); - // Add converse mapping reference - if (map != null) - { - Mapping pmap = new Mapping(product, protEMBLCDS.getMap().getMap() - .getInverse()); - DBRefEntry ncMap = new DBRefEntry(protEMBLCDS); - ncMap.setMap(pmap); - if (map.getTo() != null) - { - dna.addDBRef(ncMap); - } - } + if (dnaToProteinMapping != null + && dnaToProteinMapping.getTo() != null) + { + DBRefEntry dnaToEmblProteinRef = new DBRefEntry( + proteinToEmblProteinRef); + dnaToEmblProteinRef.setMap(dnaToProteinMapping); + dna.addDBRef(dnaToEmblProteinRef); } } } @@ -615,26 +627,30 @@ public class EmblEntry } /** - * truncate the last exon interval to the prlength'th codon + * Truncates (if necessary) the exon intervals to match 3 times the length of + * the protein; also accepts 3 bases longer (for stop codon not included in + * protein) * - * @param prlength + * @param proteinLength * @param exon - * @return new exon + * an array of [start, end, start, end...] intervals + * @return the same array (if unchanged) or a truncated copy */ - static int[] adjustForProteinLength(int prlength, int[] exon) + static int[] adjustForProteinLength(int proteinLength, int[] exon) { - if (prlength <= 0 || exon == null) + if (proteinLength <= 0 || exon == null) { return exon; } - int desiredCdsLength = prlength * 3; + int expectedCdsLength = proteinLength * 3; int exonLength = MappingUtils.getLength(Arrays.asList(exon)); /* - * assuming here exon might include stop codon in addition to protein codons + * if exon length matches protein, or is shorter, or longer by the + * length of a stop codon (3 bases), then leave it unchanged */ - if (desiredCdsLength == exonLength - || desiredCdsLength == exonLength - 3) + if (expectedCdsLength >= exonLength + || expectedCdsLength == exonLength - 3) { return exon; } @@ -648,11 +664,11 @@ public class EmblEntry for (int x = 0; x < exon.length; x += 2) { cdspos += Math.abs(exon[x + 1] - exon[x]) + 1; - if (desiredCdsLength <= cdspos) + if (expectedCdsLength <= cdspos) { // advanced beyond last codon. sxpos = x; - if (desiredCdsLength != cdspos) + if (expectedCdsLength != cdspos) { // System.err // .println("Truncating final exon interval on region by " @@ -665,11 +681,11 @@ public class EmblEntry */ if (exon[x + 1] >= exon[x]) { - endxon = exon[x + 1] - cdspos + desiredCdsLength; + endxon = exon[x + 1] - cdspos + expectedCdsLength; } else { - endxon = exon[x + 1] + cdspos - desiredCdsLength; + endxon = exon[x + 1] + cdspos - expectedCdsLength; } break; } diff --git a/test/jalview/datamodel/xdb/embl/EmblEntryTest.java b/test/jalview/datamodel/xdb/embl/EmblEntryTest.java index 3de5e3f..d987e53 100644 --- a/test/jalview/datamodel/xdb/embl/EmblEntryTest.java +++ b/test/jalview/datamodel/xdb/embl/EmblEntryTest.java @@ -42,8 +42,7 @@ public class EmblEntryTest EmblFile ef = EmblTestHelper.getEmblFile(); /* - * parse two CDS features, one with two Uniprot cross-refs, - * the other with one + * parse three CDS features, with two/one/no Uniprot cross-refs */ EmblEntry testee = new EmblEntry(); for (EmblFeature feature : ef.getEntries().get(0).getFeatures()) @@ -57,9 +56,10 @@ public class EmblEntryTest /* * peptides should now have five entries: * EMBL product and two Uniprot accessions for the first CDS / translation - * EMBL product and one Uniprot accession for the second CDS / translation + * EMBL product and one Uniprot accession for the second CDS / " + * EMBL product and synthesized EMBLCDSPROTEINaccession for the third */ - assertEquals(5, peptides.size()); + assertEquals(6, peptides.size()); assertEquals("CAA30420.1", peptides.get(0).getName()); assertEquals("MLCF", peptides.get(0).getSequenceAsString()); assertEquals("UNIPROT|B0BCM4", peptides.get(1).getName()); @@ -70,12 +70,14 @@ public class EmblEntryTest assertEquals("MSSS", peptides.get(3).getSequenceAsString()); assertEquals("UNIPROT|B0BCM3", peptides.get(4).getName()); assertEquals("MSSS", peptides.get(4).getSequenceAsString()); + assertEquals("CAA12345.6", peptides.get(5).getName()); + assertEquals("MSS", peptides.get(5).getSequenceAsString()); /* * verify dna sequence has dbrefs with mappings to the peptide 'products' */ DBRefEntry[] dbrefs = dna.getDBRefs(); - assertEquals(3, dbrefs.length); + assertEquals(4, dbrefs.length); DBRefEntry dbRefEntry = dbrefs[0]; assertEquals("UNIPROT", dbRefEntry.getSource()); assertEquals("B0BCM4", dbRefEntry.getAccessionId()); @@ -114,5 +116,51 @@ public class EmblEntryTest assertEquals(1, toRanges.size()); assertEquals(1, toRanges.get(0)[0]); assertEquals(4, toRanges.get(0)[1]); + + dbRefEntry = dbrefs[3]; + assertEquals("EMBLCDSPROTEIN", dbRefEntry.getSource()); + assertEquals("CAA12345.6", dbRefEntry.getAccessionId()); + assertSame(peptides.get(5), dbRefEntry.getMap().getTo()); + fromRanges = dbRefEntry.getMap().getMap().getFromRanges(); + assertEquals(2, fromRanges.size()); + assertEquals(4, fromRanges.get(0)[0]); + assertEquals(6, fromRanges.get(0)[1]); + assertEquals(10, fromRanges.get(1)[0]); + assertEquals(15, fromRanges.get(1)[1]); + toRanges = dbRefEntry.getMap().getMap().getToRanges(); + assertEquals(1, toRanges.size()); + assertEquals(1, toRanges.get(0)[0]); + assertEquals(3, toRanges.get(0)[1]); + } + + @Test(groups = "Functional") + public void testAdjustForProteinLength() + { + int[] exons = new int[] { 11, 15, 21, 25, 31, 38 }; // 18 bp + + // exact length match: + assertSame(exons, EmblEntry.adjustForProteinLength(6, exons)); + + // match if we assume exons include stop codon not in protein: + assertSame(exons, EmblEntry.adjustForProteinLength(5, exons)); + + // truncate last exon by 6bp + int[] truncated = EmblEntry.adjustForProteinLength(4, exons); + assertEquals("[11, 15, 21, 25, 31, 32]", + Arrays.toString(truncated)); + + // remove last exon and truncate preceding by 1bp + truncated = EmblEntry.adjustForProteinLength(3, exons); + assertEquals("[11, 15, 21, 24]", Arrays.toString(truncated)); + + // exact removal of exon case: + exons = new int[] { 11, 15, 21, 27, 33, 38 }; // 18 bp + truncated = EmblEntry.adjustForProteinLength(4, exons); + assertEquals("[11, 15, 21, 27]", Arrays.toString(truncated)); + + // what if exons are too short for protein? + truncated = EmblEntry.adjustForProteinLength(7, exons); + assertSame(exons, truncated); + // assertEquals("[11, 15, 21, 27]", Arrays.toString(truncated)); } } diff --git a/test/jalview/datamodel/xdb/embl/EmblFileTest.java b/test/jalview/datamodel/xdb/embl/EmblFileTest.java index 6955833..906436f 100644 --- a/test/jalview/datamodel/xdb/embl/EmblFileTest.java +++ b/test/jalview/datamodel/xdb/embl/EmblFileTest.java @@ -80,9 +80,9 @@ public class EmblFileTest assertEquals("0", dbref.getVersion()); /* - * two sequence features for CDS + * three sequence features for CDS */ - assertEquals(2, entry.getFeatures().size()); + assertEquals(3, entry.getFeatures().size()); /* * first CDS */ @@ -140,6 +140,23 @@ public class EmblFileTest assertEquals("MSSS", q.getValues()[0]); /* + * third CDS + */ + ef = entry.getFeatures().get(2); + assertEquals("CDS", ef.getName()); + assertEquals("join(4..6,10..15)", ef.getLocation()); + assertNull(ef.getDbRefs()); + assertEquals(2, ef.getQualifiers().size()); + q = ef.getQualifiers().get(0); + assertEquals("protein_id", q.getName()); + assertEquals(1, q.getValues().length); + assertEquals("CAA12345.6", q.getValues()[0]); + q = ef.getQualifiers().get(1); + assertEquals("translation", q.getName()); + assertEquals(1, q.getValues().length); + assertEquals("MSS", q.getValues()[0]); + + /* * Sequence - verify newline not converted to space (JAL-2029) */ EmblSequence seq = entry.getSequence(); diff --git a/test/jalview/datamodel/xdb/embl/EmblTestHelper.java b/test/jalview/datamodel/xdb/embl/EmblTestHelper.java index 9957c72..6349164 100644 --- a/test/jalview/datamodel/xdb/embl/EmblTestHelper.java +++ b/test/jalview/datamodel/xdb/embl/EmblTestHelper.java @@ -38,6 +38,14 @@ public class EmblTestHelper + "MSSS" + "" /* + * third CDS is made up - has no xref - code should synthesize + * one to an assumed EMBLCDSPROTEIN accession + */ + + "" + + "CAA12345.6" + + "MSS" + + "" + /* * sequence (modified for test purposes) * emulates EMBL XML 1.2 which splits sequence data every 60 characters * see EmblSequence.setSequence -- 1.7.10.2