Snakefile_lane

"""
Snakefile for one lane of single-cell fibroblast project (Raghd Rostom)

Author: Davis McCarthy
Affiliation: EMBL-EBI

Run: snakemake -s ../../../Snakefile_lane --jobs 400 --latency-wait 30 --cluster-config ../../../cluster.json --cluster 'bsub -J {cluster.name} -q {cluster.queue} -n {cluster.n} -R "select[singularity] rusage[mem={cluster.memory}]" -M {cluster.memory}  -o {cluster.output} -e {cluster.error}'

Davis McCarthy, 03 November 2017
"""

import glob
import os
from subprocess import run
import subprocess
import pandas as pd
import re
import h5py

shell.prefix("set -euo pipefail;") 

## reference files
human_tx_fasta = '../../../references/Homo_sapiens.GRCh37.rel75.cdna.all.ERCC92.fa.gz'
human_gx_fasta = '../../../references/GRCh37.p13.genome.ERCC92.fa'
fasta = human_tx_fasta
fasta_unzipped = human_gx_fasta
fasta_dict = fasta_unzipped.replace('fa', 'dict')
fasta_idx = fasta_unzipped + '.fai'
gencode_gtf = '../../../references/gencode.v19.annotation_ERCC.gtf'
dbSnpVcf = '../../../references/dbsnp_138.hg19.vcf.gz'
dbSnpVcfSmall = '../../../references/dbsnp_138.hg19.biallelicSNPs.HumanCoreExome12.Top1000ExpressedIpsGenes.Maf0.01.HWE0.0001.HipSci.vcf.gz'
reAlignmentIntervals = '../../../references/knownIndels.intervals'
knownIndelsMills = '../../../references/1000G_phase1.indels.hg19.sites.vcf.gz'
knownIndels100G = '../../../references/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz'
HIPSCI_VCF = '../../../references/hipsci.wec.gtarray.HumanCoreExome.imputed_phased.20170327.genotypes.allchr.fibro_samples_v2_filt_vars_sorted_oa.vcf.gz'

star_genome_files = ['chrLength.txt', 'chrNameLength.txt', 'chrName.txt', 'chrStart.txt', 'exonInfo.tab', 'Genome', 'genomeParameters.txt', 'SA', 'SAindex', 'sjdbInfo.txt', 'sjdbList.fromGTF.out.tab', 'sjdbList.out.tab', 'transcriptInfo.tab']
salmon_idx = '../../../references/Homo_sapiens.GRCh37.rel75.cdna.all.ERCC92.salmon_v0.8.2_idx'

## define commands
gatk_cmd = '/usr/bin/java -Xmx50g -Xms8g -Djava.io.tmpdir=tmp -jar /usr/GenomeAnalysisTK.jar'
read_salmon_to_scesets_cmd = '../../../src/R/read_salmon_to_scesets.R'
# rscript_cmd = 'singularity exec /hps/nobackup/hipsci/scratch/biocplus.img Rscript'

## parameter objects and samples
GENOME = 'genome.fa'
STAR_GENOME_DIR = '../../../references'
star_genome_files = ['chrLength.txt', 'chrNameLength.txt', 'chrName.txt', 'chrStart.txt', 'exonInfo.tab', 'Genome', 'genomeParameters.txt', 'SA', 'SAindex', 'sjdbInfo.txt', 'sjdbList.fromGTF.out.tab', 'sjdbList.out.tab', 'transcriptInfo.tab']

## targets
star_genome_output = expand('{genome_dir}/{genome_files}', genome_dir=STAR_GENOME_DIR, genome_files=star_genome_files)

## read in crams from SS2 run
crams_all = glob.glob('cram/*.cram')
## read in fastq from SS2 run
fastq_all = glob.glob('fastq/*_1.fastq')

## define sample names
# SAMPLES = [os.path.basename(w).replace('.cram', '') for w in crams_all]
SAMPLES = [os.path.basename(w).replace('_1.fastq', '') for w in fastq_all]
fastqc_html_reports = expand('fastqc/{sample}.2pass.Aligned.sortedByCoord.out_fastqc.html', sample=SAMPLES)
cell_vcf_files = expand('vcf/{sample}/{sample}.filtered.vcf.gz', sample=SAMPLES)
merged_bams = expand('bam/merged/{sample}.bam', sample=SAMPLES)
salmon_results_GRCh37 = expand('quant_salmon/{sample}/quant.sf', sample=SAMPLES)
sces = expand('sces/sce.salmon.{genome}.preqc_tx.rds', genome=['GRCh37'])
sces.append(expand('sces/sce.salmon.{genome}.preqc_gene.rds', genome=['GRCh37']))
scater_first_html_reports = expand('salmon.{genome}.first_qc.html', genome=['GRCh37'])

    
rule all:
    input:
        sces,
        scater_first_html_reports,
        cell_vcf_files
    
    
