src/jalview/io/vcf/VCFLoader.java

   1 package jalview.io.vcf;
   2
   3 import htsjdk.samtools.util.CloseableIterator;
   4 import htsjdk.variant.variantcontext.Allele;
   5 import htsjdk.variant.variantcontext.VariantContext;
   6 import htsjdk.variant.vcf.VCFHeader;
   7 import htsjdk.variant.vcf.VCFHeaderLine;
   8 import htsjdk.variant.vcf.VCFHeaderLineCount;
   9
  10 import jalview.analysis.AlignmentUtils;
  11 import jalview.analysis.Dna;
  12 import jalview.api.AlignViewControllerGuiI;
  13 import jalview.datamodel.AlignmentI;
  14 import jalview.datamodel.DBRefEntry;
  15 import jalview.datamodel.GeneLociI;
  16 import jalview.datamodel.Mapping;
  17 import jalview.datamodel.SequenceFeature;
  18 import jalview.datamodel.SequenceI;
  19 import jalview.ext.ensembl.EnsemblMap;
  20 import jalview.ext.htsjdk.VCFReader;
  21 import jalview.io.gff.Gff3Helper;
  22 import jalview.io.gff.SequenceOntologyI;
  23 import jalview.util.MapList;
  24 import jalview.util.MappingUtils;
  25 import jalview.util.MessageManager;
  26
  27 import java.io.IOException;
  28 import java.util.ArrayList;
  29 import java.util.HashMap;
  30 import java.util.List;
  31 import java.util.Map;
  32 import java.util.Map.Entry;
  33
  34 /**
  35  * A class to read VCF data (using the htsjdk) and add variants as sequence
  36  * features on dna and any related protein product sequences
  37  *
  38  * @author gmcarstairs
  39  */
  40 public class VCFLoader
  41 {
  42   private static final String ALLELE_FREQUENCY_KEY = "AF";
  43
  44   private static final String COMMA = ",";
  45
  46   private static final boolean FEATURE_PER_ALLELE = true;
  47
  48   private static final String FEATURE_GROUP_VCF = "VCF";
  49
  50   private static final String EXCL = "!";
  51
  52   /*
  53    * the alignment we are associated VCF data with
  54    */
  55   private AlignmentI al;
  56
  57   /*
  58    * mappings between VCF and sequence reference assembly regions, as
  59    * key = "species!chromosome!fromAssembly!toAssembly
  60    * value = Map{fromRange, toRange}
  61    */
  62   private Map<String, Map<int[], int[]>> assemblyMappings;
  63
  64   /*
  65    * holds details of the VCF header lines (metadata)
  66    */
  67   private VCFHeader header;
  68
  69   /**
  70    * Constructor given an alignment context
  71    *
  72    * @param alignment
  73    */
  74   public VCFLoader(AlignmentI alignment)
  75   {
  76     al = alignment;
  77
  78     // map of species!chromosome!fromAssembly!toAssembly to {fromRange, toRange}
  79     assemblyMappings = new HashMap<String, Map<int[], int[]>>();
  80   }
  81
  82   /**
  83    * Starts a new thread to query and load VCF variant data on to the alignment
  84    * <p>
  85    * This method is not thread safe - concurrent threads should use separate
  86    * instances of this class.
  87    *
  88    * @param filePath
  89    * @param gui
  90    */
  91   public void loadVCF(final String filePath,
  92           final AlignViewControllerGuiI gui)
  93   {
  94     if (gui != null)
  95     {
  96       gui.setStatus(MessageManager.getString("label.searching_vcf"));
  97     }
  98
  99     new Thread()
 100     {
 101
 102       @Override
 103       public void run()
 104       {
 105         VCFLoader.this.doLoad(filePath, gui);
 106       }
 107
 108     }.start();
 109   }
 110
 111   /**
 112    * Loads VCF on to an alignment - provided it can be related to one or more
 113    * sequence's chromosomal coordinates
 114    *
 115    * @param filePath
 116    * @param gui
 117    *          optional callback handler for messages
 118    */
 119   protected void doLoad(String filePath, AlignViewControllerGuiI gui)
 120   {
 121     VCFReader reader = null;
 122     try
 123     {
 124       // long start = System.currentTimeMillis();
 125       reader = new VCFReader(filePath);
 126
 127       header = reader.getFileHeader();
 128       VCFHeaderLine ref = header
 129               .getOtherHeaderLine(VCFHeader.REFERENCE_KEY);
 130       // check if reference is wrt assembly19 (GRCh37)
 131       // todo may need to allow user to specify reference assembly?
 132       boolean isRefGrch37 = (ref != null && ref.getValue().contains(
 133               "assembly19"));
 134
 135       int varCount = 0;
 136       int seqCount = 0;
 137
 138       /*
 139        * query for VCF overlapping each sequence in turn
 140        */
 141       for (SequenceI seq : al.getSequences())
 142       {
 143         int added = loadVCF(seq, reader, isRefGrch37);
 144         if (added > 0)
 145         {
 146           seqCount++;
 147           varCount += added;
 148           transferAddedFeatures(seq);
 149         }
 150       }
 151       if (gui != null)
 152       {
 153         // long elapsed = System.currentTimeMillis() - start;
 154         String msg = MessageManager.formatMessage("label.added_vcf",
 155                 varCount, seqCount);
 156         gui.setStatus(msg);
 157         if (gui.getFeatureSettingsUI() != null)
 158         {
 159           gui.getFeatureSettingsUI().discoverAllFeatureData();
 160         }
 161       }
 162     } catch (Throwable e)
 163     {
 164       System.err.println("Error processing VCF: " + e.getMessage());
 165       e.printStackTrace();
 166       if (gui != null)
 167       {
 168         gui.setStatus("Error occurred - see console for details");
 169       }
 170     } finally
 171     {
 172       if (reader != null)
 173       {
 174         try
 175         {
 176           reader.close();
 177         } catch (IOException e)
 178         {
 179           // ignore
 180         }
 181       }
 182     }
 183   }
 184
 185   /**
 186    * Transfers VCF features to sequences to which this sequence has a mapping.
 187    * If the mapping is 1:3, computes peptide variants from nucleotide variants.
 188    *
 189    * @param seq
 190    */
 191   protected void transferAddedFeatures(SequenceI seq)
 192   {
 193     DBRefEntry[] dbrefs = seq.getDBRefs();
 194     if (dbrefs == null)
 195     {
 196       return;
 197     }
 198     for (DBRefEntry dbref : dbrefs)
 199     {
 200       Mapping mapping = dbref.getMap();
 201       if (mapping == null || mapping.getTo() == null)
 202       {
 203         continue;
 204       }
 205
 206       SequenceI mapTo = mapping.getTo();
 207       MapList map = mapping.getMap();
 208       if (map.getFromRatio() == 3)
 209       {
 210         /*
 211          * dna-to-peptide product mapping
 212          */
 213         AlignmentUtils.computeProteinFeatures(seq, mapTo, map);
 214       }
 215       else
 216       {
 217         /*
 218          * nucleotide-to-nucleotide mapping e.g. transcript to CDS
 219          */
 220         // TODO no DBRef to CDS is added to transcripts
 221         List<SequenceFeature> features = seq.getFeatures()
 222                 .getPositionalFeatures(SequenceOntologyI.SEQUENCE_VARIANT);
 223         for (SequenceFeature sf : features)
 224         {
 225           if (FEATURE_GROUP_VCF.equals(sf.getFeatureGroup()))
 226           {
 227             transferFeature(sf, mapTo, map);
 228           }
 229         }
 230       }
 231     }
 232   }
 233
 234   /**
 235    * Tries to add overlapping variants read from a VCF file to the given
 236    * sequence, and returns the number of variant features added. Note that this
 237    * requires the sequence to hold information as to its chromosomal positions
 238    * and reference, in order to be able to map the VCF variants to the sequence.
 239    *
 240    * @param seq
 241    * @param reader
 242    * @param isVcfRefGrch37
 243    * @return
 244    */
 245   protected int loadVCF(SequenceI seq, VCFReader reader,
 246           boolean isVcfRefGrch37)
 247   {
 248     int count = 0;
 249     GeneLociI seqCoords = seq.getGeneLoci();
 250     if (seqCoords == null)
 251     {
 252       return 0;
 253     }
 254
 255     List<int[]> seqChromosomalContigs = seqCoords.getMap().getToRanges();
 256     for (int[] range : seqChromosomalContigs)
 257     {
 258       count += addVcfVariants(seq, reader, range, isVcfRefGrch37);
 259     }
 260
 261     return count;
 262   }
 263
 264   /**
 265    * Queries the VCF reader for any variants that overlap the given chromosome
 266    * region of the sequence, and adds as variant features. Returns the number of
 267    * overlapping variants found.
 268    *
 269    * @param seq
 270    * @param reader
 271    * @param range
 272    *          start-end range of a sequence region in its chromosomal
 273    *          coordinates
 274    * @param isVcfRefGrch37
 275    *          true if the VCF is with reference to GRCh37
 276    * @return
 277    */
 278   protected int addVcfVariants(SequenceI seq, VCFReader reader,
 279           int[] range, boolean isVcfRefGrch37)
 280   {
 281     GeneLociI seqCoords = seq.getGeneLoci();
 282
 283     String chromosome = seqCoords.getChromosomeId();
 284     String seqRef = seqCoords.getAssemblyId();
 285     String species = seqCoords.getSpeciesId();
 286
 287     /*
 288      * map chromosomal coordinates from GRCh38 (sequence) to
 289      * GRCh37 (VCF) if necessary
 290      */
 291     // TODO generalise for other assemblies and species
 292     int offset = 0;
 293     String fromRef = "GRCh38";
 294     if (fromRef.equalsIgnoreCase(seqRef) && isVcfRefGrch37)
 295     {
 296       String toRef = "GRCh37";
 297       int[] newRange = mapReferenceRange(range, chromosome, "human",
 298               fromRef, toRef);
 299       if (newRange == null)
 300       {
 301         System.err.println(String.format(
 302                 "Failed to map %s:%s:%s:%d:%d to %s", species, chromosome,
 303                 fromRef, range[0], range[1], toRef));
 304         return 0;
 305       }
 306       offset = newRange[0] - range[0];
 307       range = newRange;
 308     }
 309
 310     boolean forwardStrand = range[0] <= range[1];
 311
 312     /*
 313      * query the VCF for overlaps
 314      * (convert a reverse strand range to forwards)
 315      */
 316     int count = 0;
 317     MapList mapping = seqCoords.getMap();
 318
 319     int fromLocus = Math.min(range[0], range[1]);
 320     int toLocus = Math.max(range[0], range[1]);
 321     CloseableIterator<VariantContext> variants = reader.query(chromosome,
 322             fromLocus, toLocus);
 323     while (variants.hasNext())
 324     {
 325       /*
 326        * get variant location in sequence chromosomal coordinates
 327        */
 328       VariantContext variant = variants.next();
 329
 330       /*
 331        * we can only process SNP variants (which can be reported
 332        * as part of a MIXED variant record
 333        */
 334       if (!variant.isSNP() && !variant.isMixed())
 335       {
 336         continue;
 337       }
 338
 339       int start = variant.getStart() - offset;
 340       int end = variant.getEnd() - offset;
 341
 342       /*
 343        * convert chromosomal location to sequence coordinates
 344        * - null if a partially overlapping feature
 345        */
 346       int[] seqLocation = mapping.locateInFrom(start, end);
 347       if (seqLocation != null)
 348       {
 349         count += addVariantFeature(seq, variant, seqLocation[0],
 350                 seqLocation[1], forwardStrand);
 351       }
 352     }
 353
 354     variants.close();
 355
 356     return count;
 357   }
 358
 359   /**
 360    * Inspects the VCF variant record, and adds variant features to the sequence.
 361    * Only SNP variants are added, not INDELs. Returns the number of features
 362    * added.
 363    * <p>
 364    * If the sequence maps to the reverse strand of the chromosome, reference and
 365    * variant bases are recorded as their complements (C/G, A/T).
 366    *
 367    * @param seq
 368    * @param variant
 369    * @param featureStart
 370    * @param featureEnd
 371    * @param forwardStrand
 372    */
 373   protected int addVariantFeature(SequenceI seq, VariantContext variant,
 374           int featureStart, int featureEnd, boolean forwardStrand)
 375   {
 376     byte[] reference = variant.getReference().getBases();
 377     if (reference.length != 1)
 378     {
 379       /*
 380        * sorry, we don't handle INDEL variants
 381        */
 382       return 0;
 383     }
 384
 385     if (FEATURE_PER_ALLELE)
 386     {
 387       return addAlleleFeatures(seq, variant, featureStart, featureEnd,
 388               forwardStrand);
 389     }
 390
 391     /*
 392      * for now we extract allele frequency as feature score; note
 393      * this attribute is String for a simple SNP, but List<String> if
 394      * multiple alleles at the locus; we extract for the simple case only
 395      */
 396     float score = getAlleleFrequency(variant, 0);
 397
 398     StringBuilder sb = new StringBuilder();
 399     sb.append(forwardStrand ? (char) reference[0] : complement(reference));
 400
 401     /*
 402      * inspect alleles and record SNP variants (as the variant
 403      * record could be MIXED and include INDEL and SNP alleles)
 404      * warning: getAlleles gives no guarantee as to the order
 405      * in which they are returned
 406      */
 407     for (Allele allele : variant.getAlleles())
 408     {
 409       if (!allele.isReference())
 410       {
 411         byte[] alleleBase = allele.getBases();
 412         if (alleleBase.length == 1)
 413         {
 414           sb.append(COMMA).append(
 415                   forwardStrand ? (char) alleleBase[0]
 416                           : complement(alleleBase));
 417         }
 418       }
 419     }
 420     String alleles = sb.toString(); // e.g. G,A,C
 421
 422     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 423
 424     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 425             featureEnd, score, FEATURE_GROUP_VCF);
 426
 427     sf.setValue(Gff3Helper.ALLELES, alleles);
 428
 429     Map<String, Object> atts = variant.getAttributes();
 430     for (Entry<String, Object> att : atts.entrySet())
 431     {
 432       sf.setValue(att.getKey(), att.getValue());
 433     }
 434     seq.addSequenceFeature(sf);
 435
 436     return 1;
 437   }
 438
 439   /**
 440    * A convenience method to get the AF value for the given alternate allele
 441    * index
 442    *
 443    * @param variant
 444    * @param alleleIndex
 445    * @return
 446    */
 447   protected float getAlleleFrequency(VariantContext variant, int alleleIndex)
 448   {
 449     float score = 0f;
 450     String attributeValue = getAttributeValue(variant,
 451             ALLELE_FREQUENCY_KEY, alleleIndex);
 452     if (attributeValue != null)
 453     {
 454       try
 455       {
 456         score = Float.parseFloat(attributeValue);
 457       } catch (NumberFormatException e)
 458       {
 459         // leave as 0
 460       }
 461     }
 462
 463     return score;
 464   }
 465
 466   /**
 467    * A convenience method to get an attribute value for an alternate allele
 468    *
 469    * @param variant
 470    * @param attributeName
 471    * @param alleleIndex
 472    * @return
 473    */
 474   protected String getAttributeValue(VariantContext variant,
 475           String attributeName, int alleleIndex)
 476   {
 477     Object att = variant.getAttribute(attributeName);
 478
 479     if (att instanceof String)
 480     {
 481       return (String) att;
 482     }
 483     else if (att instanceof ArrayList)
 484     {
 485       return ((List<String>) att).get(alleleIndex);
 486     }
 487
 488     return null;
 489   }
 490
 491   /**
 492    * Adds one variant feature for each SNP allele in the VCF variant record, and
 493    * returns the number of features added.
 494    *
 495    * @param seq
 496    * @param variant
 497    * @param featureStart
 498    * @param featureEnd
 499    * @param forwardStrand
 500    * @return
 501    */
 502   protected int addAlleleFeatures(SequenceI seq, VariantContext variant,
 503           int featureStart, int featureEnd, boolean forwardStrand)
 504   {
 505     int added = 0;
 506
 507     /*
 508      * Javadoc says getAlternateAlleles() imposes no order on the list returned
 509      * so we proceed defensively to get them in strict order
 510      */
 511     int altAlleleCount = variant.getAlternateAlleles().size();
 512     for (int i = 0; i < altAlleleCount; i++)
 513     {
 514       added += addAlleleFeature(seq, variant, i, featureStart, featureEnd,
 515               forwardStrand);
 516     }
 517     return added;
 518   }
 519
 520   /**
 521    * Inspects one allele and attempts to add a variant feature for it to the
 522    * sequence. Only SNP variants are added as features. We extract as much as
 523    * possible of the additional data associated with this allele to store in the
 524    * feature's key-value map. Answers the number of features added (0 or 1).
 525    *
 526    * @param seq
 527    * @param variant
 528    * @param altAlleleIndex
 529    * @param featureStart
 530    * @param featureEnd
 531    * @param forwardStrand
 532    * @return
 533    */
 534   protected int addAlleleFeature(SequenceI seq, VariantContext variant,
 535           int altAlleleIndex, int featureStart, int featureEnd,
 536           boolean forwardStrand)
 537   {
 538     byte[] reference = variant.getReference().getBases();
 539     Allele alt = variant.getAlternateAllele(altAlleleIndex);
 540     byte[] allele = alt.getBases();
 541     if (allele.length != 1)
 542     {
 543       /*
 544        * not a SNP variant
 545        */
 546       return 0;
 547     }
 548
 549     /*
 550      * build the ref,alt allele description e.g. "G,A"
 551      */
 552     StringBuilder sb = new StringBuilder();
 553     sb.append(forwardStrand ? (char) reference[0] : complement(reference));
 554     sb.append(COMMA);
 555     sb.append(forwardStrand ? (char) allele[0] : complement(allele));
 556     String alleles = sb.toString(); // e.g. G,A
 557
 558     String type = SequenceOntologyI.SEQUENCE_VARIANT;
 559     float score = getAlleleFrequency(variant, altAlleleIndex);
 560
 561     SequenceFeature sf = new SequenceFeature(type, alleles, featureStart,
 562             featureEnd, score, FEATURE_GROUP_VCF);
 563
 564     sf.setValue(Gff3Helper.ALLELES, alleles);
 565
 566     addAlleleProperties(variant, sf, altAlleleIndex);
 567
 568     seq.addSequenceFeature(sf);
 569
 570     return 1;
 571   }
 572
 573   /**
 574    * Add any allele-specific VCF key-value data to the sequence feature
 575    *
 576    * @param variant
 577    * @param sf
 578    * @param altAlelleIndex
 579    */
 580   protected void addAlleleProperties(VariantContext variant,
 581           SequenceFeature sf, final int altAlelleIndex)
 582   {
 583     Map<String, Object> atts = variant.getAttributes();
 584
 585     /*
 586      * process variant data, extracting values which are allele-specific
 587      * these may be per alternate allele (INFO[key].Number = 'A')
 588      * or per allele including reference (INFO[key].Number = 'R')
 589      */
 590     for (Entry<String, Object> att : atts.entrySet())
 591     {
 592       String key = att.getKey();
 593       VCFHeaderLineCount number = header.getInfoHeaderLine(key)
 594               .getCountType();
 595       int index = altAlelleIndex;
 596       if (number == VCFHeaderLineCount.R)
 597       {
 598         /*
 599          * one value per allele including reference, so bump index
 600          * e.g. the 3rd value is for the  2nd alternate allele
 601          */
 602         index++;
 603       }
 604       /*
 605        * CSQ behaves as if Number=A but declares as Number=.
 606        * so give it special treatment
 607        */
 608       else if (!"CSQ".equals(key) && number != VCFHeaderLineCount.A)
 609       {
 610         /*
 611          * don't save other values as not allele-related
 612          */
 613         continue;
 614       }
 615
 616       /*
 617        * take the index'th value
 618        */
 619       String value = getAttributeValue(variant, key, index);
 620       if (value != null)
 621       {
 622         sf.setValue(key, value);
 623       }
 624     }
 625   }
 626
 627   /**
 628    * A convenience method to complement a dna base and return the string value
 629    * of its complement
 630    *
 631    * @param reference
 632    * @return
 633    */
 634   protected String complement(byte[] reference)
 635   {
 636     return String.valueOf(Dna.getComplement((char) reference[0]));
 637   }
 638
 639   /**
 640    * Determines the location of the query range (chromosome positions) in a
 641    * different reference assembly.
 642    * <p>
 643    * If the range is just a subregion of one for which we already have a mapping
 644    * (for example, an exon sub-region of a gene), then the mapping is just
 645    * computed arithmetically.
 646    * <p>
 647    * Otherwise, calls the Ensembl REST service that maps from one assembly
 648    * reference's coordinates to another's
 649    *
 650    * @param queryRange
 651    *          start-end chromosomal range in 'fromRef' coordinates
 652    * @param chromosome
 653    * @param species
 654    * @param fromRef
 655    *          assembly reference for the query coordinates
 656    * @param toRef
 657    *          assembly reference we wish to translate to
 658    * @return the start-end range in 'toRef' coordinates
 659    */
 660   protected int[] mapReferenceRange(int[] queryRange, String chromosome,
 661           String species, String fromRef, String toRef)
 662   {
 663     /*
 664      * first try shorcut of computing the mapping as a subregion of one
 665      * we already have (e.g. for an exon, if we have the gene mapping)
 666      */
 667     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
 668             species, fromRef, toRef);
 669     if (mappedRange != null)
 670     {
 671       return mappedRange;
 672     }
 673
 674     /*
 675      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
 676      */
 677     EnsemblMap mapper = new EnsemblMap();
 678     int[] mapping = mapper.getMapping(species, chromosome, fromRef, toRef,
 679             queryRange);
 680
 681     if (mapping == null)
 682     {
 683       // mapping service failure
 684       return null;
 685     }
 686
 687     /*
 688      * save mapping for possible future re-use
 689      */
 690     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 691     if (!assemblyMappings.containsKey(key))
 692     {
 693       assemblyMappings.put(key, new HashMap<int[], int[]>());
 694     }
 695
 696     assemblyMappings.get(key).put(queryRange, mapping);
 697
 698     return mapping;
 699   }
 700
 701   /**
 702    * If we already have a 1:1 contiguous mapping which subsumes the given query
 703    * range, this method just calculates and returns the subset of that mapping,
 704    * else it returns null. In practical terms, if a gene has a contiguous
 705    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
 706    * subsidiary exons occupy unchanged relative positions, and just compute
 707    * these as offsets, rather than do another lookup of the mapping.
 708    * <p>
 709    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
 710    * simply remove this method or let it always return null.
 711    * <p>
 712    * Warning: many rapid calls to the /map service map result in a 429 overload
 713    * error response
 714    *
 715    * @param queryRange
 716    * @param chromosome
 717    * @param species
 718    * @param fromRef
 719    * @param toRef
 720    * @return
 721    */
 722   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
 723           String species, String fromRef, String toRef)
 724   {
 725     String key = makeRangesKey(chromosome, species, fromRef, toRef);
 726     if (assemblyMappings.containsKey(key))
 727     {
 728       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
 729       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
 730       {
 731         int[] fromRange = mappedRange.getKey();
 732         int[] toRange = mappedRange.getValue();
 733         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
 734         {
 735           /*
 736            * mapping is 1:1 in length, so we trust it to have no discontinuities
 737            */
 738           if (MappingUtils.rangeContains(fromRange, queryRange))
 739           {
 740             /*
 741              * fromRange subsumes our query range
 742              */
 743             int offset = queryRange[0] - fromRange[0];
 744             int mappedRangeFrom = toRange[0] + offset;
 745             int mappedRangeTo = mappedRangeFrom + (queryRange[1] - queryRange[0]);
 746             return new int[] { mappedRangeFrom, mappedRangeTo };
 747           }
 748         }
 749       }
 750     }
 751     return null;
 752   }
 753
 754   /**
 755    * Transfers the sequence feature to the target sequence, locating its start
 756    * and end range based on the mapping. Features which do not overlap the
 757    * target sequence are ignored.
 758    *
 759    * @param sf
 760    * @param targetSequence
 761    * @param mapping
 762    *          mapping from the feature's coordinates to the target sequence
 763    */
 764   protected void transferFeature(SequenceFeature sf,
 765           SequenceI targetSequence, MapList mapping)
 766   {
 767     int[] mappedRange = mapping.locateInTo(sf.getBegin(), sf.getEnd());
 768
 769     if (mappedRange != null)
 770     {
 771       String group = sf.getFeatureGroup();
 772       int newBegin = Math.min(mappedRange[0], mappedRange[1]);
 773       int newEnd = Math.max(mappedRange[0], mappedRange[1]);
 774       SequenceFeature copy = new SequenceFeature(sf, newBegin, newEnd,
 775               group, sf.getScore());
 776       targetSequence.addSequenceFeature(copy);
 777     }
 778   }
 779
 780   /**
 781    * Formats a ranges map lookup key
 782    *
 783    * @param chromosome
 784    * @param species
 785    * @param fromRef
 786    * @param toRef
 787    * @return
 788    */
 789   protected static String makeRangesKey(String chromosome, String species,
 790           String fromRef, String toRef)
 791   {
 792     return species + EXCL + chromosome + EXCL + fromRef + EXCL
 793             + toRef;
 794   }
 795 }