2 * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3 * Copyright (C) $$Year-Rel$$ The Jalview Authors
5 * This file is part of Jalview.
7 * Jalview is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation, either version 3
10 * of the License, or (at your option) any later version.
12 * Jalview is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty
14 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 * PURPOSE. See the GNU General Public License for more details.
17 * You should have received a copy of the GNU General Public License
18 * along with Jalview. If not, see <http://www.gnu.org/licenses/>.
19 * The Jalview Authors are detailed in the 'AUTHORS' file.
23 import java.io.IOException;
24 import java.text.ParseException;
25 import java.util.ArrayList;
26 import java.util.Arrays;
27 import java.util.HashMap;
28 import java.util.Hashtable;
29 import java.util.List;
30 import java.util.Locale;
32 import java.util.Map.Entry;
33 import java.util.TreeMap;
35 import jalview.bin.Console;
36 import jalview.datamodel.DBRefEntry;
37 import jalview.datamodel.DBRefSource;
38 import jalview.datamodel.FeatureProperties;
39 import jalview.datamodel.Mapping;
40 import jalview.datamodel.Sequence;
41 import jalview.datamodel.SequenceFeature;
42 import jalview.datamodel.SequenceI;
43 import jalview.util.DBRefUtils;
44 import jalview.util.DnaUtils;
45 import jalview.util.MapList;
46 import jalview.util.MappingUtils;
49 * A base class to support parsing of GenBank, EMBL or DDBJ flat file format
50 * data. Example files (rather than formal specifications) are provided at
53 * https://ena-docs.readthedocs.io/en/latest/submit/fileprep/flat-file-example.html
54 * https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html
57 * or to compare the same entry, see
60 * https://www.ebi.ac.uk/ena/browser/api/embl/X81322.1
61 * https://www.ncbi.nlm.nih.gov/nuccore/X81322.1
64 * The feature table part of the file has a common definition, only the start of
65 * each line is formatted differently in GenBank and EMBL. See
66 * http://www.insdc.org/files/feature_table.html#7.1.
68 public abstract class EMBLLikeFlatFile extends AlignFile
70 protected static final String LOCATION = "location";
72 protected static final String QUOTE = "\"";
74 protected static final String DOUBLED_QUOTE = QUOTE + QUOTE;
76 protected static final String WHITESPACE = "\\s+";
79 * Removes leading or trailing double quotes (") unless doubled, and changes
80 * any 'escaped' (doubled) double quotes to single characters. As per the
81 * Feature Table specification for Qualifiers, Free Text.
86 protected static String removeQuotes(String value)
92 if (value.startsWith(QUOTE) && !value.startsWith(DOUBLED_QUOTE))
94 value = value.substring(1);
96 if (value.endsWith(QUOTE) && !value.endsWith(DOUBLED_QUOTE))
98 value = value.substring(0, value.length() - 1);
100 value = value.replace(DOUBLED_QUOTE, QUOTE);
105 * Truncates (if necessary) the exon intervals to match 3 times the length of
106 * the protein(including truncation for stop codon included in exon)
108 * @param proteinLength
110 * an array of [start, end, start, end...] intervals
111 * @return the same array (if unchanged) or a truncated copy
113 protected static int[] adjustForProteinLength(int proteinLength,
116 if (proteinLength <= 0 || exon == null)
120 int expectedCdsLength = proteinLength * 3;
121 int exonLength = MappingUtils.getLength(Arrays.asList(exon));
124 * if exon length matches protein, or is shorter, then leave it unchanged
126 if (expectedCdsLength >= exonLength)
134 origxon = new int[exon.length];
135 System.arraycopy(exon, 0, origxon, 0, exon.length);
137 for (int x = 0; x < exon.length; x += 2)
139 cdspos += Math.abs(exon[x + 1] - exon[x]) + 1;
140 if (expectedCdsLength <= cdspos)
142 // advanced beyond last codon.
144 if (expectedCdsLength != cdspos)
147 // .println("Truncating final exon interval on region by "
148 // + (cdspos - cdslength));
152 * shrink the final exon - reduce end position if forward
153 * strand, increase it if reverse
155 if (exon[x + 1] >= exon[x])
157 endxon = exon[x + 1] - cdspos + expectedCdsLength;
161 endxon = exon[x + 1] + cdspos - expectedCdsLength;
169 // and trim the exon interval set if necessary
170 int[] nxon = new int[sxpos + 2];
171 System.arraycopy(exon, 0, nxon, 0, sxpos + 2);
172 nxon[sxpos + 1] = endxon; // update the end boundary for the new exon
180 * when true, interpret the mol_type 'source' feature attribute
181 * and generate an RNA sequence from the DNA record
183 protected boolean produceRna = true;
186 * values parsed from the data file
188 protected String sourceDb;
190 protected String accession;
192 protected String version;
194 protected String description;
196 protected int length = 128;
198 protected List<DBRefEntry> dbrefs;
200 protected boolean sequenceStringIsRNA = false;
202 protected String sequenceString;
204 protected Map<String, CdsData> cds;
211 * @throws IOException
213 public EMBLLikeFlatFile(FileParse fp, String sourceId) throws IOException
215 super(false, fp); // don't parse immediately
216 this.sourceDb = sourceId;
217 dbrefs = new ArrayList<>();
220 * using TreeMap gives CDS sequences in alphabetical, so readable, order
222 cds = new TreeMap<>(String.CASE_INSENSITIVE_ORDER);
228 * process attributes for 'source' until the next FT feature entry only
229 * interested in 'mol_type'
233 * @throws IOException
235 private String parseSourceQualifiers(String[] tokens) throws IOException
237 if (!"source".equals(tokens[0]))
239 throw (new RuntimeException("Not given a 'source' qualifier line"));
241 // search for mol_type attribute
243 StringBuilder sb = new StringBuilder().append(tokens[1]); // extent of
246 String line = parseFeatureQualifier(sb, false);
249 if (!line.startsWith("FT ")) // four spaces, end of this feature table
256 int p = line.indexOf("\\mol_type");
257 int qs = line.indexOf("\"", p);
258 int qe = line.indexOf("\"", qs + 1);
259 String qualifier = line.substring(qs, qe).toLowerCase(Locale.ROOT);
260 if (qualifier.indexOf("rna") > -1)
262 sequenceStringIsRNA = true;
264 if (qualifier.indexOf("dna") > -1)
266 sequenceStringIsRNA = false;
268 line = parseFeatureQualifier(sb, false);
274 * Parses one (GenBank or EMBL format) CDS feature, saves the parsed data, and
275 * returns the next line
279 * @throws IOException
281 protected String parseCDSFeature(String location) throws IOException
286 * parse location, which can be over >1 line e.g. EAW51554
288 CdsData data = new CdsData();
289 StringBuilder sb = new StringBuilder().append(location);
290 line = parseFeatureQualifier(sb, false);
291 data.cdsLocation = sb.toString();
295 if (!isFeatureContinuationLine(line))
297 // e.g. start of next feature "FT source..."
302 * extract qualifier, e.g. FT /protein_id="CAA37824.1"
303 * - the value may extend over more than one line
304 * - if the value has enclosing quotes, these are removed
305 * - escaped double quotes ("") are reduced to a single character
307 int slashPos = line.indexOf('/');
310 Console.error("Unexpected EMBL line ignored: " + line);
314 int eqPos = line.indexOf('=', slashPos + 1);
317 // can happen, e.g. /ribosomal_slippage
321 String qualifier = line.substring(slashPos + 1, eqPos);
322 String value = line.substring(eqPos + 1);
323 value = removeQuotes(value);
324 sb = new StringBuilder().append(value);
325 boolean asText = !"translation".equals(qualifier);
326 line = parseFeatureQualifier(sb, asText);
327 String featureValue = sb.toString();
329 if ("protein_id".equals(qualifier))
331 data.proteinId = featureValue;
333 else if ("codon_start".equals(qualifier))
337 data.codonStart = Integer.parseInt(featureValue.trim());
338 } catch (NumberFormatException e)
340 Console.error("Invalid codon_start in XML for " + this.accession
341 + ": " + e.getMessage());
344 else if ("db_xref".equals(qualifier))
346 String[] parts = featureValue.split(":");
347 if (parts.length == 2)
349 String db = parts[0].trim();
350 db = DBRefUtils.getCanonicalName(db);
351 DBRefEntry dbref = new DBRefEntry(db, "0", parts[1].trim());
352 data.xrefs.add(dbref);
355 else if ("product".equals(qualifier))
357 data.proteinName = featureValue;
359 else if ("translation".equals(qualifier))
361 data.translation = featureValue;
363 else if (!"".equals(featureValue))
365 // throw anything else into the additional properties hash
366 data.cdsProps.put(qualifier, featureValue);
370 if (data.proteinId != null)
372 this.cds.put(data.proteinId, data);
376 Console.error("Ignoring CDS feature with no protein_id for "
377 + sourceDb + ":" + accession);
383 protected abstract boolean isFeatureContinuationLine(String line);
386 * Output (print) is not (yet) implemented for flat file format
389 public String print(SequenceI[] seqs, boolean jvsuffix)
395 * Constructs and saves the sequence from parsed components
397 protected void buildSequence()
399 if (this.accession == null || this.sequenceString == null)
401 Console.error("Failed to parse data from EMBL");
405 String name = this.accession;
406 if (this.sourceDb != null)
408 name = this.sourceDb + "|" + name;
411 if (produceRna && sequenceStringIsRNA)
413 sequenceString = sequenceString.replace('T', 'U').replace('t', 'u');
416 SequenceI seq = new Sequence(name, this.sequenceString);
417 seq.setDescription(this.description);
420 * add a DBRef to itself
422 DBRefEntry selfRef = new DBRefEntry(sourceDb, version, accession);
423 int[] startEnd = new int[] { 1, seq.getLength() };
424 selfRef.setMap(new Mapping(null, startEnd, startEnd, 1, 1));
425 seq.addDBRef(selfRef);
427 for (DBRefEntry dbref : this.dbrefs)
432 processCDSFeatures(seq);
434 seq.deriveSequence();
440 * Process the CDS features, including generation of cross-references and
441 * mappings to the protein products (translation)
445 protected void processCDSFeatures(SequenceI seq)
448 * record protein products found to avoid duplication i.e. >1 CDS with
449 * the same /protein_id [though not sure I can find an example of this]
451 Map<String, SequenceI> proteins = new HashMap<>();
452 for (CdsData data : cds.values())
454 processCDSFeature(seq, data, proteins);
459 * Processes data for one parsed CDS feature to
461 * <li>create a protein product sequence for the translation</li>
462 * <li>create a cross-reference to protein with mapping from dna</li>
463 * <li>add a CDS feature to the sequence for each CDS start-end range</li>
464 * <li>add any CDS dbrefs to the sequence and to the protein product</li>
470 * map of protein products so far derived from CDS data
472 void processCDSFeature(SequenceI dna, CdsData data,
473 Map<String, SequenceI> proteins)
476 * parse location into a list of [start, end, start, end] positions
478 int[] exons = getCdsRanges(this.accession, data.cdsLocation);
480 MapList maplist = buildMappingToProtein(dna, exons, data);
484 for (int xint = 0; exons != null && xint < exons.length - 1; xint += 2)
486 int exonStart = exons[xint];
487 int exonEnd = exons[xint + 1];
488 int begin = Math.min(exonStart, exonEnd);
489 int end = Math.max(exonStart, exonEnd);
491 String desc = String.format("Exon %d for protein EMBLCDS:%s",
492 exonNumber, data.proteinId);
494 SequenceFeature sf = new SequenceFeature("CDS", desc, begin, end,
496 for (Entry<String, String> val : data.cdsProps.entrySet())
498 sf.setValue(val.getKey(), val.getValue());
501 sf.setEnaLocation(data.cdsLocation);
502 boolean forwardStrand = exonStart <= exonEnd;
503 sf.setStrand(forwardStrand ? "+" : "-");
504 sf.setPhase(String.valueOf(data.codonStart - 1));
505 sf.setValue(FeatureProperties.EXONPOS, exonNumber);
506 sf.setValue(FeatureProperties.EXONPRODUCT, data.proteinName);
508 dna.addSequenceFeature(sf);
511 boolean hasUniprotDbref = false;
512 for (DBRefEntry xref : data.xrefs)
515 if (xref.getSource().equals(DBRefSource.UNIPROT))
518 * construct (or find) the sequence for (data.protein_id, data.translation)
520 SequenceI protein = buildProteinProduct(dna, xref, data, proteins);
521 Mapping map = new Mapping(protein, maplist);
522 map.setMappedFromId(data.proteinId);
526 * add DBRefs with mappings from dna to protein and the inverse
528 DBRefEntry db1 = new DBRefEntry(sourceDb, version, accession);
529 db1.setMap(new Mapping(dna, maplist.getInverse()));
530 protein.addDBRef(db1);
532 hasUniprotDbref = true;
537 * if we have a product (translation) but no explicit Uniprot dbref
538 * (example: EMBL M19487 protein_id AAB02592.1)
539 * then construct mappings to an assumed EMBLCDSPROTEIN accession
541 if (!hasUniprotDbref)
543 SequenceI protein = proteins.get(data.proteinId);
546 protein = new Sequence(data.proteinId, data.translation);
547 protein.setDescription(data.proteinName);
548 proteins.put(data.proteinId, protein);
550 // assuming CDSPROTEIN sequence version = dna version (?!)
551 DBRefEntry db1 = new DBRefEntry(DBRefSource.EMBLCDSProduct,
552 this.version, data.proteinId);
553 protein.addDBRef(db1);
555 DBRefEntry dnaToEmblProteinRef = new DBRefEntry(
556 DBRefSource.EMBLCDSProduct, this.version, data.proteinId);
557 Mapping map = new Mapping(protein, maplist);
558 map.setMappedFromId(data.proteinId);
559 dnaToEmblProteinRef.setMap(map);
560 dna.addDBRef(dnaToEmblProteinRef);
564 * comment brought forward from EmblXmlSource, lines 447-451:
565 * TODO: if retrieved from EMBLCDS, add a DBRef back to the parent EMBL
566 * sequence with the exon map; if given a dataset reference, search
567 * dataset for parent EMBL sequence if it exists and set its map;
568 * make a new feature annotating the coding contig
573 * Computes a mapping from CDS positions in DNA sequence to protein product
574 * positions, with allowance for stop codon or incomplete start codon
581 MapList buildMappingToProtein(final SequenceI dna, final int[] exons,
584 MapList dnaToProteinMapping = null;
585 int peptideLength = data.translation.length();
587 int[] proteinRange = new int[] { 1, peptideLength };
588 if (exons != null && exons.length > 0)
591 * We were able to parse 'location'; do a final
592 * product length truncation check
594 int[] cdsRanges = adjustForProteinLength(peptideLength, exons);
595 dnaToProteinMapping = new MapList(cdsRanges, proteinRange, 3, 1);
600 * workaround until we handle all 'location' formats fully
601 * e.g. X53828.1:60..1058 or <123..>289
603 Console.error(String.format(
604 "Implementation Notice: EMBLCDS location '%s'not properly supported yet"
605 + " - Making up the CDNA region of (%s:%s)... may be incorrect",
606 data.cdsLocation, sourceDb, this.accession));
608 int completeCodonsLength = 1 - data.codonStart + dna.getLength();
609 int mappedDnaEnd = dna.getEnd();
610 if (peptideLength * 3 == completeCodonsLength)
612 // this might occur for CDS sequences where no features are marked
613 Console.warn("Assuming no stop codon at end of cDNA fragment");
614 mappedDnaEnd = dna.getEnd();
616 else if ((peptideLength + 1) * 3 == completeCodonsLength)
618 Console.warn("Assuming stop codon at end of cDNA fragment");
619 mappedDnaEnd = dna.getEnd() - 3;
622 if (mappedDnaEnd != -1)
624 int[] cdsRanges = new int[] {
625 dna.getStart() + (data.codonStart - 1), mappedDnaEnd };
626 dnaToProteinMapping = new MapList(cdsRanges, proteinRange, 3, 1);
630 return dnaToProteinMapping;
634 * Constructs a sequence for the protein product for the CDS data (if there is
635 * one), and dbrefs with mappings from CDS to protein and the reverse
643 SequenceI buildProteinProduct(SequenceI dna, DBRefEntry xref,
644 CdsData data, Map<String, SequenceI> proteins)
647 * check we have some data to work with
649 if (data.proteinId == null || data.translation == null)
655 * Construct the protein sequence (if not already seen)
657 String proteinSeqName = xref.getSource() + "|" + xref.getAccessionId();
658 SequenceI protein = proteins.get(proteinSeqName);
661 protein = new Sequence(proteinSeqName, data.translation, 1,
662 data.translation.length());
663 protein.setDescription(data.proteinName != null ? data.proteinName
664 : "Protein Product from " + sourceDb);
665 proteins.put(proteinSeqName, protein);
672 * Returns the CDS location as a single array of [start, end, start, end...]
673 * positions. If on the reverse strand, these will be in descending order.
679 protected int[] getCdsRanges(String accession, String location)
681 if (location == null)
688 List<int[]> ranges = DnaUtils.parseLocation(location);
689 return MappingUtils.rangeListToArray(ranges);
690 } catch (ParseException e)
693 String.format("Not parsing inexact CDS location %s in ENA %s",
694 location, accession));
700 * Reads the value of a feature (FT) qualifier from one or more lines of the
701 * file, and returns the next line after that. Values are appended to the
702 * string buffer, which should be already primed with the value read from the
703 * first line for the qualifier (with any leading double quote removed).
704 * Enclosing double quotes are removed, and escaped (repeated) double quotes
705 * reduced to one only. For example for
708 * FT /note="gene_id=hCG28070.3
709 * FT ""foobar"" isoform=CRA_b"
710 * the returned value is
711 * gene_id=hCG28070.3 "foobar" isoform=CRA_b
714 * Note the side-effect of this method, to advance data reading to the next
715 * line after the feature qualifier (which could be another qualifier, a
716 * different feature, a non-feature line, or null at end of file).
719 * a string buffer primed with the first line of the value
722 * @throws IOException
724 String parseFeatureQualifier(StringBuilder sb, boolean asText)
728 while ((line = nextLine()) != null)
730 if (!isFeatureContinuationLine(line))
732 break; // reached next feature or other input line
734 String[] tokens = line.split(WHITESPACE);
735 if (tokens.length < 2)
737 Console.error("Ignoring bad EMBL line for " + this.accession + ": "
741 if (tokens[1].startsWith("/"))
743 break; // next feature qualifier
747 * if text (e.g. /product), add a word separator for a new line,
748 * else (e.g. /translation) don't
756 * remove trailing " and unescape doubled ""
758 String data = removeQuotes(tokens[1]);
766 * Reads and saves the sequence, read from the lines following the ORIGIN
767 * (GenBank) or SQ (EMBL) line. Whitespace and position counters are
768 * discarded. Returns the next line following the sequence data (the next line
769 * that doesn't start with whitespace).
771 * @throws IOException
773 protected String parseSequence() throws IOException
775 StringBuilder sb = new StringBuilder(this.length);
776 String line = nextLine();
777 while (line != null && line.startsWith(" "))
780 String[] blocks = line.split(WHITESPACE);
783 * the first or last block on each line might be a position count - omit
785 for (int i = 0; i < blocks.length; i++)
789 Long.parseLong(blocks[i]);
790 // position counter - ignore it
791 } catch (NumberFormatException e)
793 // sequence data - append it
794 sb.append(blocks[i]);
799 this.sequenceString = sb.toString();
805 * Processes a feature line. If it declares a feature type of interest
806 * (currently, only CDS is processed), processes all of the associated lines
807 * (feature qualifiers), and returns the next line after that, otherwise
808 * simply returns the next line.
811 * the first line for the feature (with initial FT omitted for EMBL
814 * @throws IOException
816 protected String parseFeature(String line) throws IOException
818 String[] tokens = line.trim().split(WHITESPACE);
819 if (tokens.length < 2
820 || (!"CDS".equals(tokens[0]) && (!"source".equals(tokens[0]))))
824 if (tokens[0].equals("source"))
826 return parseSourceQualifiers(tokens);
828 return parseCDSFeature(tokens[1]);
833 * A data bean class to hold values parsed from one CDS Feature
837 String translation; // from /translation qualifier
839 String cdsLocation; // the raw value e.g. join(1..1234,2012..2837)
841 int codonStart = 1; // from /codon_start qualifier
843 String proteinName; // from /product qualifier; used for protein description
845 String proteinId; // from /protein_id qualifier
847 List<DBRefEntry> xrefs = new ArrayList<>(); // from /db_xref qualifiers
849 Map<String, String> cdsProps = new Hashtable<>(); // other qualifiers