rule build_star_genome_indexes:
    input:
        fasta=fasta_unzipped,
        annotation=gencode_gtf
    output:
        star_genome_output
    threads: 8
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'STAR --runMode genomeGenerate --genomeDir {STAR_GENOME_DIR} '
        '--genomeFastaFiles {input.fasta} --runThreadN {threads} '
        '--sjdbGTFfile {input.annotation} --sjdbOverhang 100'


rule create_fasta_index:
    input:
        fasta_unzipped
    output:
        fasta_idx
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'samtools faidx {input}'


rule create_fasta_dict:
    input:
        fasta_unzipped
    output:
        fasta_dict
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'picard CreateSequenceDictionary R={input} O={output}'


rule build_salmon_index:
    input:
        fasta
    output:
        salmon_idx
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'salmon index -t {input} -i {output} --type quasi -k 31 '
        '--perfectHash'


rule cram2fastq:
    input:
        'cram/{sample}.cram'
    output:
        fq1=temp('fastq/{sample}_1.fastq'),
        fq2=temp('fastq/{sample}_2.fastq')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        """
        samtools view -u {input} | \
        samtools collate -uOn 128 - tmp/tmp-prefix-{wildcards.sample} | \
        samtools fastq -F 0xB00 -1 {output.fq1} -2 {output.fq2} -
        """


rule trim_fastq_SS2:
    input:
        fq1="fastq/{sample}_1.fastq",
        fq2="fastq/{sample}_2.fastq"
    output:
        fq1='fastq/{sample}_1_val_1.fq.gz',
        fq2='fastq/{sample}_2_val_2.fq.gz'
    priority: 6
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'trim_galore --gzip --output_dir fastq --length 40 --paired '
        '{input.fq1} {input.fq2}'


rule salmon_quant:
    input:
        sidx=salmon_idx,
        fq1='fastq/{sample}_1_val_1.fq.gz',
        fq2='fastq/{sample}_2_val_2.fq.gz'
    output:
        'quant_salmon/{sample}/quant.sf'
    threads: 4
    params:
        folder='quant_salmon/{sample}/'
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'salmon quant -i {input.sidx} -l IU '
        '-1 {input.fq1} -2 {input.fq2} '
        '--seqBias --gcBias --threads {threads} --useVBOpt ' 
        '-o {params.folder}'


rule align_with_star_2pass_SS2:
    input:
        star_genome_output,
        genome_dir=STAR_GENOME_DIR,
        fq1='fastq/{sample}_1_val_1.fq.gz',
        fq2='fastq/{sample}_2_val_2.fq.gz'
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.out.bam')
    params: 
        prefix='star/{sample}/{sample}.2pass.'
    threads: 8
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'STAR --genomeDir {input.genome_dir} '
        '--readFilesIn {input.fq1} {input.fq2} '
        '--outFileNamePrefix {params.prefix} '
        '--outSAMtype BAM Unsorted ' 
        '--alignSJoverhangMin 8 ' 
        '--alignSJDBoverhangMin 1 --genomeLoad NoSharedMemory '
        '--alignIntronMin 20 --alignIntronMax 1000000 '
        '--alignMatesGapMax 1000000 --sjdbScore 2 '
        '--outFilterType BySJout '
        '--outFilterMultimapNmax 20 --outFilterMismatchNmax 999 '
        '--outFilterMismatchNoverLmax 0.04 '
        '--outFilterScoreMinOverLread 0.33 --outFilterMatchNminOverLread 0.33 '
        '--outSAMstrandField intronMotif '
        '--outFilterIntronMotifs RemoveNoncanonical '
        '--outSAMattributes NH HI NM MD AS XS --outSAMunmapped Within '
        '--runThreadN {threads} --twopassMode Basic '
        '--readFilesCommand zcat ' 


rule sort_alignment:
    input:
        'star/{sample}/{sample}.2pass.Aligned.out.bam'
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.out.bam')
    threads: 4
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'samtools sort -m 12G -@ 4 -O bam -o {output} {input}'


rule index_star_bams:
    input:
        'star/{sample}/{sample}.2pass.Aligned.sortedByCoord.out.bam'
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.out.bam.bai')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'samtools index {input} '


rule picard_read_groups:
    input:
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.out.bam'
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.rgadded.bam')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'picard AddOrReplaceReadGroups I={input.bam} O={output} SO=coordinate '
        'RGID={wildcards.sample} RGLB={wildcards.sample} '
        'RGPL=ILLUMINA RGPU=MACHINE1 RGSM={wildcards.sample}'


rule picard_mark_dups:
    input:
        'star/{sample}/{sample}.2pass.Aligned.sortedByCoord.rgadded.bam'
    output:
        bam=temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.dedup.bam'),
        metrics='star/{sample}/{sample}.output.metrics'
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        'picard MarkDuplicates I={input} O={output.bam} CREATE_INDEX=true '
        'VALIDATION_STRINGENCY=SILENT M={output.metrics} '


