Set up my Snake and Config files

Rules to download and index genome

This is from John’s cardioQTL snake file

rule download_genome:
    output: dir_genome + "Homo_sapiens." + ensembl_genome + \
            ".dna_sm.chromosome.{chr}.fa.gz"
    params: chr = "{chr}", build = ensembl_genome
    shell: "wget -O {output} {ensembl_ftp}Homo_sapiens.{params.build}.dna_sm.chromosome.{params.chr}.fa.gz"

rule unzip_chromosome_fasta:
    input: dir_genome + "Homo_sapiens." + ensembl_genome + \
           ".dna_sm.chromosome.{chr}.fa.gz"
    output: temp(dir_genome + "Homo_sapiens." + ensembl_genome + \
                 ".dna_sm.chromosome.{chr}.fa")
    shell: "zcat {input} > {output}"

rule subread_index:
    input: expand(dir_genome + "Homo_sapiens." + ensembl_genome + \
                  ".dna_sm.chromosome.{chr}.fa", \
                  chr = chr_genes)
    output: dir_genome + ensembl_genome + ".reads"
    params: prefix = dir_genome + ensembl_genome
    shell: "subread-buildindex -o {params.prefix} {input}"

Fastqc rule

First add a path to the fastqc dir.

fastqc_dir = output + "fastqc/

rule fastqc:
    input:
        fastq + "{samples}.fastq"
    output:
        fastqc_dir + "{samples}_fastqc.html",
        fastqc_dir + "{samples}_fastqc.zip"
    params:
        outdir = fastqc_dir
    shell:
        "fastqc -o {params.outdir} {input}"

Unzip fastqc

Maybe I dont need to do this right now: This is code to maybe add later:

rule unzip:
    input:
        fastqc_dir + "{sample}_fastqc.zip"
    output:
        fastqc_dir + "{sample}_fastqc.zip/fastqc_data.txt"
    params:
        outdir = fastqc_dir
    shell:
        "unzip -d {params.outdir} {input}"

Extract umi data

This is the next step to figure out. I need to trim the primer and extract the barcode from the linker.

Barcode DNA linker /5rApp/(N)6CTGTAGGCACCATCAAT/3ddC

Use UMI tools to extract first 6 for each read

rule umi_extract_umi:
    input:
        fastq_dir + "{samples}.fastq",
    output:
        fastq_extr + "{samples}_extracted.fastq",
    params:
        bc = "NNNNNN""
    shell:
        "umi_tools extract --bc-pattern={params.bc} --stdin {input} --stdout {output}"

Need to tell it which samples to use:

samples = set(glob_wildcards(fastq_dir + "{samples}.fastq").samples)

Align with Subread-subjunc

rule subjunc:
    input: 
        read = fastq_extr + "{samples}_extracted.fastq",
        index = dir_genome + ensembl_genome + ".reads"
    output: 
        dir_bam + "{samples}.bam"
    params: 
        prefix = dir_genome + ensembl_genome
    threads: 8
    shell: "subjunc -i {params.prefix} -r {input.read} -T {threads} > {output}"

rule sort_bam:
    input: 
        dir_bam + "{samples}.bam"
    output:
        dir_sort + "{samples}-sort.bam"
    shell: "samtools sort -o {output} {input}"   
    
    
rule index_bam:
    input: 
        dir_sort + "{samples}-sort.bam"
    output:
        dir_sort + "{samples}-sort.bam.bai"
    shell: "samtools index {input}"
    
#add folder
dir_bam= dir_data + "bam/"
dir_sort= dir_data + "sort/"

Rule dedup

rule dedup_bam:
    input: 
        bam = dir_sort + "{samples}-sort.bam",
        index = dir_sort + "{samples}-sort.bam.bai"
    output:
        dir_dedup + "{samples}-sort.dedup.bam"
    shell: "umi_tools dedup -I {input.bam} -S {output}"
    
    
#dir to add
dir_dedup=dir_data + "dedup/"

Rule genomecov

I will use this bedtools command to visualize the densities of netseq reads.

http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html

bedtools genomecov -bg -5 -i fastq -g indexed genome

Flags to use:

• -bg : genome wide coverage in a bedgraph format

• -5 calculate coverage for the 5` end of the read

rule genome_cov:
    input:
       bamdedup = dir_dedup + "{samples}-sort.dedup.bam",
       genome = dir_genome + ensembl_genome + ".reads"
    output:
       dir_cov + "{samples}-sort.dedup.cov.bed"
    shell: "bedtools genomecov -bg -5 -ibam {input.bamdedup} -g {input.genome}" > {output}



#dir to add
dir_cov= dir_data + "cov/"

Works when i run it for just the file:

bedtools genomecov -bg -5 -ibam data/dedup/SRR1575922_CUT-sort.dedup.bam -g genome/GRCh37.75.reads

Problem:

This snake rule doesnt know where to write the output. I am not sure what file type this should be.

Rule all:

rule all:
    input:
       expand(fastqc_dir + "{samples}_fastqc.html", samples = samples), 
       dir_genome + ensembl_genome + ".reads",
       expand(fastq_extr + "{samples}_extracted.fastq", samples = samples),
       expand(dir_bam + "{samples}.bam", samples = samples),
       expand(dir_sort + "{samples}-sort.bam", samples= samples),
       expand(dir_sort + "{samples}-sort.bam.bai", samples=samples),
       expand(dir_dedup + "{samples}-sort.dedup.bam", samples=samples),
       expand(dir_cov + "{samples}-sort.dedup.cov.bed", samples=samples)
       

Cluster

Will need to create a cluster.json file where I specify the memory for each of my rules

Create for first 2 rules:

{
    "__default__" :
    {
        "mem" : 2000,
        "n" : 1,
    "tasks" : 1,
    "name" : "{rule}-{wildcards}",
    "logfile": "log/snake-{rule}-{wildcards}-%j.out"
    },

"subread_index" :
    {
            "mem" : 8000
    },
"fastqc" :
    {
              "mem" : 8000
    }
"subjunc" :
    {
             "mem": 16000,
             "tasks" : 8
    },
"sort_bam" :
    {
             "mem": 8000
    }

}

Full files are in my pipeline folder on midway

Run

To check the snakefile:
snakefile -np
To run the snakefile, in the pipleline directory:
scripts/submit-snakemake.sh
to check running jobs:
squeue -user=brimittleman

Things to Change

• add a get exons rule/download exon file