Snakefile_fastq_multilane.MSI1

include:
    'configs/config_Apr_04_2013_MSI1_RiboSeq.human.py'

workdir: OUT_DIR

from itertools import chain
from os.path import join
import glob
import re

SAMPLES = glob.glob('{}**/*.fastq.bz2'.format(RAWDATA_DIR), recursive=True)
SAMPLE_LANE = []
PROJECT_ID = []
SAMPLE_ID = []
for sample in SAMPLES:
    sample = sample.replace('{}'.format(RAWDATA_DIR),'')
    project_id, sample_id, filename = sample.split('/')
    #sample_name = re.split(r'_L\d\d\d_', filename)[0]
    lane_name = re.search(r'R1_\d\d\d', filename).group()
    sample_name = re.split(r'_R1_\d\d\d', filename)[0]
    SAMPLE_LANE.append((project_id, sample_id, sample_name, lane_name))

SAMPLE_LANE = list(set(SAMPLE_LANE))
SAMPLE_LANE = list(set(sorted(SAMPLE_LANE, key=lambda tup: tup[0])))
PROJECT_ID, SAMPLE_ID, SAMPLE_NAMES, LANE_NAMES = list(set(zip(*SAMPLE_LANE)))
PROJECT_ID = [x[0] for x in SAMPLE_LANE]
SAMPLE_ID = [x[1] for x in SAMPLE_LANE]
SAMPLE_NAMES = [x[2] for x in SAMPLE_LANE]
LANE_NAMES = [x[3] for x in SAMPLE_LANE]
#SAMPLE_NAMES = sorted(set(SAMPLE_NAMES))

print(len(SAMPLE_LANE))
rule all:
    input:
        STAR_INDEX,
        expand('preprocessed/gzip/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.gz', zip,
                project_id = PROJECT_ID,
                sample_id =  SAMPLE_ID,
                sample_name = SAMPLE_NAMES,
                lane = LANE_NAMES),
        #expand('qc/{project_id}/{sample_id}/{sample_name}_{lane}_fastqc.html', zip, 
        #        project_id = PROJECT_ID,
        #        sample_id =  SAMPLE_ID,
        #        sample_name = SAMPLE_NAMES,
        #        lane = LANE_NAMES),
        #expand('inferred_experiment/{sample}.txt', sample=SAMPLE_NAMES),
        #expand('preprocessed/{project_id}/{sample_id}/{sample_name}_{lane}_trimmed.fq.gz', zip,
        #        project_id = PROJECT_ID,
        #        sample_id =  SAMPLE_ID,
        #        sample_name = SAMPLE_NAMES,
        #        lane = LANE_NAMES),
        #expand('mapped/bams/{project_id}/{sample_id}/{sample_name}_{lane}.bam', zip,
        #        project_id = PROJECT_ID,
        #        sample_id =  SAMPLE_ID,
        #        sample_name = SAMPLE_NAMES,
        #        lane = LANE_NAMES),

        #expand('mapped/merged_bams/{sample_id}.bam', sample_id = SAMPLE_ID),
        #expand('mapped/HTSeq-counts/byExon/{sample}.exon.counts.tsv', sample=SAMPLE_NAMES),
        #expand('mapped/HTSeq-counts/byCDS/{sample}.CDS.counts.tsv', sample=SAMPLE_NAMES),
        #'mapped/DE_analysis/HTSeq/'+GENOME_BUILD+'.HTSeq.DESeq2.all.tsv',
        #'mapped/DE_analysis/featureCounts/'+GENOME_BUILD+'.featureCounts.DESeq2.sig.tsv',
        #'multiqc_report/multiqc_report.html',

rule convert_to_bzip2:
    ## NOt sure bzip2 is standard input for all tools
    input:
        R1 = RAWDATA_DIR + '/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.bz2'
    output:
        'preprocessed/gzip/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.gz'
    shell:
        '''bunzip2 -c < {input} | gzip -c > {output}'''