rule split_n_trim_gatk:
    input:
        fasta=fasta_unzipped,
        fai=fasta_idx,
        fa_dict=fasta_dict,
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.dedup.bam'
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.bam')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T SplitNCigarReads -R {input.fasta} -I {input.bam} -o {output} '
        '-rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 '
        '-U ALLOW_N_CIGAR_READS '


rule indel_realignment_gatk:
    input:
        fasta=fasta_unzipped,
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.bam',
        targetIntervals=reAlignmentIntervals,
        known1=knownIndelsMills,
        known2=knownIndels100G
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bam')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T IndelRealigner -R {input.fasta} -I {input.bam} '
        '-targetIntervals {input.targetIntervals} -known {input.known1} -known {input.known2} '
        '-U ALLOW_N_CIGAR_READS --consensusDeterminationModel KNOWNS_ONLY --LODThresholdForCleaning 0.4  '
        '-o {output}' 


rule base_recalibrator_gatk:
    input:
        fasta=fasta_unzipped,
        dbSnp=dbSnpVcf,
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bam',
        known1=knownIndelsMills,
        known2=knownIndels100G
    output:
        temp('star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bqsr')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T BaseRecalibrator -R {input.fasta} -I {input.bam} '
        '-knownSites {input.known1} -knownSites {input.known2} -knownSites {input.dbSnp} '
        '-nct 2 '
        '-o {output}' 


rule recalibrated_writer_gatk:
    input:
        fasta=fasta_unzipped,
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bam',
        bqsr='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bqsr'
    output:
        'star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bqsr.bam'
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T PrintReads -R {input.fasta} -I {input.bam} '
        '-BQSR {input.bqsr} -nct 2 '
        '-o {output}' 


rule call_variants_gatk:
    input:
        fasta=fasta_unzipped,
        fai=fasta_idx,
        fa_dict=fasta_dict,
        dbSnp=dbSnpVcf,
        dbSnpSmall= dbSnpVcfSmall,
        bam='star/{sample}/{sample}.2pass.Aligned.sortedByCoord.split.realigned.bqsr.bam'
    output:
        temp('vcf/{sample}/{sample}.unfiltered.vcf')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T HaplotypeCaller -R {input.fasta} -I {input.bam} '
        '-dontUseSoftClippedBases '
        '-D {input.dbSnp} -gt_mode GENOTYPE_GIVEN_ALLELES '
        '-alleles {input.dbSnpSmall} -L {input.dbSnpSmall} -o {output}' 


rule filter_variants_gatk:
    input:
        vcf='vcf/{sample}/{sample}.unfiltered.vcf',
        fasta=fasta_unzipped,
        fai=fasta_idx,
        fa_dict=fasta_dict
    output:
        temp('vcf/{sample}/{sample}.filtered.vcf')
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://broadinstitute/gatk3:3.8-1"
    shell:
        '{gatk_cmd} -T VariantFiltration -R {input.fasta} -V {input.vcf} '
        '-window 35 -cluster 3 -filterName FS -filter "FS > 30.0" '
        '-filterName QD -filter "QD < 2.0" -o {output}'


rule bgzip_vcf:
    input:
        'vcf/{sample}/{name}.filtered.vcf'
    output:
        f1='vcf/{sample}/{name}.filtered.vcf.gz',
        f2='vcf/{sample}/{name}.filtered.vcf.gz.csi'
    conda:
        "envs/myenv.yaml"
    singularity:
        "docker://davismcc/fibroblast-clonality"
    shell:
        """
        vcf-sort {input} | bgzip -c > {output.f1}
        bcftools index {output.f1}
        """


rule salmon_to_sceset_GRCh37:
    input:
        files=salmon_results_GRCh37
    output:
        'sces/sce.salmon.GRCh37.preqc_tx.rds',
        'sces/sce.salmon.GRCh37.preqc_gene.rds',
        #'sces/sce.salmon.GRCh37.preqc.feather'
    params:
        input_dir='quant_salmon',
        output_prefix='sces/sce.salmon.GRCh37.preqc'
    singularity:
        "docker://davismcc/r-singlecell-img:0.1"
    shell:
        'Rscript {read_salmon_to_scesets_cmd} '
        '--input_dir {params.input_dir} '
        '--output_prefix {params.output_prefix} '
        '--biomart feb2014.archive.ensembl.org'


rule rough_qc:
    input:
        'sces/sce.salmon.{genome}.preqc_gene.rds'
    output:
        'salmon.{genome}.first_qc.html'
    singularity:
        "docker://davismcc/r-singlecell-img:0.1"
    shell:
        'Rscript ../../../src/R/compile_report.R '
        '-i {input} -o {output} '
        '--template ../../../src/Rmd/rough_qc_template.Rmd '