cellphy.snake

"""
Explanation of wildcards:

fineTuneTypeCellPhy: types of fine-tuning to perform for CellPhy defined in the configuration file, e.g., single thread, multi thread, etc.
dataType: input data types, e.g., no CNVs, including CNVs, excluding CNVs, etc.
datasetName: simulated dataset names, e.g., 0001, 0002, 0003...
"""


include: "scripts/constants.py"
include: "scripts/utils.py"


#########################################################
#              Get categories of wildcards              #
#########################################################

SIMUDATASETNAMES = get_dataset_names(SIMUSNVSITESDIR, config["fineTune"]["indelsieve_simulator"]["subsample"])
FINETUNETYPESCELLPHY = get_fine_tune_types(
    config["fineTune"]["cellphy"],
    config["benchmark"]["fineTuneTypes"]
)
SIMUDATATYPES = get_data_types(
    config["benchmark"]["simulatedDataTypes"],
    config["benchmark"]["simulation"]["dataTypes"]["CNVTypes"]
)


########################################################
#                CELLPHY - From Monovar                #
########################################################

# Rule: run CELLPHY with filtered SNVs called by Monovar working with TRUE parameters (VCF input in EP mode)
rule runCellphyWithMonovarSNVsTrueParamsVCFEP:
    input:
        MONOVARTRUEPARAMSDIR + "{dataType}/{datasetName}/filtered.vcf"
    output:
        tree=protected(CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestTree"),
        model=protected(CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestModel")
    params:
        cellphy=config["tools"]["cellphy"],
        mode=lambda wildcards: get_dict_values(
            config,
            ["fineTune", "cellphy", wildcards.fineTuneTypeCellPhy, "mode"]
        ),
        seed=lambda wildcards: get_dict_values(
            config,
            ["fineTune", "cellphy", wildcards.fineTuneTypeCellPhy, "seed"]
        ),
        parsTrees=lambda wildcards: get_dict_values(
            config,
            ["fineTune", "cellphy", wildcards.fineTuneTypeCellPhy, "parsTrees"]
        ),
        prefix=CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}"
    threads:
        lambda wildcards, input: get_thread_num_for_cellphy(
            input[0],
            get_dict_values(
                config,
                ["fineTune", "cellphy", wildcards.fineTuneTypeCellPhy, "nrOfSitesPerThread"]
            )
        )
    benchmark:
        repeat(
            EFFICIENCYCELLPHYDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}.tsv",
            config["benchmark"]["efficiency"]["repetition"]
        )
    shell:
        "{params.cellphy} RAXML "
        "--{params.mode} "
        "--msa {input} "
        "--model GT10+FO+E+G "
        "--tree pars{{{params.parsTrees}}} "
        "--seed {params.seed} "
        "--thread {threads} "
        "--prefix {params.prefix} "
        "--prob-msa off "
        "--force model_lh_impr "
        "--force "
        "--redo "
        "&>{params.prefix}.out"