rule perform_qc:
    input:
        R1='preprocessed/gzip/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.gz',
    params:
        out_dir = 'qc'
    output:
       'qc/{project_id}/{sample_id}/{sample_name}_{lane}_fastqc.html',
       'qc/{project_id}/{sample_id}/{sample_name}_{lane}_fastqc.zip',
    shell:
        r'''
            fastqc -o {params.out_dir} -f fastq {input.R1}
        '''

rule perfom_trimming:
    input:
        R1='preprocessed/gzip/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.gz',
    params:
        out_dir='preprocessed/{project_id}/{sample_id}',
        phred_cutoff=5
    output:
        'preprocessed/{project_id}/{sample_id}/{sample_name}_{lane}_trimmed.fq.gz'
    shell:
        r'''
            trim_galore -o {params.out_dir} -q {params.phred_cutoff} {input.R1}
        '''

rule map_star:
    input:
        R1='preprocessed/gzip/{project_id}/{sample_id}/{sample_name}_{lane}.fastq.gz',
        index=STAR_INDEX
    output: 'mapped/bams/{project_id}/{sample_id}/{sample_name}_{lane}.bam'
    params:
        prefix = 'mapped/bams/{project_id}/{sample_id}',
        unmapped = 'unmapped/fastq/',
        starlogs = 'mapped/starlogs'
    threads: 16
    shell:
        r'''
        STAR --runThreadN {threads}\
             --genomeDir {input.index}\
             --outFilterMismatchNmax 2\
             --outFileNamePrefix {params.prefix} --readFilesIn {input.R1}\
             --outSAMtype BAM SortedByCoordinate\
             --readFilesCommand zcat\
             --outReadsUnmapped {params.unmapped} && mv {params.prefix}Aligned.sortedByCoord.out.bam {output} && mkdir -p {params.starlogs} && mv {params.prefix}Log.final.out {params.prefix}Log.out {params.prefix}Log.progress.out {params.starlogs}
        '''

## See: https://software.broadinstitute.org/gatk/guide/article?id=3060
## Merging should happen post alignment
rule merge_bams:
    input: expand('mapped/bams/{project_id}/{sample_id}/{sample_name}_{lane}.bam', project_id=PROJECT_ID,sample_id=SAMPLE_ID,sample_name=SAMPLE_NAMES,lane=LANE_NAMES)
    #input: 'mapped/bams/{project_id}/{sample_id}/{sample_name}_{lane}.bam'
    output: expand('mapped/merged_bams/{sample_id}.bam', sample_id=SAMPLE_ID),
    run:
        projects = {}
        for i in input:
            filepath = i.replace('mapped/bams/', '')
            project_id, sample_id, filename = filepath.split('/')
            if project_id not in projects:
                projects[project_id] = {}
            if sample_id not in projects[project_id]:
                projects[project_id][sample_id] = []
            projects[project_id][sample_id].append(filename)
        for key, value in projects.iteritems():
            cmd = ' -in '.join(value)
            outfile = 'mapped/merged_bams/'+ value+ '.bam'
            shell(r'''bamtools merge -in {cmd} -out {outfile}''')

rule sort_by_name:
    input: 'mapped/merged_bams/{sample_id}.bam'
    output: 'mapped/merged_bams/{sample_id}.sortedByName.bam'
    shell:
        r'''
            samtools sort -on {input} -T /tmp/ -o {output}
        '''

rule count_exon:
    input: 'mapped/merged_bams/{sample_id}.sortedByName.bam'
    params:
        annotation=GTF,
        phred_cutoff=5
    output: 'mapped/HTSeq-counts/byExon/{sample_id}.exon.counts.tsv'
    shell:
        r'''
        source activate clipseq2 && htseq-count --order=name --format=bam --mode={HTSEQ_MODE} --stranded={HTSEQ_STRANDED} --minaqual={params.phred_cutoff} --type=exon --idattr=gene_id {input} {params.annotation} > {output}
        '''

