BSMAP_methdiff.py problem

[Sunwoo Lee]

Hi, everybody,
I'm doing data analysis with my WGBS data with BSMAP.

I finished BSMAP(2.90)-alignment, methratio.py.
But when I run methdiff.py, I got the error as below;

* methratio_1 (sample 1)
chr pos strand context ratio eff_CT_count C_count CT_count rev_G_count rev_GA_count CI_lower CI_upper
CHR1 360 - CHH 0 3 0 3 8 8 0 0.562
CHR1 363 + CHH 0 8 0 8 4 4 0 0.324
CHR1 365 + CHH 0 6.4 0 8 4 5 0 0.375
CHR1 367 + CHH 0.15 6.67 1 8 5 6 0.027 0.529
CHR1 370 + CHH 0 6 0 8 6 8 0 0.39
CHR1 371 + CHH 0 6.67 0 8 5 6 0 0.366
CHR1 373 + CHG 0 6.75 0 9 6 8 0 0.363
CHR1 374 + CG 0.508 7.88 4 9 7 8 0.219 0.792
CHR1 375 - CG 0.6 3.33 2 6 5 9 0.184 0.909
CHR1 376 + CHH 0 6.3 0 9 7 10 0 0.379
CHR1 377 + CHH 0 9 0 9 13 13 0 0.299
CHR1 380 + CHG 0 9 0 9 16 16 0 0.299
CHR1 381 + CG 1 9 9 9 17 17 0.701 1

* methratio_4 (sample 2)

chr pos strand context ratio eff_CT_count C_count CT_count rev_G_count rev_GA_count CI_lower CI_upper
CHR1 360 - CHH 0.25 4 1 4 7 7 0.046 0.699
CHR1 363 + CHH 0 5.83 0 7 5 6 0 0.397
CHR1 365 + CHH 0 6 0 7 6 7 0 0.39
CHR1 367 + CHH 0 5 0 7 5 7 0 0.434
CHR1 370 + CHH 0 5.25 0 7 6 8 0 0.423
CHR1 371 + CHH 0 6 0 7 6 7 0 0.39
CHR1 373 + CHG 0 6.22 0 7 8 9 0 0.382
CHR1 374 + CG 1 5.83 7 7 10 12 0.603 1
CHR1 375 - CG 1 10 10 10 7 7 0.722 1
CHR1 376 + CHH 0 5.6 0 7 12 15 0 0.407
CHR1 377 + CHH 0 6.18 0 7 15 17 0 0.383
CHR1 380 + CHG 0 6.59 0 7 16 17 0 0.368
CHR1 381 + CG 1 6.61 7 7 17 18 0.632 1

.....When I run methdiff.py as below~
$ python methdiff.py -o /home/lsw/NGS_result_data/PNU_flatfish/Bisul-Seq/methdiff_out_b1 -d /home/lsw/NGS_result_data/PNU_flatfish/genome.fasta -b 1 -p 1 -r 0 /home/lsw/NGS_result_data/PNU_flatfish/Bisul-Seq/methratio_1 /home/lsw/NGS_result_data/PNU_flatfish/Bisul-Seq/methratio_4
The result is..... (
cr start end p_value diff_ratio group0_ratio group0_depth group0_meth group1_ratio group1_depth group1_meth
CHR1 360 360 1 -0.25 0 3 0 0.25 4 1
CHR1 363 363 1 0 0 8 0 0 5.8 0
CHR1 365 365 1 0 0 6.4 0 0 6 0
CHR1 367 367 1 0.15 0.15 6.7 1 0 5 0
CHR1 370 370 1 0 0 6 0 0 5.2 0
CHR1 371 371 1 0 0 6.7 0 0 6 0
CHR1 373 373 1 0 0 6.8 0 0 6.2 0

It means, there is some error in "CHR1 374".
Also, further test shows there is some error in "CHR1 381".
with the error below~

Traceback (most recent call last):
File "methdiff.py", line 133, in <module>
cmp_chrom(cr)
File "methdiff.py", line 112, in cmp_chrom
pval = get_pval(m[0], d[0], m[1], d[1])
File "methdiff.py", line 87, in get_pval
l1, u1 = conf_intv(m1, d1, z0)
File "methdiff.py", line 81, in conf_intv
span = z * (p * (1 - p) / d + z2 / (4 * d * d)) ** 0.5
ValueError: negative number cannot be raised to a fractional power
When I remove the chr1 374 in methratio result, there was no error.
so was when I remove the chr1 381.

Could you tell me some clue or comment for me?

Thank you very much, everyone.