src/jalview/datamodel/GenomicAssemblies.java

   1 package jalview.datamodel;
   2
   3 import java.util.HashMap;
   4 import java.util.Map;
   5 import java.util.Map.Entry;
   6
   7 import jalview.ext.ensembl.EnsemblMap;
   8 import jalview.io.vcf.VCFLoader;
   9 import jalview.util.MappingUtils;
  10
  11 public class GenomicAssemblies
  12 {
  13
  14   /*
  15    * mappings between VCF and sequence reference assembly regions, as
  16    * key = "species!chromosome!fromAssembly!toAssembly
  17    * value = Map{fromRange, toRange}
  18    */
  19   private static Map<String, Map<int[], int[]>> assemblyMappings = new HashMap<>();
  20
  21   /*
  22    * internal delimiter used to build keys for assemblyMappings
  23    *
  24    */
  25   private static final String EXCL = "!";
  26
  27   /**
  28    * Formats a ranges map lookup key
  29    *
  30    * @param chromosome
  31    * @param species
  32    * @param fromRef
  33    * @param toRef
  34    * @return
  35    */
  36   public static String makeRangesKey(String chromosome, String species,
  37           String fromRef, String toRef)
  38   {
  39     return species + EXCL + chromosome + EXCL + fromRef + EXCL + toRef;
  40   }
  41
  42   /**
  43    * Determines the location of the query range (chromosome positions) in a
  44    * different reference assembly.
  45    * <p>
  46    * If the range is just a subregion of one for which we already have a mapping
  47    * (for example, an exon sub-region of a gene), then the mapping is just
  48    * computed arithmetically.
  49    * <p>
  50    * Otherwise, calls the Ensembl REST service that maps from one assembly
  51    * reference's coordinates to another's
  52    *
  53    * @param queryRange
  54    *          start-end chromosomal range in 'fromRef' coordinates
  55    * @param chromosome
  56    * @param species
  57    * @param fromRef
  58    *          assembly reference for the query coordinates
  59    * @param toRef
  60    *          assembly reference we wish to translate to
  61    * @return the start-end range in 'toRef' coordinates
  62    */
  63   public static int[] mapReferenceRange(int[] queryRange, String chromosome,
  64           String species, String fromRef, String toRef)
  65   {
  66     /*
  67      * first try shorcut of computing the mapping as a subregion of one
  68      * we already have (e.g. for an exon, if we have the gene mapping)
  69      */
  70     int[] mappedRange = findSubsumedRangeMapping(queryRange, chromosome,
  71             species, fromRef, toRef);
  72     if (mappedRange != null)
  73     {
  74       return mappedRange;
  75     }
  76
  77     /*
  78      * call (e.g.) http://rest.ensembl.org/map/human/GRCh38/17:45051610..45109016:1/GRCh37
  79      */
  80     EnsemblMap mapper = new EnsemblMap();
  81     int[] mapping = mapper.getAssemblyMapping(species, chromosome, fromRef,
  82             toRef, queryRange);
  83
  84     if (mapping == null)
  85     {
  86       // mapping service failure
  87       return null;
  88     }
  89
  90     /*
  91      * save mapping for possible future re-use
  92      */
  93     String key = GenomicAssemblies.makeRangesKey(chromosome, species, fromRef, toRef);
  94     if (!assemblyMappings.containsKey(key))
  95     {
  96       assemblyMappings.put(key, new HashMap<int[], int[]>());
  97     }
  98
  99     assemblyMappings.get(key).put(queryRange, mapping);
 100
 101     return mapping;
 102   }
 103
 104   /**
 105    * If we already have a 1:1 contiguous mapping which subsumes the given query
 106    * range, this method just calculates and returns the subset of that mapping,
 107    * else it returns null. In practical terms, if a gene has a contiguous
 108    * mapping between (for example) GRCh37 and GRCh38, then we assume that its
 109    * subsidiary exons occupy unchanged relative positions, and just compute
 110    * these as offsets, rather than do another lookup of the mapping.
 111    * <p>
 112    * If in future these assumptions prove invalid (e.g. for bacterial dna?!),
 113    * simply remove this method or let it always return null.
 114    * <p>
 115    * Warning: many rapid calls to the /map service map result in a 429 overload
 116    * error response
 117    *
 118    * @param queryRange
 119    * @param chromosome
 120    * @param species
 121    * @param fromRef
 122    * @param toRef
 123    * @return
 124    */
 125   protected static int[] findSubsumedRangeMapping(int[] queryRange,
 126           String chromosome, String species, String fromRef, String toRef)
 127   {
 128     String key = GenomicAssemblies.makeRangesKey(chromosome, species, fromRef, toRef);
 129     if (assemblyMappings.containsKey(key))
 130     {
 131       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
 132       for (Entry<int[], int[]> mappedRange : mappedRanges.entrySet())
 133       {
 134         int[] fromRange = mappedRange.getKey();
 135         int[] toRange = mappedRange.getValue();
 136         if (fromRange[1] - fromRange[0] == toRange[1] - toRange[0])
 137         {
 138           /*
 139            * mapping is 1:1 in length, so we trust it to have no discontinuities
 140            */
 141           if (MappingUtils.rangeContains(fromRange, queryRange))
 142           {
 143             /*
 144              * fromRange subsumes our query range
 145              */
 146             int offset = queryRange[0] - fromRange[0];
 147             int mappedRangeFrom = toRange[0] + offset;
 148             int mappedRangeTo = mappedRangeFrom
 149                     + (queryRange[1] - queryRange[0]);
 150             return new int[] { mappedRangeFrom, mappedRangeTo };
 151           }
 152         }
 153       }
 154     }
 155     return null;
 156   }
 157
 158 }