feat: add TIN score calculation (#18)

Author: balajtimate

README.md

@@ -66,6 +66,7 @@ Outputs:
 - `analysis`: to only run the postprocessing part of the workflow (quantification of IPA usage and intron retention) - requires `input_bam`
 - `tectool`: to only run the IPA usage quantification, using TECtool - requires `input_bam`
 - `intron`: to only run the intron retention quantification subworkflow - requires `input_bam`
+- `tin_score`: to only run the TIN score calculation subworkflow - requires `input_bam`
 
 In the case of `full` and `preprocessing` modes, the `input_fastq` is required, using a wildcard character, e.g.: `--input_fastq='test_data\*{1,2}.fastq'`
 

conf/envs/slurm.config

@@ -28,11 +28,17 @@ process {
         memory = { 8.GB * task.attempt }
     }
     withLabel: TECtool {
-        clusterOptions = "--qos=6hours"
+        clusterOptions = { task.attempt > 1 ? "--qos=1day" : "--qos=6hours" }
         cpus = 8
-        time = { 1.h * task.attempt }
+        time = { 6.h * task.attempt }
         memory = { 64.GB * task.attempt }
     }
+    withLabel: calculate_tin {
+        clusterOptions = "--qos=6hours"
+        cpus = 8
+        time = { 3.h * task.attempt }
+        memory = 32.GB
+    }
     withLabel: salmon {
         clusterOptions = "--qos=6hours"
         cpus = 8
@@ -47,7 +53,9 @@ process {
     // Specific requirements
 
     withName: SAMTOOLS_GET_LOW_DUP_READS {
+        clusterOptions = "--qos=6hours"
+        cpus = 8
         memory = { 40.GB * task.attempt }
-        time = 30.m
+        time = { 2.h * task.attempt }
     }
 }

conf/envs/slurm_med.config

@@ -3,10 +3,10 @@ process {
     // standard settings if not specified otherwise
 
     executor = 'slurm'
-    clusterOptions = "--partition=dynamic-8cores-16g"
-    errorStrategy = { task.exitStatus in 137..255 ? 'retry' : 'terminate' }
+    clusterOptions = "--partition=dynamic-8cores-16g,dynamic-16cores-32g,dynamic-16cores-64g,dynamic-16cores-128g"
+    errorStrategy = { task.exitStatus in 137..140 ? 'retry' : 'terminate' }
     maxRetries = 2
-    time = { 20.m * task.attempt }
+    time = { 30.m * task.attempt }
     memory = { 4.GB * task.attempt }
     cpus = 4
     queueSize = 10

install/environment.yml

@@ -16,5 +16,7 @@ dependencies:
   - bedops
   - gffread
   - gtfparse
-  - pip:
-    - git+https://github.com/balajtimate/TECtool@master
+  - zgtf
+  - tin-score-calculation
+  # - pip:
+  #   - git+https://github.com/balajtimate/TECtool@master

main.nf

@@ -30,6 +30,8 @@ include { SALMON_TRANSCRIPTOME } from './modules/salmon.nf'
 include { SALMON_INDEX } from './modules/salmon.nf'
 include { SALMON_QUANTIFY } from './modules/salmon.nf'
 include { INTRON_RETENTION } from './modules/intron_retention.nf'
+include { TIN_GTF2BED } from './modules/tin_score.nf'
+include { CALCULATE_TIN_SCORES } from './modules/tin_score.nf'
 
 genome_index_ch = Channel.fromPath(params.genome_index, checkIfExists: true).collect()
 annotation_gtf_ch = Channel.fromPath(params.annotation_gtf, checkIfExists: true).collect()
@@ -55,6 +57,22 @@ workflow preprocessing {
         bam_low_dupl_tupl
 }
 
+// Subworkflow for TECtool analysis and downstream steps
+workflow calculate_tin_scores {
+    take:
+        input_bam
+
+    main:
+        // Convert BAM to 2 FASTQ files
+        TIN_GTF2BED(annotation_gtf_ch)
+        transcripts_bed12 = TIN_GTF2BED.out.transcripts_bed12
+        CALCULATE_TIN_SCORES(input_bam, transcripts_bed12)
+        tin_scores_tsv = CALCULATE_TIN_SCORES.out.tin_scores_tsv
+
+    emit:
+        tin_scores_tsv
+}
+
 // Subworkflow for TECtool analysis and downstream steps
 workflow tectool_analysis {
     take:
@@ -113,6 +131,7 @@ workflow {
             preprocessing(input_fastq_ch)
             tectool_analysis(preprocessing.out.bam_low_dupl_tupl)
             intron_retention(preprocessing.out.bam_low_dupl_tupl)
+            calculate_tin_scores(preprocessing.out.bam_low_dupl_tupl)
         }
     }
     if (params.run_mode == 'preprocessing') {
@@ -134,6 +153,12 @@ workflow {
             tectool_analysis(input_bam_ch)
         }
     }
+    if (params.run_mode == 'tin_score') {
+        input_bam_ch = Channel.fromPath(params.input_bam, checkIfExists: true).map { bam_path -> tuple(bam_path.baseName, bam_path) }
+        input_bam_ch.each {
+            calculate_tin_scores(input_bam_ch)
+        }
+    }
     if (params.run_mode == 'intron') {
         input_bam_ch = Channel.fromPath(params.input_bam, checkIfExists: true).map { bam_path -> tuple(bam_path.baseName, bam_path) }
         input_bam_ch.each {

modules/salmon.nf

@@ -33,20 +33,22 @@ process SALMON_INDEX {
     tag { library }
 
     // publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: '*'
-    
+    publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: '*.log'
+
     input:
     tuple val(library), file(fasta)
 
     output:
-    tuple val(library), path('*'), emit: salmon_index
+    tuple val(library), path('*_transcripts_index'), emit: salmon_index
+    path '*.log', emit: log
 
     script:
     """
     salmon index \
         --transcripts ${fasta} \
         --index ${library}_transcripts_index \
         --keepDuplicates \
-        --threads ${params.threads_pe}
+        --threads ${params.threads_pe} &> ${library}_salmon_index.log
     """
 }
 
@@ -57,14 +59,16 @@ process SALMON_QUANTIFY {
     tag { library }
 
     publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: '*_quant.tsv'
-    
+    publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: '*.log'
+
     input:
     tuple val(library), path(index)
     tuple val(library_1), path(fastq_1)
     tuple val(library_2), path(fastq_2)
 
     output:
     tuple val(library), path('*_quant.tsv'), emit: salmon_counts
+    path '*.log', emit: log
 
     script:
     """
@@ -75,7 +79,7 @@ process SALMON_QUANTIFY {
         --validateMappings \
         --seqBias \
         --output ${library}_transcript_quant \
-        --threads ${params.threads_pe}
+        --threads ${params.threads_pe} &> ${library}_salmon_quant.log
     mv ${library}_transcript_quant/quant.sf ${library}_quant.tsv
     """
 }

modules/samtools.nf

@@ -2,26 +2,6 @@
 
 nextflow.enable.dsl=2
 
-// process SAMTOOLS_INDEX {
-// 
-//     label "samtools"
-//     
-//     tag { library }
-// 
-//     // publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*.bai"
-// 
-//     input:
-//     tuple val(library), path(bam)
-// 
-//     output:
-//     tuple val(library), path('*.bai'), emit: index
-// 
-//     script:
-//     """
-//     samtools index -@ ${params.threads_se} -M ${bam}
-//     """
-// }
-
 process SAMTOOLS_GET_UNIQUE_MAPPERS {
 
     label "samtools"
@@ -81,8 +61,8 @@ process SAMTOOLS_BAM2FASTQ {
 
     tag { library }
 
-    // publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_1.fastq"
-    // publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_2.fastq"
+    publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_1.fastq"
+    publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_2.fastq"
     publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: '*.log'
 
     input:
@@ -91,6 +71,7 @@ process SAMTOOLS_BAM2FASTQ {
     output:  
     tuple val("${library}_1"), path("${library}_1.fastq"), emit: fastq1_tuple
     tuple val("${library}_2"), path("${library}_2.fastq"), emit: fastq2_tuple
+    path '*.log', emit: log
 
     script:
     """

modules/tectool.nf

@@ -9,6 +9,7 @@ process TECTOOL {
     tag { library_split }
 
     publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*/*.tsv"
+    publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*.bai"
     publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_enriched_annotation.gtf"
     publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: '*.log'
 
@@ -22,6 +23,7 @@ process TECTOOL {
     output:
     tuple val(library_split), path('*.gtf'), emit: enriched_gtf
     path '*/*.tsv', emit: tsv
+    path '*.log', emit: log
 
 
     script:
@@ -33,7 +35,8 @@ process TECTOOL {
         --bam ${bam_split} \
         --genome ${genome_fa} \
         --num_cores ${params.threads_se} \
-        --output_dir ${library_split}_tectool &> ${library_split}_tectool.log
+        --output_dir ${library_split}_tectool \
+        &> ${library_split}_tectool.log
     mv ${library_split}_tectool/enriched_annotation.gtf ${library_split}_enriched_annotation.gtf  
     """
 }
@@ -42,9 +45,10 @@ process TECTOOL_MERGE {
 
     label 'bedtools'
 
-    tag { library }
+    tag { "${library}: ${library_1}, ${library_2}" }
 
     publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*_merged.gtf"
+    publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: '*.log'
 
     input:
     tuple val(library), file(bam)
@@ -53,13 +57,15 @@ process TECTOOL_MERGE {
 
     output:
     tuple val(library), path('*_merged.gtf'), emit: merged_gtf
+    path '*.log', emit: log
 
     script:
     """
     echo -e "${gtf_files_1}\n${gtf_files_2}" > ${library}_tectool_annotation_files.tsv
     tectool_add_novel_transcripts_to_gtf_file \
         --list_of_gtf_files ${library}_tectool_annotation_files.tsv \
-        --out-dir ${library}_tectool_merged_annotations
+        --out-dir ${library}_tectool_merged_annotations \
+        &> ${library}_tectool_merged.log
     mv ${library}_tectool_merged_annotations/merged_annotation.gtf ${library}_tectool_annotation_merged.gtf
     """
 }

modules/tin_score.nf

@@ -0,0 +1,47 @@
+#!/usr/bin/env nextflow
+
+nextflow.enable.dsl=2
+
+process TIN_GTF2BED {
+
+    label 'gtf2bed'
+
+    publishDir "${params.out_dir}", mode: 'copy', pattern: "*.bed"
+    publishDir "${params.log_dir}", mode: 'copy', pattern: "*.log"
+
+    input:
+    path annotation_gtf
+
+    output:
+    path ('*.bed'), emit: transcripts_bed12
+
+    script:
+    """
+    sed 's/transcript_type/transcript_biotype/g' ${annotation_gtf} > ${annotation_gtf}.biotype
+    sort -k1,1 -k4,4n -k5,5nr ${annotation_gtf}.biotype > ${annotation_gtf}.sorted
+    gtf2bed12 --gtf ${annotation_gtf}.sorted --bed12 full_transcripts_protein_coding.bed 2> ${annotation_gtf}_gtf2bed.log
+    """
+}
+
+process CALCULATE_TIN_SCORES {
+
+    label 'calculate_tin'
+
+    tag { library }
+
+    publishDir "${params.out_dir}/${library}_results", mode: 'copy', pattern: "*.tsv"
+    publishDir "${params.log_dir}/${library}_logs", mode: 'copy', pattern: "*.log"
+
+    input:
+    tuple val(library), path(bam)
+    path transcripts_bed12
+
+    output:
+    path ('*.tsv'), emit: tin_scores_tsv
+
+    script:
+    """
+    samtools index -@ ${params.threads_se} -M ${bam}
+    calculate-tin.py -i ${bam} -r ${transcripts_bed12} --names ${library} -p ${params.threads_se} > ${library}_TIN_score.tsv 2> ${library}_tin_scores.log
+    """
+}

nextflow.config

@@ -47,7 +47,6 @@ profiles {
   }
   conda {
     conda.enabled = true
-    conda.channels = 'bioconda,conda-forge,defaults'
-    process.conda = 'install/environment.yml'
+    process.conda = "${HOME}/miniconda3/envs/ipa-immune"
   }
 }