在使用homer 做motif 分析时,报如下错误
Illegal division by zero at /data/zhangyu/bin/findKnownMotifs.pl line 154.
然后 knownResults 的结果为空,denovo的结果正常。
经过排查,发现是contig 问题。
Chr01
Chr02
Chr03
Chr04
Chr05
Chr06
Chr07
Chr08
Chr09
Chr10
Chr11
Chr12
Contig07
Contig14
Contig30
Contig31
去掉contig 的部分,即可正常运行。
具体原始不是很清楚,单独做每一个contig的分析也是OK的,所有contig 的放在一起做也是OK的,但是Chr 和contig 的放在一起做就有问题,不清楚是什么原因。
R语言版本的说是染色体字符的问题,但是perl 版本的测试了不是这个问题。
https://github.com/r3fang/SnapATAC/issues/202
2024.3.13
最后测试发现程序调用了 findKnownMotifs.pl,这一步出的问题。
findKnownMotifs.pl -s homer/seq.tsv -g homer/group.adj -o homer -pvalue 0.01 -m /data/zhangyu/miniconda3/share/homer/.//data/knownTFs/all/known.motifs -homer2 -p 24 -stat binomial -cache 500
这个命令没有什么问题,但是产生的结果有问题,第一个数据的最后一个值,背景的百分比是0.00%,对于我的数据,猜测原因是背景序列中有很多被mask 掉了。
Motif Name Consensus p-value Log p-value Log2 enrichment ratio # target sequences with Motif(total=28232) % of target sequences with Motif # background sequences with Motif(total=27111) % of background sequences with Motif
AT4G12670(MYBrelated)/col-AT4G12670-DAP-Seq(GSE60143)/Homer AGGGTTTAGGGTTTA 1e-315 -726.818 7.393 175 0.62% 0 0.00%
ETV4(ETS)/HepG2-ETV4-ChIP-Seq(ENCODE)/Homer ACCGGAAGTG 1e-51 -119.008 0.747 991 3.51% 567 2.09%
SCL(bHLH)/HPC7-Scl-ChIP-Seq(GSE13511)/Homer ANCAGCTG 1e-46 -108.105 0.268 5347 18.94% 4263 15.73%
……
程序中对这一步的判断是直接对最后一个值去掉百分号,再*0.01,所以这里是0的话,就会报错:
151 $percentTargets =~ s/\%//;
152 $percentBackground =~ s/\%//;
153 $totalNumTargets = floor($numTargets/($percentTargets*0.01)+0.5);
154 $totalNumBackground = floor($numBackground/($percentBackground*0.01)+0.5);
如果 $percentBackground == 0,这里可以做两种处理:
1. 直接跳过这个motif的分析,分析下一个motif。
2. 将$percentBackground = 0.01 ,给一个很小的值。
这里我选第2种,修改后如下:
151 $percentTargets =~ s/\%//;
152 $percentBackground =~ s/\%//;
153 $percentBackground = 0.01 if $percentBackground == 0;
154 $totalNumTargets = floor($numTargets/($percentTargets*0.01)+0.5);
155 $totalNumBackground = floor($numBackground/($percentBackground*0.01)+0.5)
结果正常运行。