main.nf

#!/usr/bin/env nextflow
nextflow.enable.dsl=2


process get_images {
  stageInMode 'symlink'
  stageOutMode 'move'

  script:
    """

    if [[ "${params.containers}" == "singularity" ]] ; 

      then

        cd ${params.image_folder}

        if [[ ! -f rnaseq.python-3.8-4.sif ]] ;
          then
            singularity pull rnaseq.python-3.8-4.sif docker://index.docker.io/mpgagebioinformatics/rnaseq.python:3.8-4
        fi

        if [[ ! -f crisprcleanr-3.0.0.sif ]] ;
          then
            singularity pull crisprcleanr-3.0.0.sif docker://index.docker.io/mpgagebioinformatics/crisprcleanr:3.0.0
        fi

    fi


    if [[ "${params.containers}" == "docker" ]] ; 

      then
        docker pull mpgagebioinformatics/rnaseq.python:3.8-4
        docker pull mpgagebioinformatics/crisprcleanr:3.0.0
    fi

    """

}

process lib_file {
  stageInMode 'symlink'
  stageOutMode 'move'

  when:
    ( ! file("${params.project_folder}/library_cleanR.tsv").exists() )

  script:
  """
#!/usr/local/bin/python
import pandas as pd
import numpy as np
import os

print("Starting")
import sys
sys.stdout.flush()

EXC=pd.ExcelFile("${params.reference_file}", engine="openpyxl")

# process library information
df=EXC.parse("library")
#df=df.dropna()
df=df.rename(columns={"gene_ID":"GENES", "UID":"CODE"})
df=df[["CODE","GENES","EXON","CHRM","STRAND","STARTpos","ENDpos","seq"]]
#df["STARTpos"]=df["STARTpos"].astype(int).astype(str)
#df["ENDpos"]=df["ENDpos"].astype(int).astype(str)
df["STARTpos"]=[str(int(s)) if str(s) != str(np.nan) else np.nan for s in df["STARTpos"]]
df["ENDpos"]=[str(int(s)) if str(s) != str(np.nan) else np.nan for s in df["ENDpos"]]
df.to_csv("${params.project_folder}"+"/library_cleanR.tsv", sep="\\t", index=None)
  """
}

process subsetting_counts_file{
  stageInMode 'symlink'
  stageOutMode 'move'

  input:
    val label
    val control
    val treatment

  output:
    val label

  when:
    ( ! file("${params.cleanR_output}/${label}.counts.txt").exists() ) 

  script:
  """
#!/usr/local/bin/python
import pandas as pd
import os

counts=pd.read_csv("${params.cleanR_input_mageck}", sep="\\t")
counts=counts.rename(columns={"Gene":"gene"})

control="${control}"
control_cols=control.split(",")

treatment="${treatment}"
treartment_cols=treatment.split(",")

counts_sub=counts[["sgRNA","gene"]+control_cols+treartment_cols]
counts_sub.to_csv("${params.cleanR_output}/${label}.counts.txt", sep="\\t", index=None, header=True)

  """

}

