Phasing using Beagle with a map file

Question

[michael.flower.14]

I'd like to phase the SNPs in a vcf file and output consensus files for each haplotype, as suggested in this post:

https://www.biostars.org/p/298635/

I've managed to install beagle in a conda environment:

conda create -n beagle -c conda-forge -c bioconda beagle
conda activate beagle

When I run beagle using this command:

java -Xmx2g -jar "/Users/michaelflower/opt/anaconda3/envs/beagle/share/beagle-5.2_21Apr21.304-0/beagle.jar" gt="$VCF" out="$DIR"/beagle/phased_beagle

I get the error message:

No genetic map is specified: using 1 cM = 1 Mb

The full output of beagle is below. Also the vcf file produced contains no SNPs.

So as suggested here, I downloaded the plink.GRCh37.map.zip map file from here:

wget -qO- http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh37.map.zip | bsdtar -xvf- -C ~/refs/hap/hg19/

The SNPs I'm interested in phasing are on chr4, so I set the map shortcut as follows:

MAP=~/refs/hap/hg19/plink.chr4.GRCh37.map

But when I try and run beagle like this:

java -Xmx2g -jar "/Users/michaelflower/opt/anaconda3/envs/beagle/share/beagle-5.2_21Apr21.304-0/beagle.jar" gt="$VCF" map="$MAP" out="$DIR"/beagle/phased_beagle

I get this error (see full error 2 below):

java.lang.IllegalArgumentException: missing genetic map for chromosome chr4

Ful beagle error 1:

$ java -Xmx2g -jar "/Users/michaelflower/opt/anaconda3/envs/beagle/share/beagle-5.2_21Apr21.304-0/beagle.jar" gt="$VCF" out="$DIR"/beagle/phased_beagle
beagle.21Apr21.304.jar (version 5.2)
Copyright (C) 2014-2021 Brian L. Browning
Enter "java -jar beagle.21Apr21.304.jar" to list command line argument
Start time: 04:44 PM GMT on 21 Nov 2021

Command line: java -Xmx2048m -jar beagle.21Apr21.304.jar
  gt=/Users/michaelflower/Documents/ACL/Research/Projects/HTT ASO/2021.11.20 125Q SNP phasing/gatk/125QiPSC_chr4_HaplotypeCallerPGT.vcf
  out=/Users/michaelflower/Documents/ACL/Research/Projects/HTT ASO/2021.11.20 125Q SNP phasing/beagle/phased_beagle

No genetic map is specified: using 1 cM = 1 Mb

Reference samples:                    0
Study     samples:                    1

Window 1 [chr4:1-39999997]
Reference markers:            2,463,043
Study     markers:            2,463,043
Exception in thread "main" java.lang.IllegalArgumentException: java.lang.IllegalArgumentException: bound must be positive
    at java.base/jdk.internal.reflect.NativeConstructorAccessorImpl.newInstance0(Native Method)
    at java.base/jdk.internal.reflect.NativeConstructorAccessorImpl.newInstance(NativeConstructorAccessorImpl.java:62)
    at java.base/jdk.internal.reflect.DelegatingConstructorAccessorImpl.newInstance(DelegatingConstructorAccessorImpl.java:45)
    at java.base/java.lang.reflect.Constructor.newInstance(Constructor.java:490)
    at java.base/java.util.concurrent.ForkJoinTask.getThrowableException(ForkJoinTask.java:600)
    at java.base/java.util.concurrent.ForkJoinTask.reportException(ForkJoinTask.java:678)
    at java.base/java.util.concurrent.ForkJoinTask.invoke(ForkJoinTask.java:737)
    at java.base/java.util.stream.Nodes.collect(Nodes.java:333)
    at java.base/java.util.stream.ReferencePipeline.evaluateToNode(ReferencePipeline.java:109)
    at java.base/java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:545)
    at java.base/java.util.stream.AbstractPipeline.evaluateToArrayNode(AbstractPipeline.java:260)
    at java.base/java.util.stream.ReferencePipeline.toArray(ReferencePipeline.java:517)
    at phase.PbwtPhaser.pbwtPhasers(PbwtPhaser.java:183)
    at phase.PbwtPhaser.initPhase(PbwtPhaser.java:72)
    at phase.EstPhase.<init>(EstPhase.java:46)
    at phase.PhaseData.<init>(PhaseData.java:52)
    at main.Main.phaseStage1Variants(Main.java:188)
    at main.Main.phaseTarg(Main.java:181)
    at main.Main.phaseAndImpute(Main.java:171)
    at main.Main.main(Main.java:126)
