src/jalview/ws/dbsources/EmblXmlSource.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   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.
  11  *
  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.
  16  *
  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.
  20  */
  21 package jalview.ws.dbsources;
  22
  23 import jalview.analysis.SequenceIdMatcher;
  24 import jalview.bin.Cache;
  25 import jalview.datamodel.Alignment;
  26 import jalview.datamodel.AlignmentI;
  27 import jalview.datamodel.DBRefEntry;
  28 import jalview.datamodel.DBRefSource;
  29 import jalview.datamodel.FeatureProperties;
  30 import jalview.datamodel.Mapping;
  31 import jalview.datamodel.Sequence;
  32 import jalview.datamodel.SequenceFeature;
  33 import jalview.datamodel.SequenceI;
  34 import jalview.util.DBRefUtils;
  35 import jalview.util.DnaUtils;
  36 import jalview.util.MapList;
  37 import jalview.util.MappingUtils;
  38 import jalview.util.MessageManager;
  39 import jalview.ws.ebi.EBIFetchClient;
  40 import jalview.xml.binding.embl.EntryType;
  41 import jalview.xml.binding.embl.EntryType.Feature;
  42 import jalview.xml.binding.embl.EntryType.Feature.Qualifier;
  43 import jalview.xml.binding.embl.ROOT;
  44 import jalview.xml.binding.embl.XrefType;
  45
  46 import java.io.File;
  47 import java.io.FileInputStream;
  48 import java.io.InputStream;
  49 import java.text.ParseException;
  50 import java.util.ArrayList;
  51 import java.util.Arrays;
  52 import java.util.Hashtable;
  53 import java.util.List;
  54 import java.util.Map;
  55 import java.util.Map.Entry;
  56 import java.util.regex.Pattern;
  57
  58 import javax.xml.bind.JAXBContext;
  59 import javax.xml.bind.JAXBElement;
  60 import javax.xml.bind.JAXBException;
  61 import javax.xml.stream.FactoryConfigurationError;
  62 import javax.xml.stream.XMLInputFactory;
  63 import javax.xml.stream.XMLStreamException;
  64 import javax.xml.stream.XMLStreamReader;
  65
  66 public abstract class EmblXmlSource extends EbiFileRetrievedProxy
  67 {
  68   /*
  69    * JAL-1856 Embl returns this text for query not found
  70    */
  71   private static final String EMBL_NOT_FOUND_REPLY = "ERROR 12 No entries found.";
  72
  73   private static final Pattern SPACE_PATTERN = Pattern.compile(" ");
  74
  75   public EmblXmlSource()
  76   {
  77     super();
  78   }
  79
  80   /**
  81    * Retrieves and parses an emblxml file, and returns an alignment containing
  82    * the parsed sequences, or null if none were found
  83    *
  84    * @param emprefx
  85    *          "EMBL" or "EMBLCDS" - anything else will not retrieve emblxml
  86    * @param query
  87    * @return
  88    * @throws Exception
  89    */
  90   protected AlignmentI getEmblSequenceRecords(String emprefx, String query)
  91           throws Exception
  92   {
  93     startQuery();
  94     EBIFetchClient dbFetch = new EBIFetchClient();
  95     File reply;
  96     try
  97     {
  98       reply = dbFetch.fetchDataAsFile(
  99               emprefx.toLowerCase() + ":" + query.trim(), "display=xml",
 100               "xml");
 101     } catch (Exception e)
 102     {
 103       stopQuery();
 104       throw new Exception(MessageManager.formatMessage(
 105               "exception.ebiembl_retrieval_failed_on", new String[]
 106               { emprefx.toLowerCase(), query.trim() }), e);
 107     }
 108     return getEmblSequenceRecords(emprefx, query, reply);
 109   }
 110
 111   /**
 112    * parse an emblxml file stored locally
 113    *
 114    * @param emprefx
 115    *          either EMBL or EMBLCDS strings are allowed - anything else will
 116    *          not retrieve emblxml
 117    * @param query
 118    * @param file
 119    *          the EMBL XML file containing the results of a query
 120    * @return
 121    * @throws Exception
 122    */
 123   protected AlignmentI getEmblSequenceRecords(String emprefx, String query,
 124           File reply) throws Exception
 125   {
 126     List<EntryType> entries = null;
 127     if (reply != null && reply.exists())
 128     {
 129       file = reply.getAbsolutePath();
 130       if (reply.length() > EMBL_NOT_FOUND_REPLY.length())
 131       {
 132         InputStream is = new FileInputStream(reply);
 133         entries = getEmblEntries(is);
 134       }
 135     }
 136
 137     /*
 138      * invalid accession gets a reply with no <entry> elements, text content of
 139      * EmbFile reads something like (e.g.) this ungrammatical phrase
 140      * Entry: <acc> display type is either not supported or entry is not found.
 141      */
 142     AlignmentI al = null;
 143     List<SequenceI> seqs = new ArrayList<>();
 144     List<SequenceI> peptides = new ArrayList<>();
 145     if (entries != null)
 146     {
 147       for (EntryType entry : entries)
 148       {
 149         SequenceI seq = getSequence(emprefx, entry, peptides);
 150         if (seq != null)
 151         {
 152           seqs.add(seq.deriveSequence());
 153           // place DBReferences on dataset and refer
 154         }
 155       }
 156       if (!seqs.isEmpty())
 157       {
 158         al = new Alignment(seqs.toArray(new SequenceI[seqs.size()]));
 159       }
 160       else
 161       {
 162         System.out.println(
 163                 "No record found for '" + emprefx + ":" + query + "'");
 164       }
 165     }
 166
 167     stopQuery();
 168     return al;
 169   }
 170
 171   /**
 172    * Reads the XML reply from file and unmarshals it to Java objects. Answers a
 173    * (possibly empty) list of <code>EntryType</code> objects.
 174    *
 175    * is
 176    *
 177    * @return
 178    */
 179   List<EntryType> getEmblEntries(InputStream is)
 180   {
 181     List<EntryType> entries = new ArrayList<>();
 182     try
 183     {
 184       JAXBContext jc = JAXBContext.newInstance("jalview.xml.binding.embl");
 185       XMLStreamReader streamReader = XMLInputFactory.newInstance()
 186               .createXMLStreamReader(is);
 187       javax.xml.bind.Unmarshaller um = jc.createUnmarshaller();
 188       JAXBElement<ROOT> rootElement =  um.unmarshal(streamReader, ROOT.class);
 189       ROOT root = rootElement.getValue();
 190
 191       /*
 192        * document root contains either "entry" or "entrySet"
 193        */
 194       if (root == null)
 195       {
 196         return entries;
 197       }
 198       if (root.getEntrySet() != null)
 199       {
 200         entries = root.getEntrySet().getEntry();
 201       }
 202       else if (root.getEntry() != null)
 203       {
 204         entries.add(root.getEntry());
 205       }
 206     } catch (JAXBException | XMLStreamException
 207             | FactoryConfigurationError e)
 208     {
 209       e.printStackTrace();
 210     }
 211     return entries;
 212   }
 213
 214   /**
 215    * A helper method to parse XML data and construct a sequence, with any
 216    * available database references and features
 217    *
 218    * @param emprefx
 219    * @param entry
 220    * @param peptides
 221    * @return
 222    */
 223   SequenceI getSequence(String sourceDb, EntryType entry,
 224           List<SequenceI> peptides)
 225   {
 226     String seqString = entry.getSequence();
 227     if (seqString == null)
 228     {
 229       return null;
 230     }
 231     seqString = seqString.replace(" ", "").replace("\n", "").replace("\t",
 232             "");
 233     String accession = entry.getAccession();
 234     SequenceI dna = new Sequence(sourceDb + "|" + accession, seqString);
 235
 236     dna.setDescription(entry.getDescription());
 237     String sequenceVersion = String.valueOf(entry.getVersion().intValue());
 238     DBRefEntry selfRref = new DBRefEntry(sourceDb, sequenceVersion,
 239             accession);
 240     dna.addDBRef(selfRref);
 241     selfRref.setMap(
 242             new Mapping(null, new int[]
 243             { 1, dna.getLength() }, new int[] { 1, dna.getLength() }, 1,
 244                     1));
 245
 246     /*
 247      * add db references
 248      */
 249     List<XrefType> xrefs = entry.getXref();
 250     if (xrefs != null)
 251     {
 252       for (XrefType xref : xrefs)
 253       {
 254         String acc = xref.getId();
 255         String source = DBRefUtils.getCanonicalName(xref.getDb());
 256         String version = xref.getSecondaryId();
 257         if (version == null || "".equals(version))
 258         {
 259           version = "0";
 260         }
 261         dna.addDBRef(new DBRefEntry(source, version, acc));
 262       }
 263     }
 264
 265     SequenceIdMatcher matcher = new SequenceIdMatcher(peptides);
 266     try
 267     {
 268       List<Feature> features = entry.getFeature();
 269       if (features != null)
 270       {
 271         for (Feature feature : features)
 272         {
 273           if (FeatureProperties.isCodingFeature(sourceDb,
 274                   feature.getName()))
 275           {
 276             parseCodingFeature(entry, feature, sourceDb, dna, peptides,
 277                     matcher);
 278           }
 279         }
 280       }
 281     } catch (Exception e)
 282     {
 283       System.err.println("EMBL Record Features parsing error!");
 284       System.err
 285               .println("Please report the following to help@jalview.org :");
 286       System.err.println("EMBL Record " + accession);
 287       System.err.println("Resulted in exception: " + e.getMessage());
 288       e.printStackTrace(System.err);
 289     }
 290
 291     return dna;
 292   }
 293
 294   /**
 295    * Extracts coding region and product from a CDS feature and decorates it with
 296    * annotations
 297    *
 298    * @param entry
 299    * @param feature
 300    * @param sourceDb
 301    * @param dna
 302    * @param peptides
 303    * @param matcher
 304    */
 305   void parseCodingFeature(EntryType entry, Feature feature, String sourceDb,
 306           SequenceI dna, List<SequenceI> peptides,
 307           SequenceIdMatcher matcher)
 308   {
 309     final boolean isEmblCdna = sourceDb.equals(DBRefSource.EMBLCDS);
 310     final String accession = entry.getAccession();
 311     final String sequenceVersion = entry.getVersion().toString();
 312
 313     int[] exons = getCdsRanges(entry.getAccession(), feature);
 314
 315     String translation = null;
 316     String proteinName = "";
 317     String proteinId = null;
 318     Map<String, String> vals = new Hashtable<>();
 319
 320     /*
 321      * codon_start 1/2/3 in EMBL corresponds to phase 0/1/2 in CDS
 322      * (phase is required for CDS features in GFF3 format)
 323      */
 324     int codonStart = 1;
 325
 326     /*
 327      * parse qualifiers, saving protein translation, protein id,
 328      * codon start position, product (name), and 'other values'
 329      */
 330     if (feature.getQualifier() != null)
 331     {
 332       for (Qualifier q : feature.getQualifier())
 333       {
 334         String qname = q.getName();
 335         String value = q.getValue();
 336         value = value == null ? ""
 337                 : value.trim().replace(" ", "").replace("\n", "")
 338                         .replace("\t", "");
 339         if (qname.equals("translation"))
 340         {
 341           translation = value;
 342         }
 343         else if (qname.equals("protein_id"))
 344         {
 345           proteinId = value;
 346         }
 347         else if (qname.equals("codon_start"))
 348         {
 349           try
 350           {
 351             codonStart = Integer.parseInt(value.trim());
 352           } catch (NumberFormatException e)
 353           {
 354             System.err.println("Invalid codon_start in XML for "
 355                     + entry.getAccession() + ": " + e.getMessage());
 356           }
 357         }
 358         else if (qname.equals("product"))
 359         {
 360           // sometimes name is returned e.g. for V00488
 361           proteinName = value;
 362         }
 363         else
 364         {
 365           // throw anything else into the additional properties hash
 366           if (!"".equals(value))
 367           {
 368             vals.put(qname, value);
 369           }
 370         }
 371       }
 372     }
 373
 374     DBRefEntry proteinToEmblProteinRef = null;
 375     exons = MappingUtils.removeStartPositions(codonStart - 1, exons);
 376
 377     SequenceI product = null;
 378     Mapping dnaToProteinMapping = null;
 379     if (translation != null && proteinName != null && proteinId != null)
 380     {
 381       int translationLength = translation.length();
 382
 383       /*
 384        * look for product in peptides list, if not found, add it
 385        */
 386       product = matcher.findIdMatch(proteinId);
 387       if (product == null)
 388       {
 389         product = new Sequence(proteinId, translation, 1,
 390                 translationLength);
 391         product.setDescription(((proteinName.length() == 0)
 392                 ? "Protein Product from " + sourceDb
 393                 : proteinName));
 394         peptides.add(product);
 395         matcher.add(product);
 396       }
 397
 398       // we have everything - create the mapping and perhaps the protein
 399       // sequence
 400       if (exons == null || exons.length == 0)
 401       {
 402         /*
 403          * workaround until we handle dna location for CDS sequence
 404          * e.g. location="X53828.1:60..1058" correctly
 405          */
 406         System.err.println(
 407                 "Implementation Notice: EMBLCDS records not properly supported yet - Making up the CDNA region of this sequence... may be incorrect ("
 408                         + sourceDb + ":" + entry.getAccession() + ")");
 409         int dnaLength = dna.getLength();
 410         if (translationLength * 3 == (1 - codonStart + dnaLength))
 411         {
 412           System.err.println(
 413                   "Not allowing for additional stop codon at end of cDNA fragment... !");
 414           // this might occur for CDS sequences where no features are marked
 415           exons = new int[] { dna.getStart() + (codonStart - 1),
 416               dna.getEnd() };
 417           dnaToProteinMapping = new Mapping(product, exons,
 418                   new int[]
 419                   { 1, translationLength }, 3, 1);
 420         }
 421         if ((translationLength + 1) * 3 == (1 - codonStart + dnaLength))
 422         {
 423           System.err.println(
 424                   "Allowing for additional stop codon at end of cDNA fragment... will probably cause an error in VAMSAs!");
 425           exons = new int[] { dna.getStart() + (codonStart - 1),
 426               dna.getEnd() - 3 };
 427           dnaToProteinMapping = new Mapping(product, exons,
 428                   new int[]
 429                   { 1, translationLength }, 3, 1);
 430         }
 431       }
 432       else
 433       {
 434         // Trim the exon mapping if necessary - the given product may only be a
 435         // fragment of a larger protein. (EMBL:AY043181 is an example)
 436
 437         if (isEmblCdna)
 438         {
 439           // TODO: Add a DbRef back to the parent EMBL sequence with the exon
 440           // map
 441           // if given a dataset reference, search dataset for parent EMBL
 442           // sequence if it exists and set its map
 443           // make a new feature annotating the coding contig
 444         }
 445         else
 446         {
 447           // final product length truncation check
 448           int[] cdsRanges = adjustForProteinLength(translationLength,
 449                   exons);
 450           dnaToProteinMapping = new Mapping(product, cdsRanges,
 451                   new int[]
 452                   { 1, translationLength }, 3, 1);
 453           if (product != null)
 454           {
 455             /*
 456              * make xref with mapping from protein to EMBL dna
 457              */
 458             DBRefEntry proteinToEmblRef = new DBRefEntry(DBRefSource.EMBL,
 459                     sequenceVersion, proteinId,
 460                     new Mapping(dnaToProteinMapping.getMap().getInverse()));
 461             product.addDBRef(proteinToEmblRef);
 462
 463             /*
 464              * make xref from protein to EMBLCDS; we assume here that the
 465              * CDS sequence version is same as dna sequence (?!)
 466              */
 467             MapList proteinToCdsMapList = new MapList(
 468                     new int[]
 469                     { 1, translationLength },
 470                     new int[]
 471                     { 1 + (codonStart - 1),
 472                         (codonStart - 1) + 3 * translationLength },
 473                     1, 3);
 474             DBRefEntry proteinToEmblCdsRef = new DBRefEntry(
 475                     DBRefSource.EMBLCDS, sequenceVersion, proteinId,
 476                     new Mapping(proteinToCdsMapList));
 477             product.addDBRef(proteinToEmblCdsRef);
 478
 479             /*
 480              * make 'direct' xref from protein to EMBLCDSPROTEIN
 481              */
 482             proteinToEmblProteinRef = new DBRefEntry(proteinToEmblCdsRef);
 483             proteinToEmblProteinRef.setSource(DBRefSource.EMBLCDSProduct);
 484             proteinToEmblProteinRef.setMap(null);
 485             product.addDBRef(proteinToEmblProteinRef);
 486           }
 487         }
 488       }
 489
 490       /*
 491        * add cds features to dna sequence
 492        */
 493       String cds = feature.getName(); // "CDS"
 494       for (int xint = 0; exons != null
 495               && xint < exons.length - 1; xint += 2)
 496       {
 497         int exonStart = exons[xint];
 498         int exonEnd = exons[xint + 1];
 499         int begin = Math.min(exonStart, exonEnd);
 500         int end = Math.max(exonStart, exonEnd);
 501         int exonNumber = xint / 2 + 1;
 502         String desc = String.format("Exon %d for protein '%s' EMBLCDS:%s",
 503                 exonNumber, proteinName, proteinId);
 504
 505         SequenceFeature sf = makeCdsFeature(cds, desc, begin, end, sourceDb,
 506                 vals);
 507
 508         sf.setEnaLocation(feature.getLocation());
 509         boolean forwardStrand = exonStart <= exonEnd;
 510         sf.setStrand(forwardStrand ? "+" : "-");
 511         sf.setPhase(String.valueOf(codonStart - 1));
 512         sf.setValue(FeatureProperties.EXONPOS, exonNumber);
 513         sf.setValue(FeatureProperties.EXONPRODUCT, proteinName);
 514
 515         dna.addSequenceFeature(sf);
 516       }
 517     }
 518
 519     /*
 520      * add feature dbRefs to sequence, and mappings for Uniprot xrefs
 521      */
 522     boolean hasUniprotDbref = false;
 523     List<XrefType> xrefs = feature.getXref();
 524     if (xrefs != null)
 525     {
 526       boolean mappingUsed = false;
 527       for (XrefType xref : xrefs)
 528       {
 529         /*
 530          * ensure UniProtKB/Swiss-Prot converted to UNIPROT
 531          */
 532         String source = DBRefUtils.getCanonicalName(xref.getDb());
 533         String version = xref.getSecondaryId();
 534         if (version == null || "".equals(version))
 535         {
 536           version = "0";
 537         }
 538         DBRefEntry dbref = new DBRefEntry(source, version, xref.getId());
 539         DBRefEntry proteinDbRef = new DBRefEntry(source, version,
 540                 dbref.getAccessionId());
 541         if (source.equals(DBRefSource.UNIPROT))
 542         {
 543           String proteinSeqName = DBRefSource.UNIPROT + "|"
 544                   + dbref.getAccessionId();
 545           if (dnaToProteinMapping != null
 546                   && dnaToProteinMapping.getTo() != null)
 547           {
 548             if (mappingUsed)
 549             {
 550               /*
 551                * two or more Uniprot xrefs for the same CDS -
 552                * each needs a distinct Mapping (as to a different sequence)
 553                */
 554               dnaToProteinMapping = new Mapping(dnaToProteinMapping);
 555             }
 556             mappingUsed = true;
 557
 558             /*
 559              * try to locate the protein mapped to (possibly by a
 560              * previous CDS feature); if not found, construct it from
 561              * the EMBL translation
 562              */
 563             SequenceI proteinSeq = matcher.findIdMatch(proteinSeqName);
 564             if (proteinSeq == null)
 565             {
 566               proteinSeq = new Sequence(proteinSeqName,
 567                       product.getSequenceAsString());
 568               matcher.add(proteinSeq);
 569               peptides.add(proteinSeq);
 570             }
 571             dnaToProteinMapping.setTo(proteinSeq);
 572             dnaToProteinMapping.setMappedFromId(proteinId);
 573             proteinSeq.addDBRef(proteinDbRef);
 574             dbref.setMap(dnaToProteinMapping);
 575           }
 576           hasUniprotDbref = true;
 577         }
 578         if (product != null)
 579         {
 580           /*
 581            * copy feature dbref to our protein product
 582            */
 583           DBRefEntry pref = proteinDbRef;
 584           pref.setMap(null); // reference is direct
 585           product.addDBRef(pref);
 586           // Add converse mapping reference
 587           if (dnaToProteinMapping != null)
 588           {
 589             Mapping pmap = new Mapping(dna,
 590                     dnaToProteinMapping.getMap().getInverse());
 591             pref = new DBRefEntry(sourceDb, sequenceVersion, accession);
 592             pref.setMap(pmap);
 593             if (dnaToProteinMapping.getTo() != null)
 594             {
 595               dnaToProteinMapping.getTo().addDBRef(pref);
 596             }
 597           }
 598         }
 599         dna.addDBRef(dbref);
 600       }
 601     }
 602
 603     /*
 604      * if we have a product (translation) but no explicit Uniprot dbref
 605      * (example: EMBL AAFI02000057 protein_id EAL65544.1)
 606      * then construct mappings to an assumed EMBLCDSPROTEIN accession
 607      */
 608     if (!hasUniprotDbref && product != null)
 609     {
 610       if (proteinToEmblProteinRef == null)
 611       {
 612         // assuming CDSPROTEIN sequence version = dna version (?!)
 613         proteinToEmblProteinRef = new DBRefEntry(DBRefSource.EMBLCDSProduct,
 614                 sequenceVersion, proteinId);
 615       }
 616       product.addDBRef(proteinToEmblProteinRef);
 617
 618       if (dnaToProteinMapping != null
 619               && dnaToProteinMapping.getTo() != null)
 620       {
 621         DBRefEntry dnaToEmblProteinRef = new DBRefEntry(
 622                 DBRefSource.EMBLCDSProduct, sequenceVersion,
 623                 proteinId);
 624         dnaToEmblProteinRef.setMap(dnaToProteinMapping);
 625         dnaToProteinMapping.setMappedFromId(proteinId);
 626         dna.addDBRef(dnaToEmblProteinRef);
 627       }
 628     }
 629   }
 630
 631   @Override
 632   public boolean isDnaCoding()
 633   {
 634     return true;
 635   }
 636
 637   /**
 638    * Returns the CDS positions as a single array of [start, end, start, end...]
 639    * positions. If on the reverse strand, these will be in descending order.
 640    *
 641    * @param accession
 642    * @param feature
 643    * @return
 644    */
 645   protected int[] getCdsRanges(String accession, Feature feature)
 646   {
 647     String location = feature.getLocation();
 648     if (location == null)
 649     {
 650       return new int[] {};
 651     }
 652
 653     try
 654     {
 655       // BH 2019.10.06 note: exo-location of "X53828.1:60..1058" causes
 656       // RemoteFormatTest to fail at this point.
 657       List<int[]> ranges = DnaUtils.parseLocation(location);
 658       return listToArray(ranges);
 659     } catch (ParseException e)
 660     {
 661       Cache.log.warn(
 662               String.format("Not parsing inexact CDS location %s in ENA %s",
 663                       location, accession));
 664       return new int[] {};
 665     }
 666   }
 667
 668   /**
 669    * Converts a list of [start, end] ranges to a single array of [start, end,
 670    * start, end ...]
 671    *
 672    * @param ranges
 673    * @return
 674    */
 675   int[] listToArray(List<int[]> ranges)
 676   {
 677     int[] result = new int[ranges.size() * 2];
 678     int i = 0;
 679     for (int[] range : ranges)
 680     {
 681       result[i++] = range[0];
 682       result[i++] = range[1];
 683     }
 684     return result;
 685   }
 686
 687   /**
 688    * Helper method to construct a SequenceFeature for one cds range
 689    *
 690    * @param type
 691    *          feature type ("CDS")
 692    * @param desc
 693    *          description
 694    * @param begin
 695    *          start position
 696    * @param end
 697    *          end position
 698    * @param group
 699    *          feature group
 700    * @param vals
 701    *          map of 'miscellaneous values' for feature
 702    * @return
 703    */
 704   protected SequenceFeature makeCdsFeature(String type, String desc,
 705           int begin, int end, String group, Map<String, String> vals)
 706   {
 707     SequenceFeature sf = new SequenceFeature(type, desc, begin, end, group);
 708     if (!vals.isEmpty())
 709     {
 710       StringBuilder sb = new StringBuilder();
 711       boolean first = true;
 712       for (Entry<String, String> val : vals.entrySet())
 713       {
 714         if (!first)
 715         {
 716           sb.append(";");
 717         }
 718         sb.append(val.getKey()).append("=").append(val.getValue());
 719         first = false;
 720         sf.setValue(val.getKey(), val.getValue());
 721       }
 722       sf.setAttributes(sb.toString());
 723     }
 724     return sf;
 725   }
 726
 727   /**
 728    * Truncates (if necessary) the exon intervals to match 3 times the length of
 729    * the protein; also accepts 3 bases longer (for stop codon not included in
 730    * protein)
 731    *
 732    * @param proteinLength
 733    * @param exon
 734    *          an array of [start, end, start, end...] intervals
 735    * @return the same array (if unchanged) or a truncated copy
 736    */
 737   static int[] adjustForProteinLength(int proteinLength, int[] exon)
 738   {
 739     if (proteinLength <= 0 || exon == null)
 740     {
 741       return exon;
 742     }
 743     int expectedCdsLength = proteinLength * 3;
 744     int exonLength = MappingUtils.getLength(Arrays.asList(exon));
 745
 746     /*
 747      * if exon length matches protein, or is shorter, or longer by the
 748      * length of a stop codon (3 bases), then leave it unchanged
 749      */
 750     if (expectedCdsLength >= exonLength
 751             || expectedCdsLength == exonLength - 3)
 752     {
 753       return exon;
 754     }
 755
 756     int origxon[];
 757     int sxpos = -1;
 758     int endxon = 0;
 759     origxon = new int[exon.length];
 760     System.arraycopy(exon, 0, origxon, 0, exon.length);
 761     int cdspos = 0;
 762     for (int x = 0; x < exon.length; x += 2)
 763     {
 764       cdspos += Math.abs(exon[x + 1] - exon[x]) + 1;
 765       if (expectedCdsLength <= cdspos)
 766       {
 767         // advanced beyond last codon.
 768         sxpos = x;
 769         if (expectedCdsLength != cdspos)
 770         {
 771           // System.err
 772           // .println("Truncating final exon interval on region by "
 773           // + (cdspos - cdslength));
 774         }
 775
 776         /*
 777          * shrink the final exon - reduce end position if forward
 778          * strand, increase it if reverse
 779          */
 780         if (exon[x + 1] >= exon[x])
 781         {
 782           endxon = exon[x + 1] - cdspos + expectedCdsLength;
 783         }
 784         else
 785         {
 786           endxon = exon[x + 1] + cdspos - expectedCdsLength;
 787         }
 788         break;
 789       }
 790     }
 791
 792     if (sxpos != -1)
 793     {
 794       // and trim the exon interval set if necessary
 795       int[] nxon = new int[sxpos + 2];
 796       System.arraycopy(exon, 0, nxon, 0, sxpos + 2);
 797       nxon[sxpos + 1] = endxon; // update the end boundary for the new exon
 798                                 // set
 799       exon = nxon;
 800     }
 801     return exon;
 802   }
 803
 804 }