rule count_cds:
    input: 'mapped/merged_bams/{sample_id}.sortedByName.bam'
    params:
        annotation=GTF,
        phred_cutoff=5
    output: 'mapped/HTSeq-counts/byCDS/{sample_id}.CDS.counts.tsv'
    shell:
        r'''
        source activate clipseq2 && htseq-count --order=name --format=bam --mode={HTSEQ_MODE} --stranded={HTSEQ_STRANDED} --minaqual={params.phred_cutoff} --type=CDS --idattr=gene_id {input} {params.annotation} > {output}
        '''

rule format_counts:
    input: 'mapped/HTSeq-counts/byExon/{samplei_id}.exon.counts.tsv'
    output: 'mapped/HTSeq-counts/byExon/{sample_id}.counts.noversion.tsv'
    shell:
        r'''
        cat {input} | sed -E 's/\.[0-9]+//' > {output}

        '''

rule featurecounts:
    input: expand('mapped/merged_bams/{sample_id}.sortedByName.bam', sample_id=set(SAMPLE_NAMES))
    params:
        annotation=GTF
    output: 'mapped/featureCounts/fcounts.tsv'
    threads: 16
    shell:
        r'''featureCounts {FEATURECOUNTS_S} -a {params.annotation} -o {output} -t exon -g gene_id -Q 4 -T {threads} {input}'''

rule format_fcounts:
    input: 'mapped/featureCounts/fcounts.tsv'
    output: 'mapped/featureCounts/fcounts.noversion.tsv'
    shell:
        r'''
        cat {input} | sed -E 's/\.[0-9]+//' > {output}

        '''
    output: 'mapped/featureCounts/fcounts.noversion.tsv'

rule run_picardmetrics:
    input: 'mapped/merged_bams/{sample_id}.bam'
    output: 'mapped/bam_metrics/{sample_id}.metrics'
    shell:
        r'''
        picard CollectInsertSizeMetrics I={input} H={output}.insertsize.pdf O={output}

        '''

rule run_deseq:
    input: expand('mapped/HTSeq-counts/byExon/{sample_id}.counts.noversion.tsv', sample_id=SAMPLE_ID)
    output:
        'mapped/DE_analysis/HTSeq/'+GENOME_BUILD+'.HTSeq.DESeq2.all.tsv',
        'mapped/DE_analysis/HTSeq/'+GENOME_BUILD+'.HTSeq.DESeq2.sig.tsv'

    params:
        basedir = 'mapped/HTSeq-counts/byExon',
        inprefix = 'counts.noversion',
        gene_annotations  = GENE_NAMES,
        outprefix = 'mapped/DE_analysis/HTSeq/'+GENOME_BUILD+'.HTSeq'
    shell:
        r'''
        Rscript {SRC_DIR}/do_DE_analysis.R --basedir={params.basedir} \
            --gene_annotations={params.gene_annotations} \
            --design_file={DESIGN_FILE} \
            --outprefix={params.outprefix} \
            --inprefix={params.inprefix}

        '''

rule run_deseq_featureCounts:
    input: 'mapped/featureCounts/fcounts.noversion.tsv'
    output:
        'mapped/DE_analysis/featureCounts/'+GENOME_BUILD+'.featureCounts.DESeq2.all.tsv',
        'mapped/DE_analysis/featureCounts/'+GENOME_BUILD+'.featureCounts.DESeq2.sig.tsv'
    params:
        gene_annotations  = GENE_NAMES,
        outprefix = 'mapped/DE_analysis/featureCounts/'+GENOME_BUILD+'.featureCounts'
    shell:
        r'''
        Rscript {SRC_DIR}/do_DE_analysis_featureCounts.R --counts={input} \
            --gene_annotations={params.gene_annotations} \
            --design_file={DESIGN_FILE} \
            --outprefix={params.outprefix}

        '''