Caused by: java.lang.IllegalArgumentException: bound must be positive
    at java.base/java.util.Random.nextInt(Random.java:388)
    at phase.RevPbwtPhaser.imputeAllele(RevPbwtPhaser.java:128)
    at phase.RevPbwtPhaser.finishPhasing(RevPbwtPhaser.java:118)
    at phase.RevPbwtPhaser.phase(RevPbwtPhaser.java:94)
    at phase.RevPbwtPhaser.<init>(RevPbwtPhaser.java:75)
    at phase.FwdPbwtPhaser.phase(FwdPbwtPhaser.java:86)
    at phase.FwdPbwtPhaser.<init>(FwdPbwtPhaser.java:79)
    at phase.PbwtPhaser.<init>(PbwtPhaser.java:59)
    at phase.PbwtPhaser.lambda$pbwtPhasers$1(PbwtPhaser.java:181)
    at java.base/java.util.stream.IntPipeline$1$1.accept(IntPipeline.java:180)
    at java.base/java.util.stream.Streams$RangeIntSpliterator.forEachRemaining(Streams.java:104)
    at java.base/java.util.Spliterator$OfInt.forEachRemaining(Spliterator.java:699)
    at java.base/java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:484)
    at java.base/java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:474)
    at java.base/java.util.stream.Nodes$SizedCollectorTask.compute(Nodes.java:1886)
    at java.base/java.util.concurrent.CountedCompleter.exec(CountedCompleter.java:746)
    at java.base/java.util.concurrent.ForkJoinTask.doExec(ForkJoinTask.java:290)
    at java.base/java.util.concurrent.ForkJoinPool$WorkQueue.topLevelExec(ForkJoinPool.java:1020)
    at java.base/java.util.concurrent.ForkJoinPool.scan(ForkJoinPool.java:1656)
    at java.base/java.util.concurrent.ForkJoinPool.runWorker(ForkJoinPool.java:1594)
    at java.base/java.util.concurrent.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:183)

Full beagle error 2:

$ java -Xmx2g -jar "/Users/michaelflower/opt/anaconda3/envs/beagle/share/beagle-5.2_21Apr21.304-0/beagle.jar" gt="$VCF" map="$MAP" out="$DIR"/beagle/phased_beagle
beagle.21Apr21.304.jar (version 5.2)
Copyright (C) 2014-2021 Brian L. Browning
Enter "java -jar beagle.21Apr21.304.jar" to list command line argument
Start time: 04:49 PM GMT on 21 Nov 2021

Command line: java -Xmx2048m -jar beagle.21Apr21.304.jar
  gt=/Users/michaelflower/Documents/ACL/Research/Projects/HTT ASO/2021.11.20 125Q SNP phasing/gatk/125QiPSC_chr4_HaplotypeCallerPGT.vcf
  map=/Users/michaelflower/refs/hap/hg19/plink.chr4.GRCh37.map
  out=/Users/michaelflower/Documents/ACL/Research/Projects/HTT ASO/2021.11.20 125Q SNP phasing/beagle/phased_beagle
java.lang.IllegalArgumentException: missing genetic map for chromosome chr4
    at vcf.PlinkGenMap.checkChromIndex(PlinkGenMap.java:213)
    at vcf.PlinkGenMap.genPos(PlinkGenMap.java:311)
    at vcf.PlinkGenMap.genPos(PlinkGenMap.java:306)
    at vcf.WindowIt$Reader.run(WindowIt.java:255)
    at java.base/java.lang.Thread.run(Thread.java:834)

Answer 1

[Guangyou]

I encountered the same error. The point is the prefix "chr". It's "4" in map but "chr4" in your vcf. So in the error message:

java.lang.IllegalArgumentException: missing genetic map for chromosome chr4

It actually wants to say it should be "4" instead of "chr4"

One simple solution is adding "chr" for each line of map file, like this:

sed 's/^/chr/' ~/refs/hap/hg19/plink.chr4.GRCh37.map > ~/refs/hap/hg19/plink.chr4.GRCh37.map.modified

And use that plink.chr4.GRCh37.map.modified as the map file.

Answer 2

[chaselizgr]

No genetic map is specified: using 1 cM = 1 Mb，Window and overlapping area values should be scaled down in cM.