# Rule: replace the cell names on the best tree of CELLPHY with input data from Monovar for the sake of tree conparison
rule replaceCellphyCellNamesWithMonovarSNVsTrueParams:
    input:
        refCellNames=SIMUTUMORCELLNAMES,
        tree=CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestTree"
    output:
        protected(CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestTree_formatted")
    shell:
        "python3 scripts/format_cell_names_on_tree.py "
        "--tree {input.tree} "
        "--ref {input.refCellNames} "
        "--out {output}"


#######################################################
#              Gather trees from CELLPHY              #
#######################################################

# Rule: gather inferred trees from CELLPHY with filtered SNVs called by Monovar working with TRUE parameters
rule gatherCellphyTreesWithTrueMonovarSNVsEP:
    input:
        inferredTrees=expand(CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestTree_formatted", fineTuneTypeCellPhy=FINETUNETYPESCELLPHY, dataType=SIMUDATATYPES, datasetName=SIMUDATASETNAMES),
        trueTrees=expand(SIMUCONVERTEDTRUETREEDIR + config["benchmark"]["simulation"]["treesWithoutTrunkWithoutNormalDir"] + config["benchmark"]["simulation"]["treesPre"] + ".{datasetName}", datasetName=SIMUDATASETNAMES)
    output:
        protected(ANALYSISTREECOMPDIR + "from_cellphy_" + config["benchmark"]["analysis"]["trueMonovarSNVs"] + "_" + config["benchmark"]["analysis"]["cellphyEPMode"] + ".tsv")
    params:
        datasetNames=SIMUDATASETNAMES,
        trueTreesDir=SIMUCONVERTEDTRUETREEDIR + config["benchmark"]["simulation"]["treesWithoutTrunkWithoutNormalDir"],
        trueTreesPre=config["benchmark"]["simulation"]["treesPre"],
        tool="cellphy",
        snvType=config["benchmark"]["analysis"]["trueMonovarSNVs"],
        toolSetup=config["benchmark"]["analysis"]["cellphyEPMode"],
        fineTuneType=FINETUNETYPESCELLPHY,
        cnvRate=str(config["fineTune"]["simulateCNVs"]["prob"]),
        dataType=SIMUDATATYPES,
        cellNum=config["benchmark"]["analysis"]["cellNum"],
        covMean=config["benchmark"]["analysis"]["covMean"],
        covVariance=config["benchmark"]["analysis"]["covVariance"],
        effSeqErrRate=config["benchmark"]["analysis"]["effSeqErrRate"],
        adoRate=config["benchmark"]["analysis"]["adoRate"],
        deletionRate=config["benchmark"]["analysis"]["deletionRate"],
        insertionRate=config["benchmark"]["analysis"]["insertionRate"],
        doubletRate=config["benchmark"]["analysis"]["doubletRate"]
    script:
        "scripts/gather_inferred_trees_part.py"


########################################################
#              Get estimates from CELLPHY              #
########################################################

# Rule: get parameter estimates from CELLPHY with filtered SNVs called by Monovar working with TRUE parameters
rule getCellphyParamsWithTrueMonovarSNVsEP:
    input:
        expand(CELLPHYTRUEMONOVARSNVSDIR + "{fineTuneTypeCellPhy}/{dataType}/{datasetName}/vcf_ep/{datasetName}.raxml.bestModel", fineTuneTypeCellPhy=FINETUNETYPESCELLPHY, dataType=SIMUDATATYPES, datasetName=SIMUDATASETNAMES)
    output:
        protected(ANALYSISPARAMSDIR + "from_cellphy_" + config["benchmark"]["analysis"]["trueMonovarSNVs"] + ".tsv")
    params:
        cellNum=config["benchmark"]["analysis"]["cellNum"],
        covMean=config["benchmark"]["analysis"]["covMean"],
        covVariance=config["benchmark"]["analysis"]["covVariance"],
        datasetNames=SIMUDATASETNAMES,
        tool="cellphy",
        snvType=config["benchmark"]["analysis"]["trueMonovarSNVs"],
        toolSetup=config["benchmark"]["analysis"]["cellphyEPMode"],
        fineTuneType=FINETUNETYPESCELLPHY,
        cnvRate=str(config["fineTune"]["simulateCNVs"]["prob"]),
        dataType=SIMUDATATYPES,
        doubletRate=config["benchmark"]["analysis"]["doubletRate"]
    script:
        "scripts/get_param_estimates_cellphy.py"


#######################################################
#       Gather the number of sites from CELLPHY       #
#######################################################

# Rule: gather the number of variant sites as input to CELLPHY
rule gatherCellphySitesInfo:
    input:
        expand(MONOVARTRUEPARAMSDIR + "{dataType}/{datasetName}/filtered.vcf", dataType=SIMUDATATYPES, datasetName=SIMUDATASETNAMES)
    output:
        protected(ANALYSISSITESINFODIR + "from_cellphy_" + config["benchmark"]["analysis"]["trueMonovarSNVs"] + ".tsv")
    params:
        datasetNames=SIMUDATASETNAMES,
        tool="cellphy",
        snvType=config["benchmark"]["analysis"]["trueMonovarSNVs"],
        toolSetup=config["benchmark"]["analysis"]["cellphyEPMode"],
        fineTuneType=FINETUNETYPESCELLPHY,
        cnvRate=str(config["fineTune"]["simulateCNVs"]["prob"]),
        dataType=SIMUDATATYPES,
        cellNum=config["benchmark"]["analysis"]["cellNum"],
        covMean=config["benchmark"]["analysis"]["covMean"],
        covVariance=config["benchmark"]["analysis"]["covVariance"],
        effSeqErrRate=config["benchmark"]["analysis"]["effSeqErrRate"],
        adoRate=config["benchmark"]["analysis"]["adoRate"],
        deletionRate=config["benchmark"]["analysis"]["deletionRate"],
        insertionRate=config["benchmark"]["analysis"]["insertionRate"],
        doubletRate=config["benchmark"]["analysis"]["doubletRate"]
    script:
        "scripts/gather_sites_info_part.py"