process cleanR_pipe {
  stageInMode 'symlink'
  stageOutMode 'move'

  input:
    val label
    val control
    val treatment
    val out

  when:
    ( ! file("${params.cleanR_output}/${label}.gene_signatures.tsv").exists() )

  script:
  """
#!/usr/bin/env Rscript
library(CRISPRcleanR)

setwd("${params.cleanR_output}/")

###library_file <- "${params.cleanR_lib_file}"

library_file <- "${params.project_folder}/library_cleanR.tsv"
fn <- "${params.cleanR_output}/${label}.counts.txt"

#############################################
## Pipeline
#############################################

#"FUN": ["ccr.getLibrary"]
#"desc": ["Load Library Annotation"]
lib <- ccr.getLibrary(library_file = library_file, library_builtin = NULL)

#"FUN": ["ccr.getCounts"]
#"desc": ["Load Count data"]
counts <- ccr.getCounts(file_counts = fn, 
                    verbose = TRUE, 
                    files_FASTQ_controls=NULL, 
                    files_FASTQ_samples=NULL,
                    files_BAM_controls=NULL,
                    files_BAM_samples=NULL)

write.table(counts, file = "${params.cleanR_output}/${label}.raw_counts.tsv", sep="\\t", quote = FALSE, row.names = FALSE)

#"FUN": ["ccr.checkCounts"]
#"desc": ["Check Library/Count data"]
####ccr.checkCounts(counts = counts, libraryAnnotation = lib, ncontrols = 2, min_reads = 30)
#### add if statement here that if false do not proceed

#"FUN": ["ccr.NormfoldChanges"]
#"desc": ["Run Count normalization"]
normCountsAndFCs <- ccr.NormfoldChanges(filename = fn, 
                            saveToFig=TRUE,
                            display = TRUE,
                            method = 'ScalingByTotalReads', ## One of (ScalingByTotalReads, MedRatios)
                            ncontrols = ${params.cleanR_control_reps},
                            libraryAnnotation = lib,
                            EXPname = '${label}',
                            min_reads = 30,
                            outdir = "${params.cleanR_output}/")

write.table(normCountsAndFCs\$norm_counts, file="${params.cleanR_output}/${label}.count_norm.tsv", sep="\\t", quote = FALSE, row.names = FALSE)
write.table(normCountsAndFCs\$logFCs, file="${params.cleanR_output}/${label}.logFCs.tsv", sep="\\t", quote = FALSE, row.names = FALSE)


#### CN correction

#"FUN": ["ccr.logFCs2chromPos"]
#"desc": ["Run sgRNA Sorting"]
gwSortedFCs <- ccr.logFCs2chromPos(foldchanges = normCountsAndFCs\$logFCs, libraryAnnotation = lib)

# "FUN": ["ccr.GWclean"]
#"desc": ["Run sgRNA Sorting"]
correctedFCs <- ccr.GWclean(gwSortedFCs = gwSortedFCs ,label='${label}',display=TRUE,
                        saveTO="${params.cleanR_output}/",
                        ignoredGenes=NULL,
                        min.ngenes=3,
                        alpha = 0.01,
                        nperm = 10000,
                        p.method ="hybrid", ## one of (perm, hybrid)
                        min.width=2,
                        kmax=25,
                        nmin=200, 
                        eta=0.05,
                        trim = 0.025,
                        undo.splits = "none",
                        undo.prune=0.05, 
                        undo.SD=3,
                        return.segments.unadj = TRUE,
                        return.segments.adj = TRUE)

write.table(correctedFCs\$segments, file="${params.cleanR_output}/${label}.segments.tsv", sep="\\t", quote = FALSE, row.names = FALSE)
write.table(correctedFCs\$segments_adj, file="${params.cleanR_output}/${label}.segments_adj.tsv", sep="\\t", quote = FALSE, row.names = FALSE)
write.table(correctedFCs\$corrected_logFCs, "${params.cleanR_output}/${label}.logFCs_adj.tsv", sep="\\t", quote = FALSE, row.names = FALSE)


#"FUN": ["ccr.correctCounts"]
#"desc": ["Correct counts"]
correctedCounts <- ccr.correctCounts(CL = '${label}',
                                normalised_counts = normCountsAndFCs\$norm_counts,
                                correctedFCs_and_segments = correctedFCs,
                                libraryAnnotation = lib,
                                minTargetedGenes=3,
                                OutDir="${params.cleanR_output}/",
                                ncontrols=${params.cleanR_control_reps})

write.table(correctedCounts, file="${params.cleanR_output}/${label}.counts_corrected.tsv", sep="\\t", quote = FALSE, row.names = FALSE)



#"FUN": ["ccr.sgRNAmeanFCs"]
#"desc": ["Get corrected FCs by probe"]
sgRNAmeanFCs <- ccr.sgRNAmeanFCs(foldchanges = correctedFCs\$corrected_logFCs) 

#"FUN": ["ccr.geneMeanFCs"]
#"desc": ["Get corrected FCs by gene"]
geneMeanFCs <- ccr.geneMeanFCs(sgRNA_FCprofile = sgRNAmeanFCs, libraryAnnotation = lib)

#"FUN": ["ccr.geneSummary"]
#"desc": ["Get corrected FCs by gene"]
geneSummary <- ccr.geneSummary(sgRNA_FCprofile = sgRNAmeanFCs,
                        libraryAnnotation = lib,
                        FDRth=0.05)
write.table(geneSummary, file="${params.cleanR_output}/${label}.gene_summary.tsv", sep="\\t", quote = FALSE, row.names = FALSE)

data("BAGEL_essential")
data("BAGEL_nonEssential")

BAGEL_essential_sgRNAs <- ccr.genes2sgRNAs(libraryAnnotation = lib, BAGEL_essential)
BAGEL_nonEssential_sgRNAs <- ccr.genes2sgRNAs(libraryAnnotation = lib,BAGEL_nonEssential)

#"FUN": ["ccr.ROC_Curve"]
#"desc": ["ROC curve at sgRNA level"]
pdf(paste0("${params.cleanR_output}/${label}.ROC_by_sgRNA.pdf"), width = 10, height = 10)
sgRNA_level_ROC <- ccr.ROC_Curve(FCsprofile = sgRNAmeanFCs, 
                            positives = BAGEL_essential_sgRNAs,
                            negatives = BAGEL_nonEssential_sgRNAs,
                            display = TRUE,
                            FDRth = 0.05,
                            expName = '${label}')
dev.off()

#"FUN": ["ccr.ROC_Curve"]
#"desc": ["ROC curve at gene level"]
pdf(paste0("${params.cleanR_output}/${label}.ROC_by_gene.pdf"), width = 10, height = 10)
gene_level_ROC <- ccr.ROC_Curve(FCsprofile = geneMeanFCs,
                        positives = BAGEL_essential,
                        negatives = BAGEL_nonEssential,
                        display = TRUE,
                        FDRth = 0.05,
                        expName = '${label}')
dev.off()                 

#"FUN": ["ccr.PrRc_Curve"]
#"desc": ["PrRc curve at sgRNA level"]
pdf(paste0("${params.cleanR_output}/${label}.PrRc_by_sgRNA.pdf"), width = 10, height = 10)
sgRNA_level_PrRc <- ccr.PrRc_Curve(FCsprofile = sgRNAmeanFCs,
                            positives = BAGEL_essential_sgRNAs,
                            negatives = BAGEL_nonEssential_sgRNAs,
                            display = TRUE,
                            FDRth = 0.05,
                            expName = "${label}")
dev.off()

#"FUN": ["ccr.PrRc_Curve"]
#"desc": ["PrRc curve at gene level"]
pdf(paste0("${params.cleanR_output}/${label}.PrRc_by_gene.pdf"), width = 10, height = 10)
gene_level_PrRc<-ccr.PrRc_Curve(FCsprofile = geneMeanFCs,
                            positives = BAGEL_essential,
                            negatives = BAGEL_nonEssential,
                            display = TRUE,
                            FDRth = 0.05,
                            expName = "${label}")
dev.off()


#"FUN": ["ccr.VisDepAndSig"]
#"desc": ["Visialize dependecy by gene signatures"]

data(EssGenes.ribosomalProteins)
data(EssGenes.DNA_REPLICATION_cons)
data(EssGenes.KEGG_rna_polymerase)
data(EssGenes.PROTEASOME_cons)
data(EssGenes.SPLICEOSOME_cons)

SIGNATURES <- list(Ribosomal_Proteins=EssGenes.ribosomalProteins,
            DNA_Replication = EssGenes.DNA_REPLICATION_cons,
            RNA_polymerase = EssGenes.KEGG_rna_polymerase,
            Proteasome = EssGenes.PROTEASOME_cons,
            Spliceosome = EssGenes.SPLICEOSOME_cons,
            CFE = BAGEL_essential,
            non_essential = BAGEL_nonEssential)

pdf(paste0("${params.cleanR_output}/${label}.gene_signatures.pdf"), width = 10, height = 10)
Recall_scores <- ccr.VisDepAndSig(FCsprofile = geneMeanFCs,
                            SIGNATURES = SIGNATURES,
                            TITLE = '${label}',
                            pIs = 6,
                            nIs = 7,
                            th=0.05,
                            plotFCprofile=TRUE)
dev.off()
write.table(as.data.frame(Recall_scores), file="${params.cleanR_output}/${label}.gene_signatures.tsv", sep="\\t", quote = FALSE, row.names = TRUE)



  """
}

workflow images {
  main:
    get_images()
}

workflow lib_file_formatting {
  if ( 'cleanR_output' in params.keySet() ) {
    lib_file()
  } 
}

workflow cleanR_workflow {
  if ( 'cleanR_output' in params.keySet() ) {

    if ( ! file("${params.cleanR_output}").isDirectory() ) {
      file("${params.cleanR_output}").mkdirs()
    }

    rows=Channel.fromPath("${params.samples_tsv}", checkIfExists:true).splitCsv(sep:';')
    rows=rows.filter{ ! file( "${params.cleanR_output}/${it[0]}.gene_signatures.tsv" ).exists() }
    
    label=rows.flatMap { n -> n[0] }

    control=rows.flatMap { n -> n[2] }
    control=control.map{ "$it".replace(".fastq.gz","") }
    treatment=rows.flatMap { n -> n[3] }
    treatment=treatment.map{ "$it".replace(".fastq.gz","") }
    
    subsetting_counts_file( label, control, treatment )
    cleanR_pipe( label, control, treatment, subsetting_counts_file.out.collect() )

  }
}