rule create_insertsize_tsv:
    input: 'mapped/bam_metrics/{sample}.metrics'
    output: 'mapped/bam_metrics/{sample}.insertsizes.tsv'
    shell:
        r'''
        python {SRC_DIR}/collect_picard_metrics.py {input} {output}

        '''
rule counts_to_tpm:
    input:
        count = expand('mapped/HTSeq-counts/byExon/{sample}.counts.noversion.tsv', sample=SAMPLE_NAMES),
        insert_size = expand('mapped/bam_metrics/{sample}.insertsizes.tsv', sample=SAMPLE_NAMES),
    output: 
        'mapped/tpm/masterTPM.tsv'
    params:
        gene_lengths=GENE_LENGTHS,
        name=expand('{sample}', sample=SAMPLE_NAMES),
        outprefix='mapped/tpm',
        gene_map=GENE_NAMES
    run:
        counts_input = (',').join(input.count)
        sizes_input = (',').join(input.insert_size)
        names = (',').join(params.name)
        shell('Rscript {SRC_DIR}/counts_to_tpm.R --counts={counts_input} --insert_sizes={sizes_input} --gene_lengths={params.gene_lengths} --inprefix={names} --gene_map={params.gene_map} --outprefix={params.outprefix}')

rule plot_tpm:
    input: 'mapped/tpm/masterTPM.tsv'
    output: 'mapped/plots/tpm_scatter/{sample1}_VS_{sample2}.png'
    params:
        prefix = 'mapped/plots/tpm_scatter/'
    shell:
        r'''
        python {SRC_DIR}/plot_tpm_scatter.py --master {input} --outprefix  {params.prefix}
        '''

rule perform_qualimap_qc:
    input:  'mapped/merged_bams/{sample_id}.bam',
    output: 'mapped/post_mapping_qualimap/{sample_id}/qualimapReport.html',
    params:
        outdir='mapped/post_mapping_qualimap/{samplei_id}',
        gtf=GTF
    shell:
        r'''
        qualimap rnaseq -bam {input} -gtf {params.gtf} --outdir {params.outdir} --java-mem-size=8G

        '''

rule get_duplication_estimate:
    input:  'mapped/merged_bams/{sample_id}.bam'
    output: 'mapped/post_mapping_deduplication/{sample}/output.DupRate_plot.r'
    params:
        outprefix='mapped/post_mapping_deduplication/{sample}/output'
    shell:
        r'''
        source activate clipseq2 && read_duplication.py -i {input} -o {params.outprefix}
        '''

rule infer_experiment:
    input: 'mapped/merged_bams/{sample_id}.bam'
    output: 'inferred_experiment/{sample}.txt'
    shell:
        r'''
        source activate clipseq2 && infer_experiment.py -r {GENE_BED} -i {input} 2>&1 > {output} 
        '''

rule run_multiqc:
    input:
        expand('qc/{sample_name}_{lane}_R1_001_fastqc.html', zip, sample_name=SAMPLE_NAMES, lane=LANE_NAMES),
        expand('mapped/post_mapping_deduplication/{sample}/output.DupRate_plot.r', sample=SAMPLE_NAMES),
        expand('mapped/post_mapping_qualimap/{sample}/qualimapReport.html', sample=SAMPLE_NAMES),
        expand('mapped/merged_bams/{sample_id}.bam', sample_id=SAMPLE_ID),
        expand('inferred_experiment/{sample}.txt', sample=SAMPLE_NAMES),
        'mapped/DE_analysis/HTSeq/'+GENOME_BUILD+'.HTSeq.DESeq2.sig.tsv',
        'mapped/DE_analysis/featureCounts/'+GENOME_BUILD+'.featureCounts.DESeq2.sig.tsv',
    output:
        'multiqc_report/multiqc_report.html'
    shell:
        'export LC_ALL=en_US.UTF-8 && multiqc -f --outdir multiqc_report .'