Snakefile

## Snakefile


    ## module load snakemake/5.2.2 (also loads anaconda and python3)
    ## snakemake -np TARGET


from datetime import datetime


## globals
anno = '/scratch/chd5n/aneuploidy/raw-data/annotations/'
seq = '/scratch/chd5n/aneuploidy/raw-data/sequencing/'
count = '/scratch/chd5n/aneuploidy/raw-data/counts/'
plot = '/scratch/chd5n/aneuploidy/results/plots/'
hetdat = '/scratch/chd5n/aneuploidy/hetsites-data/'
dt = '_'.join(
    [
        str(datetime.now()).split(' ')[0],
        ''.join(str(datetime.now()).split(' ')[1].split('.')[0].split(':')[0:3])
    ]
)
set_num = 13  ## touch gdc_download.py to download next batch
v_num = 'v3'  ## set_size in gdc_download 100 in v2 from 20 in v1 to get 50 pairs instead of 10; v3 is custom run to fix errors
rna_set_name = 'mss_msil'  ## valid options: simple, mss, mss_msil
sva_pass_num = 2  ## first pass (1) or second pass (2)
sva_refine = 0  ## refinement no (0) or yes (1)
sva_nsv = 2  ## preferred number of svs to model


## rules
rule write_annotations:
    input:
        'gdc_write_anno.py'
    output:
        anno + 'pheno.tsv',
        anno + 'manifest.tsv'
    shell:
        'python {input}'


rule parse_info:
    input:
        parse = 'gdc_parse_info.py',
        prior_flag = anno + 'pheno.tsv'
    output:
        anno + 'coad-read.file_info',
        anno + 'coad-read.sub_errors'
    shell:
        'python {input.parse}'


rule download:
    input:
        dld = 'gdc_download.py',
        prior_flag = anno + 'coad-read.file_info'
    output:
        seq + 'latest_download.txt'
    shell:
        'nohup python {input.dld} --set_num {set_num} > ../logs/gdc_download_set_{set_num}{v_num}_{dt}.out 2> ../logs/gdc_download_set_{set_num}{v_num}_{dt}.err'


rule assemble_files:
    input:
        assemble = 'gdc_assemble_files.py',
        prior_flag = seq + 'latest_download.txt'
    output:
        anno + 'coad-read_current.file_set',
        anno + 'coad-read_current.file_err'
    shell:
        'python {input.assemble} --set_num {set_num} --datetime {dt}'


rule call_variants:
    input:
        gatk = 'run_gatk_haplo.py',
        file_set = anno + 'coad-read_current.file_set'
    output:
        seq + 'latest_gatk_haplo.txt'
    shell:
        'python {input.gatk} --in_file {input.file_set}'


rule count_hetsites:
    ## note: hetcnts_2R.py is deprecated
    input:
        hetcount = 'run_het_counter.py',
        file_set = anno + 'coad-read_current.file_set',
        prior_flag = seq + 'latest_gatk_haplo.txt'
    output:
        seq + 'latest_het_counts.txt'
    shell:
        'python {input.hetcount} --in_file {input.file_set}'


rule parse_hetcnts:
    ## aprox 3 min on frontend for 587 subjects
    ## this rule fixes a bug in hetcnts_2R.py
    ## note: output compressed and moved to hetsites-data manually
    input:
        py = 'parse_hetcnts.py',
        info = '/scratch/chd5n/aneuploidy/coad-read.file_info'
    output:
        end = '/scratch/chd5n/aneuploidy/hetcnts/TCGA-WS-AB45_hetcnts.tsv'
    shell:
        'python {input.py} --info {input.info} > parse_hetcnts_{dt}.log'


rule summarize_hetsites:
    ## this rule extracts total sites and unique chromosomes tested
    ## a quality control step to identify problems with pipeline
    ## 3-4 min on front end
    input:
        py = 'sum_sites.py',
        R = 'sum_sites.R'
    output:
        '/scratch/chd5n/aneuploidy/plots/summary_sites.pdf'
    shell:
        '''
        module load anaconda/5.2.0-py3.6 gcc/7.1.0 openmpi/3.1.4 R/4.0.0
        python {input.py}
        Rscript {input.R}
        '''


rule make_density_plots:
    input:
        plotter = 'run_plotter.py',
        file_set = anno + 'coad-read_current.file_set',
        prior_flag = seq + 'latest_het_counts.txt'
    output:
        plot + 'coad-read_current_hetcnts_plot.pdf'
    shell:
        'python {input.plotter} --in_file {input.file_set} --set_num {set_num} --datetime {dt}'


rule plot_hetcnts:
    ## this rule fixes some problems with former density plots
    ## approx 13 min on front end
    input:
        R = 'plot_hetcnts.R'
    output:
        end = '/scratch/chd5n/aneuploidy/plots/TCGA-WS-AB45_hetcnts.pdf'
    shell:
        '''
        module load gcc/7.1.0 openmpi/3.1.4 R/4.0.0
        Rscript {input.R} > plot_hetcnts_{dt}.log 2>&1
        '''


rule store_hetsite_data:
    input:
        store = 'store_hetsite_data.py',
        file_set = anno + 'coad-read_current.file_set',
        prior_flag = plot + 'coad-read_current_hetcnts_plot.pdf'
    output:
        hetdat + 'latest_data.txt'
    shell:
        'python {input.store} --in_file {input.file_set} --set_num {set_num} --datetime {dt}'


rule rank_sd:
    input:
        rank = 'rank_sd.R'
    shell:
        'module load gcc/7.1.0 R/3.6.1; Rscript {input.rank}'


rule download_de:
    input:
        dld = 'download_de.py'
    output:
        count + 'latest_download.txt'
    shell:
        'nohup python {input.dld} --set_name {rna_set_name} > ../logs/download_de_set_{rna_set_name}_{dt}.out 2> ../logs/download_de_set_{rna_set_name}_{dt}.err'


rule unzip_de:
    input:
        uz = 'unzip_de.py'
    shell:
        'python {input.uz} --set_name {rna_set_name} > ../logs/unzip_de_set_{rna_set_name}_{dt}.out 2> ../logs/unzip_de_set_{rna_set_name}_{dt}.err'


rule run_de:
    input:
        de = 'run_de.sh'
    shell:
        'sbatch {input.de} {rna_set_name} {sva_pass_num} {sva_refine} {sva_nsv}'


rule run_cor:
    input:
        cor = 'gene_fa_cor.R'
    shell:
        'module load gcc/7.1.0 R/3.6.1; Rscript {input.cor} --args {rna_set_name}'


rule plot_per_sample:
    input:
        pps = 'plot_per_sample.py'
    shell:
        'python {input.pps}'