sergpolly
commented
6 years ago An entire Snakefile that I'd use to call dots at 5 and 10kb and combine for a bunch of samples/coolers.
Broken down and commented:
you'd typically need ~10kb and 5kb coolers to perform GM-level dot-calling also, describe your samples here - i.e. ./path(1,2)/your-sample(1,2).cool must exist
resolutions = [5000,10000]
names = ["your-sample1","your-sample2"]
paths = ["path1", "path2"]the snake-thing about wildcards ...
wildcard_constraints:
res="\d+"get the combined bedpe lists with the dots
rule all:
input:
expand("{path}/combineddots/final_cloops_{hic_name}.combined.bedpe",zip,path=paths,hic_name=names)force rebalance - in case you'd want cis only, modify tolerance, diags, whatever ...
rule balance:
input:
"{path}/{hic_name}.{res}.cool"
output:
touch("{path}/touchtmp/{hic_name}.{res}.touch")
threads: 9
shell:
"cooler balance -p {threads} --ignore-diags 1 --force {input}"expected ...
rule compute_expected:
input:
"{path}/{hic_name}.{res}.cool",
"{path}/touchtmp/{hic_name}.{res}.touch"
output:
"{path}/expected/{hic_name}.{res}.cool.cis.expected.ALL"
threads: 9
shell:
"cooltools compute_expected -p {threads} --drop-diags 1 {input[0]} > {output}"chrM and sometimes chrY screw up dot calling - remove em from expected so that dot-caller won't use them
rule clean_expected:
input:
"{path}/expected/{hic_name}.{res}.cool.cis.expected.ALL"
output:
"{path}/expected/{hic_name}.{res}.cool.cis.expected"
shell:
"grep -v -e \"^chrM\" {input} | grep -v -e \"^chrY\" > {output}"dot-calling for each sample/resolution ...
rule call_dots:
input:
cooler = "{path}/{hic_name}.{res}.cool",
expected = "{path}/expected/{hic_name}.{res}.cool.cis.expected"
params:
fdr = 0.1,
diag_width = 10000000,
tile_size = 10000000
output:
signif_dots="{path}/dots/cloops_{hic_name}.{res}.bedpe",
filtered_dots="{path}/dots/final_cloops_{hic_name}.{res}.bedpe"
threads: 17
shell:
"cooltools call_dots -n {threads} "
" -o {output.signif_dots} -v --fdr {params.fdr} "
" --max-nans-tolerated 7 "
" --max-loci-separation {params.diag_width} "
" --dots-clustering-radius 21000 "
" --tile-size {params.tile_size} "
" {input.cooler} {input.expected}"merge dots using little unfinished script that should probably end up in cooltools at some point https://github.com/sergpolly/random_tools
rule merge_dots_res:
input:
expand("{{path}}/dots/final_cloops_{{hic_name}}.{res}.bedpe",res=resolutions)
params:
radius = 10000
output:
"{path}/combineddots/final_cloops_{hic_name}.combined.bedpe"
shell:
"python ./random_tools/merge_dot_lists.py "
" --radius {params.radius} -v "
" --output {output} {input}"in case anyone needs this - to prepare multires bedpe for higlass ingestion ...
# rule prepare_multires:
# input:
# "final_cloops_{hic_name}.combined.bedpe"
# output:
# "final_cloops_{hic_name}.combined.bedpe.multires"
# params:
# genome = "hg19"
# shell:
# "clodius aggregate bedpe "
# " -a {params.genome} --has-header "
# " --chr1-col 1 --from1-col 2 --to1-col 3 "
# " --chr2-col 4 --from2-col 5 --to2-col 6 "
# " --output-file {output} {input} "
Phlya
gfudenberg
xinyangbing
This issue is going to briefly touch upon the state of the dot-calling in
cooltools.Intro
We were trying to re-implement HiCCUPS dot-calling algorithm under the
cooltoolsumbrella for some time now. It is still under active development and right now code is scattered across forks and branches.master
The initial progress that we made with dot-calling, by implementing convolution based calculation of locally adjusted expected (
donut,lowleft,vertical,horizontal) is reflected in this repository in themasterbranch. The post processing steps in thismasterare closest to the so-called BH-FDR version of dot-calling in the original HiCCUPS paper (Rao etal 2014) - in a sense that we do not do the lambda-chunking to perform multiple hypothesis testing. Moreover this implementation simply ends with the dump of the pre-calculated adjusted expected for different kernels:donut,lowleft,vertical,horizontal, and reports that post-processing in a bad shape. Thus this isn't very usable for now, not for the final dot-calling at least.new-dotfinder
Lambda-chunking procedure is implemented in the
dekkerlabfork of thecooltoolsbranchnew-dotfinder, which ispip-installable:The second command would allow to modify the source code, whereas the 1st one would simply install it. One would want to modify the source code if the enrichment threshold modification is needed, for instance - as those are not implemented as CLI options just yet. A typical run would be:
Which would produce list of dots that passed multiple hypothesis testing (the lambda-chunking step itself) but haven't been post-processed, i.e. clustered or filtered by enrichment. The post processed list of dots would show up in the same folder as
signif_dotsbut with the prefixfinal_- we'll fix this ugliness later on of course.call_dotsCLI determines resolution and pick correct kernels parametersw,paccordingly (all of the defaults kept same as HiCCUPS). This dot-caller albeit very close to HiCCUPS implementation, deviates from it in some regards, some of the most importnat aspects:lowleft)Birch, HiCCUPS implements special greedy algorithm for that - results are very close anyways.NaNsafter balancing): HiCCUPS disregard pixels that are within~5pixels away from bad-rows/cols, instead we simply check number ofNaNsin a kernel footprint--max-nans-tolerated- given the resolution/size of the kernels one can realize which pixels would be discarded.... I might add more details here later on, and edit/elaborate more on this here ...
shrink-donut-dotfinder
This branch elaborates further on top of the
new-dotfinderby dealing with the disrepancy#1- dynamic kernels resizing. As it follows from the name here we implemented only the near-diagonal kernels shrinakge - arguably the most important aspect of the dynamic kernels, which was preventing us from calling dots really close to the diagonal and driving the deviation betweencooltoolsdot-calling and HiCCUPS. There are no additional parameters that user needs to control in this case, everything is done the same way as in HiCCUPS - and this is data-independent kernel shrinkage, as opposed to the enlarging kernels based on the value of thelowleftkernels for each pixel tested.This branch eliminates another difference between HiCCUPS and
cooltoolswhich is related to the way lambda-chunking is implemented and is too technical to describe here. To give some numbers, for Rao et al 2014 GM... primary dataset we are getting ~7700 dots vs ~8050 by HiCCUPS, where overlap is ~7600.