feature: Read genetic distance data from the VCF directly. No map file necessary.

src/hapibd/HapIbdPar.java

@@ -38,6 +38,7 @@ public final class HapIbdPar {
     // data input/output parameters
     private final File gt;
     private final File map;
+    private final boolean useVcfMap;
     private final String out;
     private final File excludesamples;
 
@@ -79,7 +80,8 @@ public HapIbdPar(String[] args) {
         // data input/output parameters
         gt = Validate.getFile(Validate.stringArg("gt", argsMap, true, null,
                 null));
-        map = Validate.getFile(Validate.stringArg("map", argsMap, true, null, null));
+        useVcfMap = Validate.booleanArg("vcf-has-cm", argsMap, false, false);
+        map = Validate.getFile(Validate.stringArg("map", argsMap, false, null, null));
         out = Validate.stringArg("out", argsMap, true, null, null);
         excludesamples = Validate.getFile(Validate.stringArg("excludesamples",
                 argsMap, false, null, null));
@@ -156,6 +158,14 @@ public File map() {
         return map;
     }
 
+    /**
+     * Returns if we should get genetic position from the VCF
+     * @return useVcfMap boolean
+     */
+    public boolean useVcfMap() {
+        return useVcfMap;
+    }
+
     /**
      * Returns the out parameter.
      * @return the out parameter

src/hapibd/PbwtIbdDriver.java

@@ -62,8 +62,7 @@ public final class PbwtIbdDriver {
      * @throws NullPointerException if {@code par == null}
      */
     public static long[] detectIbd(HapIbdPar par) {
-        ChromInterval chromInt = null;
-        GeneticMap genMap = GeneticMap.geneticMap(par.map(), chromInt);
+        GeneticMap genMap = GeneticMap.geneticMap(par, null);
         long[] nSamplesAndMarkers = new long[2];
         File hbdFile = new File(par.out() + ".hbd.gz");
         File ibdFile = new File(par.out() + ".ibd.gz");

src/vcf/BasicMarker.java

@@ -48,6 +48,8 @@ public class BasicMarker implements Marker {
     private final int nGenotypes;
     private final int end;
 
+    private final double geneticPos;
+
     /**
      * Constructs a new {@code BasicMarker} instance from the specified data.
      * The {@code end()} method of the new instance will return {@code -1}.
@@ -69,7 +71,7 @@ public class BasicMarker implements Marker {
      * of {@code alleles} is {@code null}
      */
     public BasicMarker(int chrom, int pos, String[] ids, String[] alleles) {
-        this(chrom, pos, ids, alleles, -1);
+        this(chrom, pos, ids, alleles, -1, Double.NaN);
     }
 
     /**
@@ -93,7 +95,7 @@ public BasicMarker(int chrom, int pos, String[] ids, String[] alleles) {
      * of {@code alleles} is {@code null}
      */
     public BasicMarker(int chrom, int pos, String[] ids, String[] alleles,
-            int end) {
+            int end, double geneticPos) {
         if (chrom<0 || chrom>=ChromIds.instance().size()) {
             throw new IndexOutOfBoundsException(String.valueOf(chrom));
         }
@@ -106,6 +108,7 @@ public BasicMarker(int chrom, int pos, String[] ids, String[] alleles,
         this.alleles = canonicalAlleles(alleles.clone());
         this.nGenotypes = (alleles.length*(alleles.length+1))/2;
         this.end = end;
+        this.geneticPos = geneticPos;
     }
 
     private static String[] removeMissingIds(int chrom, int pos, String[] ids) {
@@ -253,6 +256,7 @@ public BasicMarker(String vcfRecord) {
         this.alleles = extractAlleles(chromIndex, pos, fields[3], fields[4]);
         this.nGenotypes = alleles.length*(alleles.length+1)/2;
         this.end = extractEnd(chromIndex, pos, infoField);
+        this.geneticPos = extractGeneticPos(chromIndex, pos, infoField);
     }
 
     private static int extractChrom(String chrom, String vcfRecord) {
@@ -317,6 +321,26 @@ private static String[] extractAlleles(int chrom, int pos, String ref,
         return canonicalAlleles(alleles);
     }
 
+    private static double extractGeneticPos(int chrom, int pos, String info) {
+        String[] fields = StringUtil.getFields(info, Const.semicolon);
+        String key = "CM=";
+        for (String field : fields) {
+            if (field.startsWith(key)) {
+                String value = field.substring(3);
+                try {
+                    return Double.parseDouble(value);
+                } catch (NumberFormatException e) {
+                    String s = "ERROR: invalid INFO:CM field at "
+                            + coordinate(chrom, pos) + " [" + key + value
+                            + "]";
+                    Utilities.exit(s);
+                }
+            }
+        }
+
+        return Double.NaN;
+    }
+
     /*
      * Returns value of first END key in the specified INFO field, or
      * returns -1 if there is no END key in INFO field.
@@ -361,6 +385,7 @@ private BasicMarker(Marker markerOnReverseStrand) {
         Marker m = markerOnReverseStrand;
         this.chromIndex = m.chromIndex();
         this.pos = m.pos();
+        this.geneticPos = m.geneticPos();
         this.ids = new String[m.nIds()];
         for (int j=0; j<this.ids.length; ++j) {
             this.ids[j] = m.id(j);
@@ -424,6 +449,11 @@ public int pos() {
         return pos;
     }
 
+    @Override
+    public double geneticPos() {
+        return geneticPos;
+    }
+
     @Override
     public int nIds() {
         return ids.length;
@@ -492,6 +522,11 @@ public String toString() {
                 sb.append(alleles[j]);
             }
         }
+        if (!Double.isNaN(geneticPos)) {
+            sb.append(Const.tab);
+            sb.append("CM=");
+            sb.append(geneticPos);
+        }
         return sb.toString();
     }
 

src/vcf/GeneticMap.java

@@ -18,6 +18,8 @@
 package vcf;
 
 import beagleutil.ChromInterval;
+import hapibd.HapIbdPar;
+
 import java.io.File;
 
 /**
@@ -109,17 +111,21 @@ public interface GeneticMap {
      * @throws IllegalArgumentException if all base positions on a chromosome
      * have the same genetic map position
      */
-    static GeneticMap geneticMap(File file, ChromInterval chromInt) {
-        if (file==null) {
+    static GeneticMap geneticMap(HapIbdPar par, ChromInterval chromInt) {
+        if (par.useVcfMap()) {
+            return new VcfGeneticMap();
+        }
+
+        if (par.map()==null) {
             double scaleFactor = 1e-6;
             return new PositionMap(scaleFactor);
         }
         else {
             if (chromInt==null) {
-                return PlinkGenMap.fromPlinkMapFile(file);
+                return PlinkGenMap.fromPlinkMapFile(par.map());
             }
             else {
-                return PlinkGenMap.fromPlinkMapFile(file, chromInt.chrom());
+                return PlinkGenMap.fromPlinkMapFile(par.map(), chromInt.chrom());
             }
         }
     }

src/vcf/Marker.java

@@ -45,6 +45,12 @@ public interface Marker extends Comparable<Marker> {
      */
     int pos();
 
+    /**
+     * Returns the genetic position
+     * @return the genetic position in CM
+     */
+    double geneticPos();
+
     /**
      * Returns the number of marker identifiers.
      * @return the number of marker identifiers

src/vcf/VcfGeneticMap.java

@@ -0,0 +1,55 @@
+package vcf;
+
+public class VcfGeneticMap implements GeneticMap {
+    /**
+     * Returns the base position corresponding to the specified genetic map
+     * position. If the genetic position is not a map position then the base
+     * position is estimated from the nearest genetic map positions using
+     * linear interpolation.
+     *
+     * @param chrom           the chromosome index
+     * @param geneticPosition the genetic position on the chromosome
+     * @return the base position corresponding to the specified genetic map
+     * position
+     * @throws UnsupportedOperationException Not implemented for VcfGeneticMap
+     */
+    @Override
+    public int basePos(int chrom, double geneticPosition) {
+        // VcfGeneticMap cannot go from genetic to physical position
+        // At this time, this method is not used
+        throw new UnsupportedOperationException("VcfGeneticMap cannot map from genetic to physical position");
+    }
+
+    /**
+     * Returns the genetic map position of the specified marker. The
+     * genetic map position is estimated using linear interpolation.
+     *
+     * @param marker a genetic marker
+     * @return the genetic map position of the specified marker
+     * @throws IllegalArgumentException if the vcf does not have a position mapped for this marker
+     * @throws NullPointerException     if {@code marker == null}
+     */
+    @Override
+    public double genPos(Marker marker) {
+        double geneticPos = marker.geneticPos();
+        if (Double.isNaN(geneticPos)) {
+            throw new IllegalArgumentException("Called with vcf-has-cm, but marker at position " + marker.pos() +
+                    " does not have genetic distance in INFO column.");
+        }
+        return marker.geneticPos();
+    }
+
+    /**
+     * Returns the genetic map position of the specified genome coordinate.
+     * The genetic map position is estimated using linear interpolation.
+     *
+     * @param chrom        the chromosome index
+     * @param basePosition the base coordinate on the chromosome
+     * @return the genetic map position of the specified genome coordinate
+     * @throws UnsupportedOperationException Not implemented for VcfGeneticMap
+     */
+    @Override
+    public double genPos(int chrom, int basePosition) {
+        throw new UnsupportedOperationException("VcfMap cannot interpolate positions.");
+    }
+}