tapestry/Snakefile at main · quinlan-lab/tapestry · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
# [WIP]

# TODO: complete this snakemake pipeline for companion paper
# Using the following scripts as "context", ask Gemini to complete the pipeline:
# - run-hiphase.sh
# - build-iht-based-haplotype-map-and-phase-variants.sh
# - aligned_bam_to_cpg_scores.sh
# - phase_meth_to_founder_haps.sh

# Resources:
# https://snakemake.readthedocs.io/en/stable/index.html
# https://github.com/quinlan-lab/Snakemake_Tutorial

# Installation Problem:
# https://uofu.service-now.com/it?sys_id=6273777683206290ba4da250ceaad3c3&view=sp&id=ticket_history&table=incident
# https://quinlangroup.slack.com/archives/D9LFRMXV3/p1743466049592919

# Installation Resolution:
# https://github.com/snakemake/snakemake/releases/tag/v9.2.1

# Command to run:
# snakemake --cores 16 --snakefile dna-methylation/Snakefile

input_dir = "/scratch/ucgd/lustre-core/UCGD_Datahub/Mosaic/1654/UCGD/GRCh38/LongRead/Data/PolishedCrams"
reference = "/scratch/ucgd/lustre-core/common/data/Reference/homo_sapiens/GRCh38/primary_assembly_decoy_phix.fa"
output_dir = "/scratch/ucgd/lustre-labs/quinlan/data-shared/dna-methylation/medgenome-model-mode-test-snakemake"
bin_dir = "/uufs/chpc.utah.edu/common/HIPAA/u6018199/pb-CpG-tools-v3.0.0-x86_64-unknown-linux-gnu/bin/"

sys.path.append(f'/scratch/ucgd/lustre-labs/quinlan/u6018199/medgenome-hifi-k1375-hackathon/dna-methylation/util')

from get_id_to_paths import get_id_to_paths_medgenome_pilot, get_uids

experiment = get_id_to_paths_medgenome_pilot()
samples = get_uids(experiment)
for sample in samples:
    print(sample)

rule all:
    input:
        expand(
            "{output_dir}/{sample}.GRCh38.haplotagged.combined.bed.gz.tbi",
            output_dir=output_dir,
            sample=samples
        )
    shell:
        """
        echo {input}
        """

rule compute_methylation_level:
    log:
        # TODO: is expand necessary, to indicate that output_dir is not a wildcard
        # or could I make output_dir a snakemake param?
        stderr = expand(
            "{output_dir}/snakemake-logs/{sample}-compute_methylation_level.log",
            output_dir=output_dir,
            sample=samples
        )
    input:
        bam = expand(
            "{input_dir}/{sample}.cram",
            input_dir=input_dir,
            sample=samples
        ),
        ref = reference
    output:
        # TODO:
        # does this output have to match the input of "rule all"
        tbi = expand(
            "{output_dir}/{sample}.GRCh38.haplotagged.combined.bed.gz.tbi",
            output_dir=output_dir,
            sample=samples
        )
    shell:
        """
        echo "hello"
        touch {output.tbi}
        """
        # TODO: make bin_dir a snakemake param?

        # echo $bin_dir
        # aligned_bam_to_cpg_scores --help
        # export PATH=${bin_dir}:$PATH
        # aligned_bam_to_cpg_scores \
        #     --bam {input.bam} \
        #     --ref {input.ref} \
        #     --output-prefix {output.prefix} \
        #     --threads 8 \
        #     --min-coverage 10 \
        #     --min-mapq 1 \
        #     --pileup-mode model \
        #     2> {log.stderr}