Update Snakefile for pilot analysis

Genome location

My old snakefile create the reference genome for every new analysis. It would be better to build the reference once and always keep it in one place. I will now be in /project2/gilad/briana/genome_anotation_data. To do this I will change the snakefile to point to this location and change the config file to point to the correct location.

#config file  

dir_genome: /project2/gilad/briana/genome_anotation_data/

#snakefile
dir_gen=config["dir_genome"]
dir_genome= dir_gen + "genome/" 

Picard tools

Picard tools will take in a flat file and a rRna file and will assess coverage at genomic regions. Before I make the rule I need to convert the gencode.v19.annotation.gtf to a flat file.

Update net-seq environment:

• ucsc-gtfToGenePred
  
• picard >= 2.9.2

#convert gtf to genepred flat  
gtfToGenePred -genePredExt -ignoreGroupsWithoutExons gencode.v19.annotation.gtf gencode.v19.annotation.refFlat

Adapt rule from Seb’s Snakefile

#config file 
ref_flat: /project2/gilad/briana/genome_anotation_data/gencode.v19.annotation.refFlat
ribosome_intervals: /project2/gilad/briana/genome_anotation_data/hg19.rRNA.intervals  

#snakefile  
ref_flat = config["ref_flat"]
ribosome_int = config["ribosome_intervals"]


#rule 
rule collect_rna_metrics:
    input:
        bam = dir_sort + "{samples}-sort.bam"
        ref_flat = ref_flat,
        ribosome_int = ribosome_int
    output:
        rna_metrics = output + "{samples}_RNAmetrics.picard.txt"
    params:
        strand = "NONE",
        memory = "-Xmx12G"
    shell:
        "picard CollectRnaSeqMetrics {params.memory} I={input.bam} O={output.rna_metrics} REF_FLAT={input.ref_flat} STRAND={params.strand} RIBOSOMAL_INTERVALS={input.ribosome_int}"

Genome Coverage

The old snakefile has a genome_cov rule for genome cov of the deduplicated files. I want to do this for the sort and dedup files. I also will change the rule so it is -d -5.

rule genome_cov_dedup:
    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 -d -5 -ibam {input.bamdedup} -g {input.genome} > {output}"
    
    
rule genome_cov_sort:
    input:
      bamsort = dir_sort + "{samples}-sort.bam"
      genome = dir_genome + ensembl_genome + ".reads"  
    output: 
      dir_cov + "{samples}-sort.cov.bed"
    shell: "bedtools genomecov -d -5 -ibam {input.bamsort} -g {input.genome} > {output}"
   
   
#in rule all:  
expand(dir_cov + "{samples}-sort.dedup.cov.bed", samples=samples),
expand(dir_cov + "{samples}-sort.cov.bed", samples=samples),

To test: snakemake -np
to run: nohup scripts/submit-snakemake.sh

awk '{ print $12 "\t" $1 "\t" $2 "\t" $3 "\t" $4 "\t" $5 "\t" $6 "\t" $7 "\t" $8 "\t" $9 "\t" $10 }' gencode.v19.annotation.refFlat > gencode.v19.annotation2.refFlat

picard_all=read.csv("../data/all_RNAmetrics.picard.none.csv", header=TRUE)

library(tidyr)
library(reshape2)

Warning: package 'reshape2' was built under R version 3.4.3

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

library(ggplot2)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(stringr)
#make libarary a factor
picard_all$line <= factor(picard_all$line)

Warning in Ops.factor(picard_all$line, factor(picard_all$line)): '<=' not
meaningful for factors

[1] NA NA NA NA NA NA NA NA

long_picard=melt(picard_all, id.vars=("line"))

pct_df=filter(long_picard,variable=='PCT_RIBOSOMAL_BASES' | variable=='PCT_CODING_BASES' | variable=='PCT_UTR_BASES' | variable=='PCT_INTRONIC_BASES' | variable=='PCT_MRNA_BASES')

pct_df$line = factor(pct_df$line)
ggplot(pct_df,aes(line, value, fill=variable )) + geom_bar(stat = "identity",position="dodge") + ggtitle("Picard read distribution") + labs(y="Percent")