beagleutil/ChromIds.java

@@ -18,10 +18,11 @@
  */
 package beagleutil;
 
+import java.util.ArrayList;
 import java.util.Arrays;
+import java.util.List;
 import java.util.concurrent.ConcurrentHashMap;
 import java.util.concurrent.ConcurrentMap;
-import java.util.concurrent.CopyOnWriteArrayList;
 
 /**
  * <p>Class {@code ChromIds} is a singleton class that represents a
@@ -38,7 +39,7 @@ public final class ChromIds {
     private final ThreadSafeIndexer<String> indexer;
 
     private final ConcurrentMap<String, Integer> map;
-    private volatile CopyOnWriteArrayList<String> ids;
+    private volatile List<String> ids;
 
     private ChromIds() {
         // private constructor to restrict instantiation.
@@ -46,7 +47,7 @@ public final class ChromIds {
         this.indexer = new ThreadSafeIndexer<>(initCapacity);
 
         this.map = new ConcurrentHashMap<>(initCapacity);
-        this.ids = new CopyOnWriteArrayList<>(new String[0]);
+        this.ids = new ArrayList<>();
     }
 
     /**
@@ -149,9 +150,7 @@ public final class ChromIds {
      */
     public String id(int index) {
         if (index >= ids.size()) {
-            for (int j=ids.size(), n=indexer.size(); j<n; ++j) {
-                ids.add(indexer.item(j));
-            }
+            ids = indexer.items();
         }
         return ids.get(index);
     }

main/Main.java

