search

stschiff / msmc

Implementation of the multiple sequential markovian coalescent
GNU General Public License v3.0
87 stars 20 forks source link

AssertError bootstrap #34

Closed bioinformaticspcj closed 6 years ago

bioinformaticspcj commented 6 years ago

Hi, I just try to run msmc bootstrap to analyze population separations, but one error occured as follows:

maxIterations: 20 mutationRate: 0.00121796 recombinationRate: 0.000304489 subpopLabels: [0, 0, 0, 0, 1, 1, 1, 1] timeSegmentPattern: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] nrThreads: 96 nrTtotSegments: 40 verbose: false outFilePrefix: /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out naiveImplementation: false hmmStrideWidth: 1000 fixedPopSize: false fixedRecombination: true initialLambdaVec: [] directedEmissions: false skipAmbiguous: true indices: [0, 1, 2, 3, 4, 5, 6, 7] logging information written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.log loop information written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.loop.txt final results written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.final.txt [11/30] estimating total branchlengths [1/20] Baumwelch iteration

The command I used was: /data/project/tools/msmc-tools-master/msmc_1.0.0_linux64bit --fixedRecombination --skipAmbiguous -P 0,0,0,0,1,1,1,1 -o /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.bootrap.out_1/*.txt

I do not know why this would happen. I am looking forward to getting your help.

stschiff commented 6 years ago

Could be overfitting. Have you tried just running on 4 haplotypes first? Also, I have a new version in the dev branch of the github repository which implements an option to constrain the coalescence rates between a minimum and a maximum, which might help avoiding under- and overflow exceptions.

Stephan

On 10 Apr 2018, at 04:16, bioinformaticspcj notifications@github.com wrote:

Hi, I just try to run msmc bootstrap to analyze population separations, but one error occured as follows:

maxIterations: 20 mutationRate: 0.00121796 recombinationRate: 0.000304489 subpopLabels: [0, 0, 0, 0, 1, 1, 1, 1] timeSegmentPattern: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] nrThreads: 96 nrTtotSegments: 40 verbose: false outFilePrefix: /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out naiveImplementation: false hmmStrideWidth: 1000 fixedPopSize: false fixedRecombination: true initialLambdaVec: [] directedEmissions: false skipAmbiguous: true indices: [0, 1, 2, 3, 4, 5, 6, 7] logging information written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.log loop information written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.loop.txt final results written to /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out.final.txt [11/30] estimating total branchlengths [1/20] Baumwelch iteration

[30/30] Expectation Step, log likelihood: -1643.79 [4/200(max)] Maximization Step, Q-function before: 220.606, after:220.403 [2/20] Baumwelch iteration [30/30] Expectation Step, log likelihood: -923.319 [2/200(max)] Maximization Step, Q-function before: 2.65402, after:2.6227 [3/20] Baumwelch iteration [29/30] Expectation Step, log likelihood: -865.246 [2/200(max)] Maximization Step, Q-function before: 1.04843, after:1.02013 [4/20] Baumwelch iteration [30/30] Expectation Step, log likelihood: -818.439 [2/200(max)] Maximization Step, Q-function before: 0.360145, after:0.344497 [5/20] Baumwelch iteration [16/30] Expectation Step, log likelihood: -794.868 [3/200(max)] Maximization Step, Q-function before: 0.10929, after:0.103583 [6/20] Baumwelch iteration [30/30] Expectation Step, log likelihood: -786.659 [85/200(max)] Maximization Step, Q-function before: 0.0316963, after:0.0312611 [7/20] Baumwelch iteration [27/30] Expectation Step, log likelihood: -785.941 [33/200(max)] Maximization Stepcore.exception.AssertError@model/time_intervals.d(147): Assertion failure ??:? _d_assert [0x5ebe93] ??:? void model.time_intervals.assert(int) [0x58c30c] ??:? _D5model14time_intervals13TimeIntervals18meanTimeWithLambdaMxFmdZ8ensureMFKxdZv [0x58bf25] ??:? const(double function(ulong, double)) model.time_intervals.TimeIntervals.meanTimeWithLambda [0x58be93] ??:? const(double function(double, ulong, ulong, double, ulong, ulong)) model.transition_rate.TransitionRate.q2IntegralSmaller [0x58cdc5] ??:? const(double function(ulong, ulong)) model.transition_rate.TransitionRate.transitionProbabilityOffDiagonal [0x58cad0] ??:? _D5model15transition_rate14TransitionRate35fillTransitionProbabilitiesParallelMFZ14__foreachbody1MFmZi [0x58c8b2] ??:? _D3std11parallelism64T15ParallelForeachTS3std5range13T4iotaTmTmZ4iotaFmmZ6ResultZ15ParallelForeach7opApplyMFMDFmZiZ4doItMFZv [0x58e39a] ??:? void std.parallelism.run!(void delegate()).run(void delegate()) [0x5fb49b] ??:? void std.parallelism.__T4TaskS213std11parallelism3runTDFZvZ.Task.impl(void) [0x5faf9b] ??:? void std.parallelism.AbstractTask.job() [0x62933e] ??:? void std.parallelism.TaskPool.doJob(std.parallelism.AbstractTask) [0x5f9c37] ??:? void std.parallelism.TaskPool.executeWorkLoop() [0x5f9d65] ??:? void std.parallelism.TaskPool.startWorkLoop() [0x5f9d0d] ??:? void core.thread.Thread.run() [0x630481] ??:? thread_entryPoint [0x6300f1] ??:? [0x6c406d73] core.exception.AssertError@model/time_intervals.d(147): Assertion failure

The command I used was: /data/project/tools/msmc-tools-master/msmc_1.0.0_linux64bit --fixedRecombination --skipAmbiguous -P 0,0,0,0,1,1,1,1 -o /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.out/hmh.boot.1.out /data/project/pcj/primates/Macaca_arctoides/msmc/msmc1/YND_YNQ/all.msmc.bootrap.out_1/*.txt

I do not know why this would happen. I am looking forward to getting your help.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/stschiff/msmc/issues/34, or mute the thread https://github.com/notifications/unsubscribe-auth/AAbQmufBPmJJV7G2UJX1HQQ5oUq7ZCMxks5tnBYHgaJpZM4TNfyX.

stschiff commented 6 years ago

Note that I have now built in an option to artificially constrain the coalescence rates, see release 1.1.0, with options --loBoundLambda and --hiBoundLambda. Lambda is the scaled coalescence rate, with typical values around 1000 or so (it's on average around 2/theta, where theta is the heterozygosity), so you could try setting a lower bound of 10, and an upper bound of 1e6 or so, which might solve this issue.