JAL-2738 extract utility method to MappingUtils

[jalview.git] / src / jalview / io / vcf / VCFLoader.java

diff --git a/src/jalview/io/vcf/VCFLoader.java b/src/jalview/io/vcf/VCFLoader.java
index 48f55d4..4adc97c 100644 (file)
--- a/src/jalview/io/vcf/VCFLoader.java
+++ b/src/jalview/io/vcf/VCFLoader.java
@@ -7,6 +7,7 @@ import htsjdk.variant.vcf.VCFHeader;
 import htsjdk.variant.vcf.VCFHeaderLine;
 
 import jalview.analysis.AlignmentUtils;
+import jalview.api.AlignViewControllerGuiI;
 import jalview.datamodel.AlignmentI;
 import jalview.datamodel.DBRefEntry;
 import jalview.datamodel.GeneLoci;
@@ -18,7 +19,9 @@ import jalview.ext.ensembl.EnsemblMap;
 import jalview.ext.htsjdk.VCFReader;
 import jalview.io.gff.SequenceOntologyI;
 import jalview.util.MapList;
+import jalview.util.MappingUtils;
 
+import java.io.IOException;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
@@ -67,14 +70,15 @@ public class VCFLoader
    * instances of this class.
    * 
    * @param filePath
+   * @param status
    */
-  public void loadVCF(String filePath)
+  public void loadVCF(String filePath, AlignViewControllerGuiI status)
   {
     VCFReader reader = null;
 
     try
     {
-      long start = System.currentTimeMillis();
+      // long start = System.currentTimeMillis();
       reader = new VCFReader(filePath);
 
       VCFHeader header = reader.getFileHeader();
@@ -101,16 +105,29 @@ public class VCFLoader
           computePeptideVariants(seq);
         }
       }
-      long elapsed = System.currentTimeMillis() - start;
-      System.out.println(String.format(
-              "Added %d VCF variants to %d sequence(s) (%dms)", varCount,
-              seqCount, elapsed));
-
-      reader.close();
+      // long elapsed = System.currentTimeMillis() - start;
+      String msg = String.format("Added %d VCF variants to %d sequence(s)",
+              varCount, seqCount);
+      if (status != null)
+      {
+        status.setStatus(msg);
+      }
     } catch (Throwable e)
     {
       System.err.println("Error processing VCF: " + e.getMessage());
       e.printStackTrace();
+    } finally
+    {
+      if (reader != null)
+      {
+        try
+        {
+          reader.close();
+        } catch (IOException e)
+        {
+          // ignore
+        }
+      }
     }
   }
 
@@ -164,8 +181,7 @@ public class VCFLoader
       return 0;
     }
 
-    MapList mapping = seqCoords.getMapping();
-    List<int[]> seqChromosomalContigs = mapping.getToRanges();
+    List<int[]> seqChromosomalContigs = seqCoords.mapping.getToRanges();
     for (int[] range : seqChromosomalContigs)
     {
       count += addVcfVariants(seq, reader, range, isVcfRefGrch37);
@@ -193,9 +209,9 @@ public class VCFLoader
   {
     GeneLoci seqCoords = ((Sequence) seq).getGeneLoci();
 
-    String chromosome = seqCoords.getChromosome();
-    String seqRef = seqCoords.getReference();
-    String species = seqCoords.getSpecies();
+    String chromosome = seqCoords.chromosome;
+    String seqRef = seqCoords.assembly;
+    String species = seqCoords.species;
 
     // TODO handle species properly
     if ("".equals(species))
@@ -228,12 +244,15 @@ public class VCFLoader
 
     /*
      * query the VCF for overlaps
+     * (convert a reverse strand range to forwards)
      */
     int count = 0;
-    MapList mapping = seqCoords.getMapping();
+    MapList mapping = seqCoords.mapping;
 
+    int fromLocus = Math.min(range[0], range[1]);
+    int toLocus = Math.max(range[0], range[1]);
     CloseableIterator<VariantContext> variants = reader.query(chromosome,
-            range[0], range[1]);
+            fromLocus, toLocus);
     while (variants.hasNext())
     {
       /*
@@ -286,6 +305,11 @@ public class VCFLoader
     String alleles = sb.toString(); // e.g. G,A,C
 
     String type = SequenceOntologyI.SEQUENCE_VARIANT;
+
+    /*
+     * extract allele frequency as feature score, but only if
+     * a simple SNP (not for >1 co-located SNPs as each has a score)
+     */
     float score = 0f;
     if (alleleCount == 1)
     {
@@ -368,8 +392,7 @@ public class VCFLoader
     /*
      * save mapping for possible future re-use
      */
-    String key = species + EXCL + chromosome + EXCL + fromRef + EXCL
-            + toRef;
+    String key = makeRangesKey(chromosome, species, fromRef, toRef);
     if (!assemblyMappings.containsKey(key))
     {
       assemblyMappings.put(key, new HashMap<int[], int[]>());
@@ -404,8 +427,7 @@ public class VCFLoader
   protected int[] findSubsumedRangeMapping(int[] queryRange, String chromosome,
           String species, String fromRef, String toRef)
   {
-    String key = species + EXCL + chromosome + EXCL + fromRef + EXCL
-            + toRef;
+    String key = makeRangesKey(chromosome, species, fromRef, toRef);
     if (assemblyMappings.containsKey(key))
     {
       Map<int[], int[]> mappedRanges = assemblyMappings.get(key);
@@ -418,7 +440,7 @@ public class VCFLoader
           /*
            * mapping is 1:1 in length, so we trust it to have no discontinuities
            */
-          if (fromRange[0] <= queryRange[0] && fromRange[1] >= queryRange[1])
+          if (MappingUtils.rangeContains(fromRange, queryRange))
           {
             /*
              * fromRange subsumes our query range
@@ -433,4 +455,20 @@ public class VCFLoader
     }
     return null;
   }
+
+  /**
+   * Formats a ranges map lookup key
+   * 
+   * @param chromosome
+   * @param species
+   * @param fromRef
+   * @param toRef
+   * @return
+   */
+  protected static String makeRangesKey(String chromosome, String species,
+          String fromRef, String toRef)
+  {
+    return species + EXCL + chromosome + EXCL + fromRef + EXCL
+            + toRef;
+  }
 }