src/jalview/io/EMBLLikeFlatFile.java

   1 package jalview.io;
   2
   3 import java.io.IOException;
   4 import java.text.ParseException;
   5 import java.util.ArrayList;
   6 import java.util.Arrays;
   7 import java.util.HashMap;
   8 import java.util.Hashtable;
   9 import java.util.List;
  10 import java.util.Locale;
  11 import java.util.Map;
  12 import java.util.Map.Entry;
  13 import java.util.TreeMap;
  14
  15 import jalview.bin.Cache;
  16 import jalview.datamodel.DBRefEntry;
  17 import jalview.datamodel.DBRefSource;
  18 import jalview.datamodel.FeatureProperties;
  19 import jalview.datamodel.Mapping;
  20 import jalview.datamodel.Sequence;
  21 import jalview.datamodel.SequenceFeature;
  22 import jalview.datamodel.SequenceI;
  23 import jalview.util.DBRefUtils;
  24 import jalview.util.DnaUtils;
  25 import jalview.util.MapList;
  26 import jalview.util.MappingUtils;
  27
  28 /**
  29  * A base class to support parsing of GenBank, EMBL or DDBJ flat file format
  30  * data. Example files (rather than formal specifications) are provided at
  31  *
  32  * <pre>
  33  * https://ena-docs.readthedocs.io/en/latest/submit/fileprep/flat-file-example.html
  34  * https://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html
  35  * </pre>
  36  *
  37  * or to compare the same entry, see
  38  *
  39  * <pre>
  40  * https://www.ebi.ac.uk/ena/browser/api/embl/X81322.1
  41  * https://www.ncbi.nlm.nih.gov/nuccore/X81322.1
  42  * </pre>
  43  *
  44  * The feature table part of the file has a common definition, only the start of
  45  * each line is formatted differently in GenBank and EMBL. See
  46  * http://www.insdc.org/files/feature_table.html#7.1.
  47  */
  48 public abstract class EMBLLikeFlatFile extends AlignFile
  49 {
  50   protected static final String LOCATION = "location";
  51
  52   protected static final String QUOTE = "\"";
  53
  54   protected static final String DOUBLED_QUOTE = QUOTE + QUOTE;
  55
  56   protected static final String WHITESPACE = "\\s+";
  57
  58   /**
  59    * Removes leading or trailing double quotes (") unless doubled, and changes
  60    * any 'escaped' (doubled) double quotes to single characters. As per the
  61    * Feature Table specification for Qualifiers, Free Text.
  62    *
  63    * @param value
  64    * @return
  65    */
  66   protected static String removeQuotes(String value)
  67   {
  68     if (value == null)
  69     {
  70       return null;
  71     }
  72     if (value.startsWith(QUOTE) && !value.startsWith(DOUBLED_QUOTE))
  73     {
  74       value = value.substring(1);
  75     }
  76     if (value.endsWith(QUOTE) && !value.endsWith(DOUBLED_QUOTE))
  77     {
  78       value = value.substring(0, value.length() - 1);
  79     }
  80     value = value.replace(DOUBLED_QUOTE, QUOTE);
  81     return value;
  82   }
  83
  84   /**
  85    * Truncates (if necessary) the exon intervals to match 3 times the length of
  86    * the protein; also accepts 3 bases longer (for stop codon not included in
  87    * protein)
  88    *
  89    * @param proteinLength
  90    * @param exon
  91    *          an array of [start, end, start, end...] intervals
  92    * @return the same array (if unchanged) or a truncated copy
  93    */
  94   protected static int[] adjustForProteinLength(int proteinLength,
  95           int[] exon)
  96   {
  97     if (proteinLength <= 0 || exon == null)
  98     {
  99       return exon;
 100     }
 101     int expectedCdsLength = proteinLength * 3;
 102     int exonLength = MappingUtils.getLength(Arrays.asList(exon));
 103
 104     /*
 105      * if exon length matches protein, or is shorter, or longer by the
 106      * length of a stop codon (3 bases), then leave it unchanged
 107      */
 108     if (expectedCdsLength >= exonLength
 109             || expectedCdsLength == exonLength - 3)
 110     {
 111       return exon;
 112     }
 113
 114     int origxon[];
 115     int sxpos = -1;
 116     int endxon = 0;
 117     origxon = new int[exon.length];
 118     System.arraycopy(exon, 0, origxon, 0, exon.length);
 119     int cdspos = 0;
 120     for (int x = 0; x < exon.length; x += 2)
 121     {
 122       cdspos += Math.abs(exon[x + 1] - exon[x]) + 1;
 123       if (expectedCdsLength <= cdspos)
 124       {
 125         // advanced beyond last codon.
 126         sxpos = x;
 127         if (expectedCdsLength != cdspos)
 128         {
 129           // System.err
 130           // .println("Truncating final exon interval on region by "
 131           // + (cdspos - cdslength));
 132         }
 133
 134         /*
 135          * shrink the final exon - reduce end position if forward
 136          * strand, increase it if reverse
 137          */
 138         if (exon[x + 1] >= exon[x])
 139         {
 140           endxon = exon[x + 1] - cdspos + expectedCdsLength;
 141         }
 142         else
 143         {
 144           endxon = exon[x + 1] + cdspos - expectedCdsLength;
 145         }
 146         break;
 147       }
 148     }
 149
 150     if (sxpos != -1)
 151     {
 152       // and trim the exon interval set if necessary
 153       int[] nxon = new int[sxpos + 2];
 154       System.arraycopy(exon, 0, nxon, 0, sxpos + 2);
 155       nxon[sxpos + 1] = endxon; // update the end boundary for the new exon
 156                                 // set
 157       exon = nxon;
 158     }
 159     return exon;
 160   }
 161
 162   /*
 163    * when true, interpret the mol_type 'source' feature attribute
 164    * and generate an RNA sequence from the DNA record
 165    */
 166   protected boolean produceRna=true;
 167
 168
 169   /*
 170    * values parsed from the data file
 171    */
 172   protected String sourceDb;
 173
 174   protected String accession;
 175
 176   protected String version;
 177
 178   protected String description;
 179
 180   protected int length = 128;
 181
 182   protected List<DBRefEntry> dbrefs;
 183
 184   protected boolean sequenceStringIsRNA=false;
 185
 186   protected String sequenceString;
 187
 188   protected Map<String, CdsData> cds;
 189
 190   /**
 191    * Constructor
 192    *
 193    * @param fp
 194    * @param sourceId
 195    * @throws IOException
 196    */
 197   public EMBLLikeFlatFile(FileParse fp, String sourceId) throws IOException
 198   {
 199     super(false, fp); // don't parse immediately
 200     this.sourceDb = sourceId;
 201     dbrefs = new ArrayList<>();
 202
 203     /*
 204      * using TreeMap gives CDS sequences in alphabetical, so readable, order
 205      */
 206     cds = new TreeMap<>(String.CASE_INSENSITIVE_ORDER);
 207
 208     parse();
 209   }
 210
 211   /**
 212    * process attributes for 'source' until the next FT feature entry
 213    * only interested in 'mol_type'
 214    * @param tokens
 215    * @return
 216    * @throws IOException
 217    */
 218   private String parseSourceQualifiers(String[] tokens) throws IOException
 219   {
 220     if (!"source".equals(tokens[0]))
 221     {
 222       throw (new RuntimeException("Not given a 'source' qualifier line"));
 223     }
 224     // search for mol_type attribute
 225
 226     StringBuilder sb = new StringBuilder().append(tokens[1]); // extent of
 227                                                               // sequence
 228
 229     String line = parseFeatureQualifier(sb, false);
 230     while (line != null)
 231     {
 232       if (!line.startsWith("FT    ")) // four spaces, end of this feature table
 233                                       // entry
 234       {
 235         return line;
 236       }
 237
 238       // case sensitive ?
 239       int p = line.indexOf("\\mol_type");
 240       int qs = line.indexOf("\"", p);
 241       int qe = line.indexOf("\"", qs + 1);
 242       String qualifier=line.substring(qs,qe).toLowerCase(Locale.ROOT);
 243       if (qualifier.indexOf("rna") > -1)
 244       {
 245         sequenceStringIsRNA = true;
 246       }
 247       if (qualifier.indexOf("dna") > -1)
 248       {
 249         sequenceStringIsRNA = false;
 250       }
 251       line=parseFeatureQualifier(sb, false);
 252     }
 253     return line;
 254   }
 255
 256
 257   /**
 258    * Parses one (GenBank or EMBL format) CDS feature, saves the parsed data, and
 259    * returns the next line
 260    *
 261    * @param location
 262    * @return
 263    * @throws IOException
 264    */
 265   protected String parseCDSFeature(String location) throws IOException
 266   {
 267     String line;
 268
 269     /*
 270      * parse location, which can be over >1 line e.g. EAW51554
 271      */
 272     CdsData data = new CdsData();
 273     StringBuilder sb = new StringBuilder().append(location);
 274     line = parseFeatureQualifier(sb, false);
 275     data.cdsLocation = sb.toString();
 276
 277     while (line != null)
 278     {
 279       if (!isFeatureContinuationLine(line))
 280       {
 281         // e.g. start of next feature "FT source..."
 282         break;
 283       }
 284
 285       /*
 286        * extract qualifier, e.g. FT    /protein_id="CAA37824.1"
 287        * - the value may extend over more than one line
 288        * - if the value has enclosing quotes, these are removed
 289        * - escaped double quotes ("") are reduced to a single character
 290        */
 291       int slashPos = line.indexOf('/');
 292       if (slashPos == -1)
 293       {
 294         Cache.error("Unexpected EMBL line ignored: " + line);
 295         line = nextLine();
 296         continue;
 297       }
 298       int eqPos = line.indexOf('=', slashPos + 1);
 299       if (eqPos == -1)
 300       {
 301         // can happen, e.g. /ribosomal_slippage
 302         line = nextLine();
 303         continue;
 304       }
 305       String qualifier = line.substring(slashPos + 1, eqPos);
 306       String value = line.substring(eqPos + 1);
 307       value = removeQuotes(value);
 308       sb = new StringBuilder().append(value);
 309       boolean asText = !"translation".equals(qualifier);
 310       line = parseFeatureQualifier(sb, asText);
 311       String featureValue = sb.toString();
 312
 313       if ("protein_id".equals(qualifier))
 314       {
 315         data.proteinId = featureValue;
 316       }
 317       else if ("codon_start".equals(qualifier))
 318       {
 319         try
 320         {
 321           data.codonStart = Integer.parseInt(featureValue.trim());
 322         } catch (NumberFormatException e)
 323         {
 324           Cache.error("Invalid codon_start in XML for " + this.accession
 325                   + ": " + e.getMessage());
 326         }
 327       }
 328       else if ("db_xref".equals(qualifier))
 329       {
 330         String[] parts = featureValue.split(":");
 331         if (parts.length == 2)
 332         {
 333           String db = parts[0].trim();
 334           db = DBRefUtils.getCanonicalName(db);
 335           DBRefEntry dbref = new DBRefEntry(db, "0", parts[1].trim());
 336           data.xrefs.add(dbref);
 337         }
 338       }
 339       else if ("product".equals(qualifier))
 340       {
 341         data.proteinName = featureValue;
 342       }
 343       else if ("translation".equals(qualifier))
 344       {
 345         data.translation = featureValue;
 346       }
 347       else if (!"".equals(featureValue))
 348       {
 349         // throw anything else into the additional properties hash
 350         data.cdsProps.put(qualifier, featureValue);
 351       }
 352     }
 353
 354     if (data.proteinId != null)
 355     {
 356       this.cds.put(data.proteinId, data);
 357     }
 358     else
 359     {
 360       Cache.error("Ignoring CDS feature with no protein_id for "
 361               + sourceDb + ":" + accession);
 362     }
 363
 364     return line;
 365   }
 366
 367   protected abstract boolean isFeatureContinuationLine(String line);
 368
 369   /**
 370    * Output (print) is not (yet) implemented for flat file format
 371    */
 372   @Override
 373   public String print(SequenceI[] seqs, boolean jvsuffix)
 374   {
 375     return null;
 376   }
 377
 378   /**
 379    * Constructs and saves the sequence from parsed components
 380    */
 381   protected void buildSequence()
 382   {
 383     if (this.accession == null || this.sequenceString == null)
 384     {
 385       Cache.error("Failed to parse data from EMBL");
 386       return;
 387     }
 388
 389     String name = this.accession;
 390     if (this.sourceDb != null)
 391     {
 392       name = this.sourceDb + "|" + name;
 393     }
 394
 395     if (produceRna && sequenceStringIsRNA)
 396     {
 397       sequenceString = sequenceString.replace('T', 'U').replace('t', 'u');
 398     }
 399
 400     SequenceI seq = new Sequence(name, this.sequenceString);
 401     seq.setDescription(this.description);
 402
 403     /*
 404      * add a DBRef to itself
 405      */
 406     DBRefEntry selfRef = new DBRefEntry(sourceDb, version, accession);
 407     int[] startEnd = new int[] { 1, seq.getLength() };
 408     selfRef.setMap(new Mapping(null, startEnd, startEnd, 1, 1));
 409     seq.addDBRef(selfRef);
 410
 411     for (DBRefEntry dbref : this.dbrefs)
 412     {
 413       seq.addDBRef(dbref);
 414     }
 415
 416     processCDSFeatures(seq);
 417
 418     seq.deriveSequence();
 419
 420     addSequence(seq);
 421   }
 422
 423   /**
 424    * Process the CDS features, including generation of cross-references and
 425    * mappings to the protein products (translation)
 426    *
 427    * @param seq
 428    */
 429   protected void processCDSFeatures(SequenceI seq)
 430   {
 431     /*
 432      * record protein products found to avoid duplication i.e. >1 CDS with
 433      * the same /protein_id [though not sure I can find an example of this]
 434      */
 435     Map<String, SequenceI> proteins = new HashMap<>();
 436     for (CdsData data : cds.values())
 437     {
 438       processCDSFeature(seq, data, proteins);
 439     }
 440   }
 441
 442   /**
 443    * Processes data for one parsed CDS feature to
 444    * <ul>
 445    * <li>create a protein product sequence for the translation</li>
 446    * <li>create a cross-reference to protein with mapping from dna</li>
 447    * <li>add a CDS feature to the sequence for each CDS start-end range</li>
 448    * <li>add any CDS dbrefs to the sequence and to the protein product</li>
 449    * </ul>
 450    *
 451    * @param SequenceI
 452    *          dna
 453    * @param proteins
 454    *          map of protein products so far derived from CDS data
 455    */
 456   void processCDSFeature(SequenceI dna, CdsData data,
 457           Map<String, SequenceI> proteins)
 458   {
 459     /*
 460      * parse location into a list of [start, end, start, end] positions
 461      */
 462     int[] exons = getCdsRanges(this.accession, data.cdsLocation);
 463
 464     MapList maplist = buildMappingToProtein(dna, exons, data);
 465
 466     int exonNumber = 0;
 467
 468     for (int xint = 0; exons != null && xint < exons.length - 1; xint += 2)
 469     {
 470       int exonStart = exons[xint];
 471       int exonEnd = exons[xint + 1];
 472       int begin = Math.min(exonStart, exonEnd);
 473       int end = Math.max(exonStart, exonEnd);
 474       exonNumber++;
 475       String desc = String.format("Exon %d for protein EMBLCDS:%s",
 476               exonNumber, data.proteinId);
 477
 478       SequenceFeature sf = new SequenceFeature("CDS", desc, begin, end,
 479               this.sourceDb);
 480       for (Entry<String, String> val : data.cdsProps.entrySet())
 481       {
 482         sf.setValue(val.getKey(), val.getValue());
 483       }
 484
 485       sf.setEnaLocation(data.cdsLocation);
 486       boolean forwardStrand = exonStart <= exonEnd;
 487       sf.setStrand(forwardStrand ? "+" : "-");
 488       sf.setPhase(String.valueOf(data.codonStart - 1));
 489       sf.setValue(FeatureProperties.EXONPOS, exonNumber);
 490       sf.setValue(FeatureProperties.EXONPRODUCT, data.proteinName);
 491
 492       dna.addSequenceFeature(sf);
 493     }
 494
 495     boolean hasUniprotDbref = false;
 496     for (DBRefEntry xref : data.xrefs)
 497     {
 498       dna.addDBRef(xref);
 499       if (xref.getSource().equals(DBRefSource.UNIPROT))
 500       {
 501         /*
 502          * construct (or find) the sequence for (data.protein_id, data.translation)
 503          */
 504         SequenceI protein = buildProteinProduct(dna, xref, data, proteins);
 505         Mapping map = new Mapping(protein, maplist);
 506         map.setMappedFromId(data.proteinId);
 507         xref.setMap(map);
 508
 509         /*
 510          * add DBRefs with mappings from dna to protein and the inverse
 511          */
 512         DBRefEntry db1 = new DBRefEntry(sourceDb, version, accession);
 513         db1.setMap(new Mapping(dna, maplist.getInverse()));
 514         protein.addDBRef(db1);
 515
 516         hasUniprotDbref = true;
 517       }
 518     }
 519
 520     /*
 521      * if we have a product (translation) but no explicit Uniprot dbref
 522      * (example: EMBL M19487 protein_id AAB02592.1)
 523      * then construct mappings to an assumed EMBLCDSPROTEIN accession
 524      */
 525     if (!hasUniprotDbref)
 526     {
 527       SequenceI protein = proteins.get(data.proteinId);
 528       if (protein == null)
 529       {
 530         protein = new Sequence(data.proteinId, data.translation);
 531         protein.setDescription(data.proteinName);
 532         proteins.put(data.proteinId, protein);
 533       }
 534       // assuming CDSPROTEIN sequence version = dna version (?!)
 535       DBRefEntry db1 = new DBRefEntry(DBRefSource.EMBLCDSProduct,
 536               this.version, data.proteinId);
 537       protein.addDBRef(db1);
 538
 539       DBRefEntry dnaToEmblProteinRef = new DBRefEntry(
 540               DBRefSource.EMBLCDSProduct, this.version, data.proteinId);
 541       Mapping map = new Mapping(protein, maplist);
 542       map.setMappedFromId(data.proteinId);
 543       dnaToEmblProteinRef.setMap(map);
 544       dna.addDBRef(dnaToEmblProteinRef);
 545     }
 546
 547     /*
 548      * comment brought forward from EmblXmlSource, lines 447-451:
 549      * TODO: if retrieved from EMBLCDS, add a DBRef back to the parent EMBL
 550      * sequence with the exon  map; if given a dataset reference, search
 551      * dataset for parent EMBL sequence if it exists and set its map;
 552      * make a new feature annotating the coding contig
 553      */
 554   }
 555
 556   /**
 557    * Computes a mapping from CDS positions in DNA sequence to protein product
 558    * positions, with allowance for stop codon or incomplete start codon
 559    *
 560    * @param dna
 561    * @param exons
 562    * @param data
 563    * @return
 564    */
 565   MapList buildMappingToProtein(final SequenceI dna, final int[] exons,
 566           final CdsData data)
 567   {
 568     MapList dnaToProteinMapping = null;
 569     int peptideLength = data.translation.length();
 570
 571     int[] proteinRange = new int[] { 1, peptideLength };
 572     if (exons != null && exons.length > 0)
 573     {
 574       /*
 575        * We were able to parse 'location'; do a final
 576        * product length truncation check
 577        */
 578       int[] cdsRanges = adjustForProteinLength(peptideLength, exons);
 579       dnaToProteinMapping = new MapList(cdsRanges, proteinRange, 3, 1);
 580     }
 581     else
 582     {
 583       /*
 584        * workaround until we handle all 'location' formats fully
 585        * e.g. X53828.1:60..1058 or <123..>289
 586        */
 587       Cache.error(String.format(
 588               "Implementation Notice: EMBLCDS location '%s'not properly supported yet"
 589                       + " - Making up the CDNA region of (%s:%s)... may be incorrect",
 590               data.cdsLocation, sourceDb, this.accession));
 591
 592       int completeCodonsLength = 1 - data.codonStart + dna.getLength();
 593       int mappedDnaEnd = dna.getEnd();
 594       if (peptideLength * 3 == completeCodonsLength)
 595       {
 596         // this might occur for CDS sequences where no features are marked
 597         Cache.warn("Assuming no stop codon at end of cDNA fragment");
 598         mappedDnaEnd = dna.getEnd();
 599       }
 600       else if ((peptideLength + 1) * 3 == completeCodonsLength)
 601       {
 602         Cache.warn("Assuming stop codon at end of cDNA fragment");
 603         mappedDnaEnd = dna.getEnd() - 3;
 604       }
 605
 606       if (mappedDnaEnd != -1)
 607       {
 608         int[] cdsRanges = new int[] {
 609             dna.getStart() + (data.codonStart - 1), mappedDnaEnd };
 610         dnaToProteinMapping = new MapList(cdsRanges, proteinRange, 3, 1);
 611       }
 612     }
 613
 614     return dnaToProteinMapping;
 615   }
 616
 617   /**
 618    * Constructs a sequence for the protein product for the CDS data (if there is
 619    * one), and dbrefs with mappings from CDS to protein and the reverse
 620    *
 621    * @param dna
 622    * @param xref
 623    * @param data
 624    * @param proteins
 625    * @return
 626    */
 627   SequenceI buildProteinProduct(SequenceI dna, DBRefEntry xref,
 628           CdsData data, Map<String, SequenceI> proteins)
 629   {
 630     /*
 631      * check we have some data to work with
 632      */
 633     if (data.proteinId == null || data.translation == null)
 634     {
 635       return null;
 636     }
 637
 638     /*
 639      * Construct the protein sequence (if not already seen)
 640      */
 641     String proteinSeqName = xref.getSource() + "|" + xref.getAccessionId();
 642     SequenceI protein = proteins.get(proteinSeqName);
 643     if (protein == null)
 644     {
 645       protein = new Sequence(proteinSeqName, data.translation, 1,
 646               data.translation.length());
 647       protein.setDescription(data.proteinName != null ? data.proteinName
 648               : "Protein Product from " + sourceDb);
 649       proteins.put(proteinSeqName, protein);
 650     }
 651
 652     return protein;
 653   }
 654
 655   /**
 656    * Returns the CDS location as a single array of [start, end, start, end...]
 657    * positions. If on the reverse strand, these will be in descending order.
 658    *
 659    * @param accession
 660    * @param location
 661    * @return
 662    */
 663   protected int[] getCdsRanges(String accession, String location)
 664   {
 665     if (location == null)
 666     {
 667       return new int[] {};
 668     }
 669
 670     try
 671     {
 672       List<int[]> ranges = DnaUtils.parseLocation(location);
 673       return MappingUtils.rangeListToArray(ranges);
 674     } catch (ParseException e)
 675     {
 676       Cache.warn(
 677               String.format("Not parsing inexact CDS location %s in ENA %s",
 678                       location, accession));
 679       return new int[] {};
 680     }
 681   }
 682
 683   /**
 684    * Reads the value of a feature (FT) qualifier from one or more lines of the
 685    * file, and returns the next line after that. Values are appended to the
 686    * string buffer, which should be already primed with the value read from the
 687    * first line for the qualifier (with any leading double quote removed).
 688    * Enclosing double quotes are removed, and escaped (repeated) double quotes
 689    * reduced to one only. For example for
 690    *
 691    * <pre>
 692    * FT      /note="gene_id=hCG28070.3
 693    * FT      ""foobar"" isoform=CRA_b"
 694    * the returned value is
 695    * gene_id=hCG28070.3 "foobar" isoform=CRA_b
 696    * </pre>
 697    *
 698    * Note the side-effect of this method, to advance data reading to the next
 699    * line after the feature qualifier (which could be another qualifier, a
 700    * different feature, a non-feature line, or null at end of file).
 701    *
 702    * @param sb
 703    *          a string buffer primed with the first line of the value
 704    * @param asText
 705    * @return
 706    * @throws IOException
 707    */
 708   String parseFeatureQualifier(StringBuilder sb, boolean asText)
 709           throws IOException
 710   {
 711     String line;
 712     while ((line = nextLine()) != null)
 713     {
 714       if (!isFeatureContinuationLine(line))
 715       {
 716         break; // reached next feature or other input line
 717       }
 718       String[] tokens = line.split(WHITESPACE);
 719       if (tokens.length < 2)
 720       {
 721         Cache.error("Ignoring bad EMBL line for " + this.accession
 722                 + ": " + line);
 723         break;
 724       }
 725       if (tokens[1].startsWith("/"))
 726       {
 727         break; // next feature qualifier
 728       }
 729
 730       /*
 731        * if text (e.g. /product), add a word separator for a new line,
 732        * else (e.g. /translation) don't
 733        */
 734       if (asText)
 735       {
 736         sb.append(" ");
 737       }
 738
 739       /*
 740        * remove trailing " and unescape doubled ""
 741        */
 742       String data = removeQuotes(tokens[1]);
 743       sb.append(data);
 744     }
 745
 746     return line;
 747   }
 748
 749   /**
 750    * Reads and saves the sequence, read from the lines following the ORIGIN
 751    * (GenBank) or SQ (EMBL) line. Whitespace and position counters are
 752    * discarded. Returns the next line following the sequence data (the next line
 753    * that doesn't start with whitespace).
 754    *
 755    * @throws IOException
 756    */
 757   protected String parseSequence() throws IOException
 758   {
 759     StringBuilder sb = new StringBuilder(this.length);
 760     String line = nextLine();
 761     while (line != null && line.startsWith(" "))
 762     {
 763       line = line.trim();
 764       String[] blocks = line.split(WHITESPACE);
 765
 766       /*
 767        * the first or last block on each line might be a position count - omit
 768        */
 769       for (int i = 0; i < blocks.length; i++)
 770       {
 771         try
 772         {
 773           Long.parseLong(blocks[i]);
 774           // position counter - ignore it
 775         } catch (NumberFormatException e)
 776         {
 777           // sequence data - append it
 778           sb.append(blocks[i]);
 779         }
 780       }
 781       line = nextLine();
 782     }
 783     this.sequenceString = sb.toString();
 784
 785     return line;
 786   }
 787
 788   /**
 789    * Processes a feature line. If it declares a feature type of interest
 790    * (currently, only CDS is processed), processes all of the associated lines
 791    * (feature qualifiers), and returns the next line after that, otherwise
 792    * simply returns the next line.
 793    *
 794    * @param line
 795    *          the first line for the feature (with initial FT omitted for EMBL
 796    *          format)
 797    * @return
 798    * @throws IOException
 799    */
 800   protected String parseFeature(String line) throws IOException
 801   {
 802     String[] tokens = line.trim().split(WHITESPACE);
 803     if (tokens.length < 2 || (!"CDS".equals(tokens[0]) && (!"source".equals(tokens[0]))))
 804     {
 805       return nextLine();
 806     }
 807     if (tokens[0].equals("source"))
 808     {
 809       return parseSourceQualifiers(tokens);
 810     }
 811     return parseCDSFeature(tokens[1]);
 812   }
 813 }
 814
 815 /**
 816  * A data bean class to hold values parsed from one CDS Feature
 817  */
 818 class CdsData
 819 {
 820   String translation; // from /translation qualifier
 821
 822   String cdsLocation; // the raw value e.g. join(1..1234,2012..2837)
 823
 824   int codonStart = 1; // from /codon_start qualifier
 825
 826   String proteinName; // from /product qualifier; used for protein description
 827
 828   String proteinId; // from /protein_id qualifier
 829
 830   List<DBRefEntry> xrefs = new ArrayList<>(); // from /db_xref qualifiers
 831
 832   Map<String, String> cdsProps = new Hashtable<>(); // other qualifiers
 833 }