@@ -68,8 +68,8 @@ public class Main {
      * The program name and version.
      */
     public static final String VERSION = "(version 5.0)";
-    public static final String PROGRAM = "beagle.28Sep18.793.jar";
-    public static final String COMMAND = "java -jar beagle.28Sep18.793.jar";
+    public static final String PROGRAM = "beagle.03Jul19.b33.jar";
+    public static final String COMMAND = "java -jar beagle.03Jul19.b33.jar";
 
     /**
      * The copyright string.
@@ -81,7 +81,7 @@ public class Main {
      */
     public static final String SHORT_HELP = Main.PROGRAM + " " + VERSION
             + Const.nl + Main.COPYRIGHT
-            + Const.nl + "Enter \"java -jar beagle.28Sep18.793.jar\" to "
+            + Const.nl + "Enter \"java -jar beagle.03Jul19.b33.jar\" to "
             + "list command line argument";
 
     private final Par par;

vcf/AllData.java

@@ -130,8 +130,7 @@ public class AllData implements Data {
 
     private static GT targGTWindow(Samples samples, GTRec[] targData,
             Pedigree ped) {
-        GT gt = new BasicGT(samples, targData);
-        return gt.isGTData() ? new XBasicGT(gt, ped) : gt;
+        return new BasicGT(samples, targData);
     }
 
     @Override
@@ -214,8 +213,7 @@ public class AllData implements Data {
         for (int j=0; j<refMarkerIndex.length; ++j) {
             restricted[j] = vma[refMarkerIndex[j]];
         }
-        GT gt = new BasicGT(samples, restricted);
-        return gt.isGTData() ? new XBasicGT(gt, ped) : gt;
+        return new BasicGT(samples, restricted);
     }
 
     private void setRefHaplotypes(RefGT refGT) {

vcf/RestrictedVcfWindow.java

@@ -98,10 +98,12 @@ public class RestrictedVcfWindow  implements Closeable {
         for (int j = fullOverlap, n=nextMarkers.nMarkers(); j<n; ++j) {
             Marker m = nextMarkers.marker(j);
             if (next!=null && next.marker().chromIndex()==m.chromIndex()) {
-                while (next != null && next.marker().pos() < m.pos()) {
+                while (next != null && next.marker().chromIndex()==m.chromIndex()
+                        && next.marker().pos()<m.pos()) {
                     next = it.hasNext() ? it.next() : null;
                 }
-                while (next != null && next.marker().pos() == m.pos()
+                while (next != null && next.marker().chromIndex()==m.chromIndex()
+                        && next.marker().pos()==m.pos()
                         && next.marker().equals(m)==false) {
                     next = it.hasNext() ? it.next() : null;
                 }

vcf/TargetData.java

@@ -88,8 +88,7 @@ public class TargetData implements Data {
 
     private static GT targGT(Samples samples, GTRec[] targData,
             Pedigree ped) {
-        GT gt = new BasicGT(samples, targData);
-        return gt.isGTData() ? new XBasicGT(gt, ped) : gt;
+        return new BasicGT(samples, targData);
     }
 
     @Override

vcf/VcfWriter.java

@@ -48,7 +48,7 @@ public final class VcfWriter {
     private static final String gtFormat = "##FORMAT=<ID=GT,Number=1,Type=String,"
             + "Description=\"Genotype\">";
     private static final String dsFormat = "##FORMAT=<ID=DS,Number=A,Type=Float,"
-            +"Description=\"estimated ALT dose [P(RA) + P(AA)]\">";
+            +"Description=\"estimated ALT dose [P(RA) + 2*P(AA)]\">";
     private static final String ap1Format = "##FORMAT=<ID=AP1,Number=A,Type=Float,"
             +"Description=\"estimated ALT dose on first haplotype\">";
     private static final String ap2Format = "##FORMAT=<ID=AP2,Number=A,Type=Float,"

vcf/XBasicGT.java

-/*
- * Copyright (C) 2014-2016 Brian L. Browning
- *
- * This file is part of Beagle
- *
- * Beagle is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Beagle is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-package vcf;
-
-import beagleutil.Samples;
-import blbutil.Const;
-import main.Pedigree;
-
-/**
- * <p>Class {@code XBasicGT} represents genotype and genotype emission
- * probabilities for a set of samples optimized by sample.
- * </p>
- * Instances of class {@code XBasicGTWindow} are immutable.
- *
- * @author Brian L. Browning {@code <browning@uw.edu>}
- */
-public final class XBasicGT implements GT {
-
-    private final Samples samples;
-    private final Markers markers;
-    private final boolean isRefData;
-    private final XBasicGT1[] gt1;
-
-    /**
-     * Constructs a {@code XBasicGT} instance from the specified data.
-     *
-     * @param gl the genotype likelihoods
-     * @param ped the pedigrees
-     *
-     * @throws IllegalArgumentException if
-     * {@code gl.samples().equals(ped.samples())==false}
-     * @throws NullPointerException if {@code gl == null || ped == null}
-     */
-    public XBasicGT(GT gl, Pedigree ped) {
-        if (gl.samples().equals(ped.samples())==false) {
-            throw new IllegalArgumentException("inconsistent samples");
-        }
-        int nSamples = gl.nSamples();
-        this.markers = gl.markers();
-        this.samples = gl.samples();
-        this.isRefData = gl.isPhased();
-        this.gt1 = new XBasicGT1[nSamples];
-        for (int s=0; s<nSamples; ++s) {
-            int father = ped.father(s);
-            int mother = ped.mother(s);
-            gt1[s] = new XBasicGT1(gl, s, father, mother);
-        }
-    }
-
-    @Override
-    public boolean isPhased() {
-        return isRefData;
-    }
-
-    @Override
-    public boolean isGTData() {
-        return true;
-    }
-
-    @Override
-    public boolean isPhased(int sample) {
-        return gt1[sample].isRefSample();
-    }
-
-    @Override
-    public float gl(int marker, int sample, int allele1, int allele2) {
-        return gt1[sample].gl(marker, allele1, allele2);
-    }
-
-    @Override
-    public boolean isPhased(int marker, int sample) {
-        return gt1[sample].isPhased(marker);
-    }
-
-    @Override
-    public int allele1(int marker, int sample) {
-        return gt1[sample].allele1(marker);
-    }
-
-    @Override
-    public int allele2(int marker, int sample) {
-        return gt1[sample].allele2(marker);
-    }
-
-    @Override
-    public int allele(int marker, int hap) {
-        int sample = hap>>1;
-        if ((hap & 1) == 0) {
-            return allele1(marker, sample);
-        }
-        else {
-            return allele2(marker, sample);
-        }
-    }
-
-    @Override
-    public int nMarkers() {
-        return markers.nMarkers();
-    }
-
-    @Override
-    public Marker marker(int markerIndex) {
-        return markers.marker(markerIndex);
-    }
-
-    @Override
-    public Markers markers() {
-        return markers;
-    }
-
-    @Override
-    public int nHaps() {
-        return 2*samples.nSamples();
-    }
-
-    @Override
-    public int nSamples() {
-        return samples.nSamples();
-    }
-
-    @Override
-    public Samples samples() {
-        return samples;
-    }
-
-    @Override
-    public String toString() {
-        StringBuilder sb = new StringBuilder();
-        sb.append("[XBasicGTWindow: nMarkers=");
-        sb.append(nMarkers());
-        sb.append(" nSamples=");
-        sb.append(nSamples());
-        sb.append(Const.nl);
-        for (int m=0, n=nMarkers(); m<n; ++m) {
-            sb.append(markers.marker(m));
-            sb.append(Const.nl);
-            sb.append(Const.MISSING_DATA_CHAR);     // QUAL
-            sb.append(Const.tab);
-            sb.append("PASS");                      // FILTER
-            sb.append(Const.tab);
-            sb.append(Const.MISSING_DATA_CHAR);     // INFO
-            sb.append(Const.tab);
-            sb.append("GT");                        // FORMAT
-            for (XGT1 s : gt1) {
-                sb.append(Const.tab);
-                sb.append(s.allele1(m));
-                sb.append(s.isPhased(m) ? Const.phasedSep : Const.unphasedSep);
-                sb.append(s.allele2(m));
-            }
-        }
-        sb.append(Const.nl);
-        sb.append(']');
-        return sb.toString();
-    }
-}

vcf/XBasicGT1.java

-/*
- * Copyright (C) 2014-2016 Brian L. Browning
- *
- * This file is part of Beagle
- *
- * Beagle is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Beagle is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-package vcf;
-
-import blbutil.Const;
-import java.util.BitSet;
-
-/**
- * <p>Class {@code XBasicGT1} represents genotype likelihoods for one sample.
- * </p>
- * <p>Instances of class {@code XBasicGT1} are immutable.
- * </p>
- *
- * @author Brian L. Browning {@code <browning@uw.edu>}
- */
-public final class XBasicGT1 implements XGT1 {
-
-    private final int idIndex;
-    private final Markers markers;
-    private final BitSet missing1;
-    private final BitSet missing2;
-    private final BitSet allele1;
-    private final BitSet allele2;
-    private final BitSet isPhased;
-    private final boolean isRefSample;
-
-    /**
-     * Constructs a {@code XBasicGL1} instance from the specified data.
-     *
-     * @param gl the genotype likelihoods
-     * @param sample the sample index
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code samples < 0 || sample >= gl.nSamples()}
-     * @throws NullPointerException if {@code gl == null}
-     */
-    public XBasicGT1(GT gl, int sample) {
-        this(gl, sample, -1, -1);
-    }
-
-    /**
-     * Constructs a {@code XBasicGL1} instance from the specified data.
-     *
-     * @param gt the genotype data
-     * @param sample the sample index
-     * @param father the sample index of the sample's father, or -1 if the
-     * father is not genotyped
-     * @param mother the sample index of the sample's mother, or -1 if the
-     * mother is not genotyped
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code samples < 0 || sample >= gl.nSamples()}
-     * @throws NullPointerException if {@code gl == null}
-     */
-    public XBasicGT1(GT gt, int sample, int father, int mother) {
-        // NB: phase information in GL is ignored if parental data is used
-        int nMarkers = gt.markers().nMarkers();
-        int sumHapBits = gt.markers().sumHaplotypeBits();
-        this.idIndex = gt.samples().idIndex(sample);
-        this.markers = gt.markers();
-        this.missing1 = new BitSet(nMarkers);
-        this.missing2 = new BitSet(nMarkers);
-        this.allele1 = new BitSet(sumHapBits);
-        this.allele2 = new BitSet(sumHapBits);
-        this.isPhased = new BitSet(nMarkers);
-//      // Removing use of pedigree constraints to preserve regression
-//      //   after changing main.Par.ped() to return pedigree
-//        if (father==-1 && mother==-1) {
-//            this.isRefSample = copyGL(gl, sample);
-//        }
-//        else {
-//            this.isRefSample = useConstraintAndGL(gl, sample, father, mother);
-//        }
-        this.isRefSample = copyGL(gt, sample);
-    }
-
-    private boolean copyGL(GT gl, int sample) {
-        int nMarkers = markers.nMarkers();
-        for (int m=0; m<nMarkers; ++m) {
-            int a1 = gl.allele1(m, sample);
-            int a2 = gl.allele2(m, sample);
-            setBits(markers, m, a1, allele1, missing1);
-            setBits(markers, m, a2, allele2, missing2);
-            if (gl.isPhased(m, sample)) {
-                isPhased.set(m);
-            }
-        }
-        return gl.isPhased(sample);
-    }
-
-    private boolean useConstraintAndGL(GT gl, int sample, int father,
-            int mother) {
-        int nMarkers = markers.nMarkers();
-        int phasedCnt = 0;
-        for (int m=0; m<nMarkers; ++m) {
-            int a1 = gl.allele1(m, sample);
-            int a2 = gl.allele2(m, sample);
-            int tr1 = transmittedAllele(gl, m, father);
-            int tr2 = transmittedAllele(gl, m, mother);
-            if ((tr1>=0 || tr2>=0) && isUnphasedConsistent(a1, a2, tr1, tr2)) {
-                if (isPhasedConsistent(a1, a2, tr1, tr2)==false) {
-                    int tmp = a1;
-                    a1 = a2;
-                    a2 = tmp;
-                }
-                if (tr1 != -1) {
-                    a1 = tr1;
-                }
-                if (tr2 != -1) {
-                    a2 = tr2;
-                }
-                if (a1>=0 && a2>=0) {
-                    ++phasedCnt;
-                    isPhased.set(m);
-                }
-            }
-            setBits(markers, m, a1, allele1, missing1);
-            setBits(markers, m, a2, allele2, missing2);
-        }
-        return phasedCnt==nMarkers;
-}
-
-    private static boolean isUnphasedConsistent(int a1, int a2, int b1, int b2) {
-        return isPhasedConsistent(a1, a2, b1, b2)
-                || isPhasedConsistent(a1, a2, b2, b1);
-    }
-
-    private static boolean isPhasedConsistent(int a1, int a2, int b1, int b2) {
-        return (isConsistent(a1, b1) && isConsistent(a2, b2));
-    }
-
-    private static boolean isConsistent(int a1, int b1) {
-        return (a1==-1 || b1==-1 || a1==b1);
-    }
-
-    private static int transmittedAllele(GT gl, int marker, int sample) {
-        if (sample==-1) {
-            return -1;
-        }
-        int a1 = gl.allele1(marker, sample);
-        return a1>=0 && a1==gl.allele2(marker, sample) ? a1 : -1;
-    }
-
-    private static void setBits(Markers markers, int marker, int allele,
-            BitSet alleles, BitSet missing) {
-        if (allele == -1) {
-            missing.set(marker);
-        }
-        else {
-            int mask = 1;
-            int start = markers.sumHaplotypeBits(marker);
-            int end = markers.sumHaplotypeBits(marker+1);
-            for (int i=start; i<end; ++i) {
-                boolean b = (allele & mask)==mask;
-                alleles.set(i, b);
-                mask <<= 1;
-            }
-        }
-    }
-
-    private int allele(BitSet bitset, int marker) {
-        int start = markers.sumHaplotypeBits(marker);
-        int end = markers.sumHaplotypeBits(marker+1);
-        if (end==(start+1)) {
-            return bitset.get(start) ? 1 : 0;
-        }
-        int allele = 0;
-        int mask = 1;
-        for (int j=start; j<end; ++j) {
-            if (bitset.get(j)) {
-                allele += mask;
-            }
-            mask <<= 1;
-        }
-        return allele;
-    }
-
-    @Override
-    public boolean isRefSample() {
-        return isRefSample;
-    }
-
-    @Override
-    public float gl(int marker, int a1, int a2) {
-        if (a1 < 0 || a1 >= markers.marker(marker).nAlleles()) {
-            throw new IndexOutOfBoundsException(String.valueOf(a1));
-        }
-        if (a2 < 0 || a2 >= markers.marker(marker).nAlleles()) {
-            throw new IndexOutOfBoundsException(String.valueOf(a2));
-        }
-        int obsA1 = allele1(marker);
-        int obsA2 = allele2(marker);
-        boolean consistent = (obsA1==-1 || obsA1==a1) && (obsA2==-1 || obsA2==a2);
-        if (consistent==false && isPhased.get(marker)==false) {
-            consistent = (obsA1==-1 || obsA1==a2) && (obsA2==-1 || obsA2==a1);
-        }
-        return consistent ? 1.0f : 0.0f;
-    }
-
-    @Override
-    public boolean isPhased(int marker) {
-        return isPhased.get(marker);
-    }
-
-    @Override
-    public int allele1(int marker) {
-        return missing1.get(marker) ? -1 : allele(allele1, marker);
-    }
-
-    @Override
-    public int allele2(int marker) {
-        return missing2.get(marker) ? -1 : allele(allele2, marker);
-    }
-
-    @Override
-    public int nMarkers() {
-        return markers.nMarkers();
-    }
-
-    @Override
-    public Marker marker(int markerIndex) {
-        return markers.marker(markerIndex);
-    }
-
-    @Override
-    public Markers markers() {
-        return markers;
-    }
-
-    @Override
-    public int idIndex() {
-        return idIndex;
-    }
-
-    @Override
-    public String toString() {
-        int nMarkers = markers.nMarkers();
-        StringBuilder sb  = new StringBuilder();
-        sb.append("[XBasicGL1: nMarkers=");
-        sb.append(nMarkers);
-        sb.append(Const.nl);
-        for (int m=0; m<nMarkers; ++m) {
-            sb.append(markers.marker(m));
-            sb.append(Const.tab);
-            sb.append(allele1(m));
-            sb.append(isPhased(m) ? Const.phasedSep : Const.unphasedSep);
-            sb.append(allele2(m));
-            sb.append(Const.nl);
-        }
-        sb.append(']');
-        return sb.toString();
-    }
-}

vcf/XGT1.java

-/*
- * Copyright (C) 2014-2016 Brian L. Browning
- *
- * This file is part of Beagle
- *
- * Beagle is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Beagle is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-package vcf;
-
-// Later consider relaxing mutability, and determine whether interface
-// can be replaced with implementing class
-
-/**
- * <p>Interface {@code XGT1} (Genotype Likelihoods) represents genotype
- * likelihoods for one sample.
- * </p>
- * <p>Instances of {@code XGT1} are required to be immutable.
- * </p>
- *
- * @author Brian L. Browning {@code <browning@uw.edu>}
- */
-public interface XGT1 {
-
-    /**
-     * Returns {@code true} if the observed data for each marker
-     * includes a phased genotype that has no missing alleles,
-     * and returns {@code false} otherwise.
-     * @return {@code true} if the observed data for each marker
-     * includes a phased genotype that has no missing alleles,
-     * and {@code false} otherwise
-     */
-    boolean isRefSample();
-
-    /**
-     * Returns the probability of the observed data for the specified marker
-     * if the specified pair of ordered alleles is the true ordered genotype.
-     * @param marker the marker index
-     * @param allele1 the first allele index
-     * @param allele2 the second allele index
-     * @return the probability of the observed data for the specified marker
-     * and sample if the specified pair of ordered alleles is the true
-     * ordered genotype
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     * @throws IndexOutOfBoundsException if
-     * {@code allele1 < 0 || allele1 >= this.marker(marker).nAlleles()}
-     * @throws IndexOutOfBoundsException if
-     * {@code allele2 < 0 || allele2 >= this.marker(marker).nAlleles()}
-     */
-    float gl(int marker, int allele1, int allele2);
-
-    /**
-     * Returns {@code true} if the observed data for the specified
-     * marker includes a phased genotype, and returns {@code false} otherwise.
-     * @param marker the marker index
-     * @return {@code true} if the observed data for the specified
-     * marker includes a phased genotype, and {@code false} otherwise
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     */
-    boolean isPhased(int marker);
-
-    /**
-     * Returns the first allele for the specified marker if the observed data
-     * include a non-missing allele, and returns -1 otherwise. Alleles are
-     * arbitrarily ordered if the genotype is unphased.
-     * @param marker the marker index
-     * @return the first allele for the specified marker if the observed data
-     * include a non-missing allele, and -1 otherwise
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     */
-    int allele1(int marker);
-
-    /**
-     * Returns the second allele for the specified marker if the observed data
-     * include a non-missing allele, and returns -1 otherwise. Alleles are
-     * arbitrarily ordered if the genotype is unphased.
-     * @param marker the marker index
-     * @return the second allele for the specified marker if the observed data
-     * include a non-missing allele, and -1 otherwise
-     *
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     */
-    int allele2(int marker);
-
-    /**
-     * Returns the number of markers.
-     * @return the number of markers
-     */
-    int nMarkers();
-
-    /**
-     * Returns the specified marker.
-     * @param marker the marker index
-     * @return the specified marker
-     * @throws IndexOutOfBoundsException if
-     * {@code marker < 0 || marker >= this.nMarkers()}
-     */
-    Marker marker(int marker);
-
-    /**
-     * Returns the list of markers.
-     * @return the list of markers
-     */
-    Markers markers();
-
-    /**
-     * Returns the sample identifier index.
-     * @return the sample identifier index
-     */
-    int idIndex();
-
-    /**
-     * Returns a string representation of {@code this}.  The exact
-     * details of the representation are unspecified and subject to change.
-     * @return a string representation of {@code this}
-     */
-    @Override
-    String toString();
-}