diff --git a/bin/awk_filter_reads.sh b/bin/awk_filter_reads.sh new file mode 100755 index 0000000..d50aa9e --- /dev/null +++ b/bin/awk_filter_reads.sh @@ -0,0 +1 @@ +awk 'BEGIN{OFS="\t"}{if($1 ~ /^\@/) {print($0)} else {$2=and($2,compl(2048)); print(substr($0,2))}}' diff --git a/bin/filter_five_end.pl b/bin/filter_five_end.pl new file mode 100644 index 0000000..6f9b4d5 --- /dev/null +++ b/bin/filter_five_end.pl @@ -0,0 +1,109 @@ +#!/usr/bin/perl +use strict; +use warnings; + +my $prev_id = ""; +my @five; +my @three; +my @unmap; +my @mid; +my @all; +my $counter = 0; + +while (){ + chomp; + if (/^@/){ + print $_."\n"; + next; + } + my ($id, $flag, $chr_from, $loc_from, $mapq, $cigar, $d1, $d2, $d3, $read, $read_qual, @rest) = split /\t/; + my $bin = reverse(dec2bin($flag)); + my @binary = split(//,$bin); + if ($prev_id ne $id && $prev_id ne ""){ + if ($counter == 1){ + if (@five == 1){ + print $five[0]."\n"; + } + else{ + my ($id_1, $flag_1, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1) = split /\t/, $all[0]; + my $bin_1 = reverse(dec2bin($flag_1)); + my @binary_1 = split(//,$bin_1); + $binary_1[2] = 1; + my $bin_1_new = reverse(join("",@binary_1)); + my $flag_1_new = bin2dec($bin_1_new); + print(join("\t",$id_1, $flag_1_new, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1)."\n"); + } + } + elsif ($counter == 2 && @five == 1){ + print $five[0]."\n"; + } + else{ + my ($id_1, $flag_1, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1) = split /\t/, $all[0]; + my $bin_1 = reverse(dec2bin($flag_1)); + my @binary_1 = split(//,$bin_1); + $binary_1[2] = 1; + my $bin_1_new = reverse(join("",@binary_1)); + my $flag_1_new = bin2dec($bin_1_new); + print(join("\t",$id_1, $flag_1_new, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1)."\n"); + } + + $counter = 0; + undef @unmap; + undef @five; + undef @three; + undef @mid; + undef @all; + } + + $counter++; + $prev_id = $id; + push @all,$_; + if ($binary[2]==1){ + push @unmap,$_; + } + elsif ($binary[4]==0 && $cigar =~ m/^[0-9]*M/ || $binary[4]==1 && $cigar =~ m/.*M$/){ + push @five, $_; + } + elsif ($binary[4]==1 && $cigar =~ m/^[0-9]*M/ || $binary[4]==0 && $cigar =~ m/.*M$/){ + push @three, $_; + } + elsif ($cigar =~ m/^[0-9]*[HS].*M.*[HS]$/){ + push @mid, $_; + } +} + +if ($counter == 1){ + if (@five == 1){ + print $five[0]."\n"; + } + else{ + my ($id_1, $flag_1, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1) = split /\t/, $all[0]; + my $bin_1 = reverse(dec2bin($flag_1)); + my @binary_1 = split(//,$bin_1); + $binary_1[2] = 1; + my $bin_1_new = reverse(join("",@binary_1)); + my $flag_1_new = bin2dec($bin_1_new); + print(join("\t",$id_1, $flag_1_new, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1)."\n"); + } +} +elsif ($counter == 2 && @five == 1){ + print $five[0]."\n"; +} +else{ + my ($id_1, $flag_1, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1) = split /\t/, $all[0]; + my $bin_1 = reverse(dec2bin($flag_1)); + my @binary_1 = split(//,$bin_1); + $binary_1[2] = 1; + my $bin_1_new = reverse(join("",@binary_1)); + my $flag_1_new = bin2dec($bin_1_new); + print(join("\t",$id_1, $flag_1_new, $chr_from_1, $loc_from_1, $mapq_1, $cigar_1, $d1_1, $d2_1, $d3_1, $read_1, $read_qual_1, @rest_1)."\n"); +} + +sub dec2bin { + my $str = unpack("B32", pack("N", shift)); + return $str; +} + +sub bin2dec { + return unpack("N", pack("B32", substr("0" x 32 . shift, -32))); +} diff --git a/bin/generate_cram_csv.sh b/bin/generate_cram_csv.sh new file mode 100755 index 0000000..9e56831 --- /dev/null +++ b/bin/generate_cram_csv.sh @@ -0,0 +1,99 @@ +#!/bin/bash + +# generate_cram_csv.sh +# ------------------- +# Generate a csv file describing the CRAM folder +# ><((((°> Y ><((((°> U ><((((°> M ><((((°> I ><((((°> +# Author = yy5 +# ><((((°> Y ><((((°> U ><((((°> M ><((((°> I ><((((°> + +# NOTE: chunk_size is the number of containers per chunk (not records/reads) + +# Function to process chunking of a CRAM file + +chunk_cram() { + local cram=$1 + local chunkn=$2 + local outcsv=$3 + local crai=$4 + local chunk_size=$5 + + local rgline=$(samtools view -H "${realcram}" | grep "@RG" | sed 's/\t/\\t/g' | sed "s/'//g") + local ncontainers=$(zcat "${realcrai}" | wc -l) + local base=$(basename "${realcram}" .cram) + + if [ $chunk_size -gt $ncontainers ]; then + chunk_size=$((ncontainers - 1)) + fi + local from=0 + local to=$((chunk_size - 1)) + + while [ $to -lt $ncontainers ]; do + #echo "chunk $chunkn: $from - $to" + echo "${realcram},${realcrai},${from},${to},${base},${chunkn},${rgline}" >> $outcsv + from=$((to + 1)) + to=$((to + chunk_size)) + ((chunkn++)) + done + + # Catch any remainder + if [ $from -lt $ncontainers ]; then + to=$((ncontainers - 1)) + #echo "chunk $chunkn: $from - $to" + echo "${realcram},${realcrai},${from},${to},${base},${chunkn},${rgline}" >> $outcsv + ((chunkn++)) + fi +} + +# Function to process a CRAM file +process_cram_file() { + local cram=$1 + local chunkn=$2 + local outcsv=$3 + local crai=$4 + local chunk_size=$5 + + local read_groups=$(samtools view -H "$cram" | grep '@RG' | awk '{for(i=1;i<=NF;i++){if($i ~ /^ID:/){print substr($i,4)}}}') + local num_read_groups=$(echo "$read_groups" | wc -w) + if [ "$num_read_groups" -gt 1 ]; then + # Multiple read groups: process each separately + for rg in $read_groups; do + local output_cram="$(basename "${cram%.cram}")_output_${rg}.cram" + samtools view -h -r "$rg" -o "$output_cram" "$cram" + #chunkn=$(chunk_cram "$output_cram" "$chunkn" "$outcsv" "$crai" "$chunk_size") + chunk_cram "$output_cram" "$chunkn" "$outcsv" "$crai" "$chunk_size" + done + else + # Single read group or no read groups + #chunkn=$(chunk_cram "$cram" "$chunkn" "$outcsv" "$crai" "$chunk_size") + chunk_cram "$cram" "$chunkn" "$outcsv" "$crai" "$chunk_size" + fi +} + +# /\_/\ /\_/\ +# ( o.o ) main ( o.o ) +# > ^ < > ^ < + +# Check if cram_path is provided +if [ $# -lt 3 ]; then + echo "Usage: $0 " + exit 1 +fi + +cram=$1 +outcsv=$2 +crai=$3 +if [ -z "$4" ]; then + chunk_size=10000 +else + chunk_size=$4 +fi +chunkn=0 + +# Operates on a single CRAM file +realcram=$(readlink -f $cram) +realcrai=$(readlink -f $crai) +if [ -f "$outcsv" ]; then + rm "$outcsv" +fi +process_cram_file $realcram $chunkn $outcsv $realcrai $chunk_size diff --git a/bin/grep_pg.sh b/bin/grep_pg.sh new file mode 100755 index 0000000..680b5ec --- /dev/null +++ b/bin/grep_pg.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +# grep_pg.sh +# ------------------- +# A shell script to exclude pg lines and label read 1 and read 2 from cram containers +# +# ------------------- +# Author = yy5 + +grep -v "^\@PG" | awk '{if($1 ~ /^\@/) {print($0)} else {if(and($2,64)>0) {print(1$0)} else {print(2$0)}}}' diff --git a/conf/base.config b/conf/base.config index 10ab224..6f3dab3 100644 --- a/conf/base.config +++ b/conf/base.config @@ -57,6 +57,12 @@ process { time = { check_max( 2.h * Math.ceil( meta.read_count / 100000000 ) * task.attempt / log_increase_cpus(2, 6*task.attempt, 1, 2), 'time' ) } } + withName: SAMTOOLS_INDEX { + cpus = { log_increase_cpus(2, 6*task.attempt, 1, 2) } + memory = { check_max( 4.GB + 850.MB * log_increase_cpus(2, 6*task.attempt, 1, 2) * task.attempt + 0.6.GB * Math.ceil( meta.read_count / 100000000 ), 'memory' ) } + time = { check_max( 2.h * Math.ceil( meta.read_count / 100000000 ) * task.attempt / log_increase_cpus(2, 6*task.attempt, 1, 2), 'time' ) } + } + withName: BLAST_BLASTN { time = { check_max( 2.hour * Math.ceil( meta.read_count / 1000000 ) * task.attempt, 'time' ) } memory = { check_max( 100.MB + 20.MB * Math.ceil( meta.read_count / 1000000 ) * task.attempt, 'memory' ) } @@ -93,6 +99,21 @@ process { time = { check_max( 1.h * Math.ceil( meta.read_count / 1000000 ) * task.attempt, 'time' ) } } + withName: CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT { + cpus = { check_max( 16 * 1 , 'cpus' ) } + memory = { check_max( 1.GB * ( reference.size() < 2e9 ? 50 : Math.ceil( ( reference.size() / 1e+9 ) * 20 ) * Math.ceil( task.attempt * 1 ) ) , 'memory') } + } + + withName: CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT { + cpus = { check_max( 16 * 1 , 'cpus' ) } + memory = { check_max( 1.GB * ( reference.size() < 2e9 ? 50 : Math.ceil( ( reference.size() / 1e+9 ) * 20 ) * Math.ceil( task.attempt * 1 ) ) , 'memory') } + } + + withName: MINIMAP2_INDEX { + cpus = { check_max( 2 * task.attempt, 'cpus' ) } + memory = { check_max( 1.GB * Math.ceil( 30 * fasta.size() / 1e+9 ) * task.attempt, 'memory' ) } + } + withName: CRUMBLE { // No correlation between memory usage and the number of reads or the genome size. // Most genomes seem happy with 1 GB, then some with 2 GB, then some with 5 GB. diff --git a/conf/modules.config b/conf/modules.config index d05fe4b..28116bd 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -12,6 +12,7 @@ process { withName: SAMTOOLS_FASTQ { + beforeScript = { "export REF_PATH=spoof"} ext.args = '-F 0x200 -nt' } @@ -28,6 +29,7 @@ process { } withName: SAMTOOLS_MERGE { + beforeScript = { "export REF_PATH=spoof"} ext.args = { "-c -p" } ext.prefix = { "${meta.id}.merge" } } @@ -39,6 +41,7 @@ process { } withName: SAMTOOLS_COLLATETOFASTA { + beforeScript = { "export REF_PATH=spoof"} ext.args = { (params.use_work_dir_as_temp ? "-T." : "") } } @@ -47,9 +50,58 @@ process { } withName: SAMTOOLS_CONVERT { + beforeScript = { "export REF_PATH=spoof"} ext.args = "-be '[rq]>=0.99' -x fi -x fp -x ri -x rp --write-index" } + withName: CONVERT_CRAM { + ext.args = "--output-fmt cram" + } + + withName: '.*:.*:ALIGN_HIC:BWAMEM2_MEM' { + ext.args = { "-5SPCp -R ${meta.read_group}" } + } + + withName: '.*:.*:ALIGN_ILLUMINA:BWAMEM2_MEM' { + ext.args = { "-p -R ${meta.read_group}" } + } + + withName: ".*:ALIGN_ILLUMINA:.*:CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT" { + ext.args = "" + ext.args1 = { "-F 0x200 -nt" } + ext.args2 = { "-p -R '${rglines}'" } + ext.args3 = "-mpu" + ext.args4 = { "--write-index -l1" } + } + + withName: ".*:ALIGN_ILLUMINA:.*:CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT" { + ext.args = "" + ext.args1 = { "-F 0x200 -nt" } + ext.args2 = { "-ax sr" } + ext.args3 = "-mpu" + ext.args4 = { "--write-index -l1" } + } + + withName: ".*:ALIGN_HIC:.*:CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT" { + ext.args = "" + ext.args1 = { "-F 0x200 -nt" } + ext.args2 = { "-5SPCp -R '${rglines}'" } + ext.args3 = "-mpu" + ext.args4 = { "--write-index -l1" } + } + + withName: ".*:ALIGN_HIC:.*:CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT" { + ext.args = "" + ext.args1 = { "" } + ext.args2 = { "-ax sr" } + ext.args3 = "-mpu" + ext.args4 = { "--write-index -l1" } + } + + withName: "MINIMAP2_INDEX" { + ext.args = { "${fasta.size() > 2.5e9 ? (" -I " + Math.ceil(fasta.size()/1e9)+"G") : ""} "} + } + // minimap2 2.24 can only work with genomes up to 4 Gbp. For larger genomes, add the -I option with the genome size in Gbp. // In fact, we can also use -I to *decrease* the memory requirements for smaller genomes // NOTE: minimap2 uses the decimal system ! 1G = 1,000,000,000 bp @@ -68,6 +120,7 @@ process { } withName: '.*:CONVERT_STATS:SAMTOOLS_CRAM' { + beforeScript = { "export REF_PATH=spoof"} ext.prefix = { "${fasta.baseName}.${meta.datatype}.${meta.id}" } ext.args = '--output-fmt cram --write-index' } @@ -85,7 +138,13 @@ process { ext.prefix = { "${bam.baseName}" } } + withName: SAMTOOLS_ADDREPLACERG { + ext.prefix = { "${input.baseName}_addRG" } + ext.args = { "-r ${meta.read_group} --no-PG" } + } + withName: SAMTOOLS_STATS { + beforeScript = { "export REF_PATH=spoof"} ext.prefix = { "${input.baseName}" } } @@ -95,6 +154,7 @@ process { } withName: '.*:CONVERT_STATS:SAMTOOLS_.*' { + beforeScript = { "export REF_PATH=spoof"} publishDir = [ path: { "${params.outdir}/read_mapping/${meta.datatype}" }, mode: params.publish_dir_mode, diff --git a/conf/test.config b/conf/test.config index d78e008..add9178 100644 --- a/conf/test.config +++ b/conf/test.config @@ -25,6 +25,9 @@ params { // Output directory outdir = "${projectDir}/results" + // Aligner + short_aligner = "minimap2" + // Fasta references fasta = "https://tolit.cog.sanger.ac.uk/test-data/Meles_meles/assembly/release/mMelMel3.1_paternal_haplotype/GCA_922984935.2.subset.fasta.gz" header = "https://tolit.cog.sanger.ac.uk/test-data/Meles_meles/assembly/release/mMelMel3.1_paternal_haplotype/GCA_922984935.2.subset.header.sam" diff --git a/modules.json b/modules.json index 522bb30..ea8b364 100644 --- a/modules.json +++ b/modules.json @@ -45,6 +45,11 @@ "git_sha": "a33ef9475558c6b8da08c5f522ddaca1ec810306", "installed_by": ["modules"] }, + "minimap2/index": { + "branch": "master", + "git_sha": "72e277acfd9e61a9f1368eafb4a9e83f5bcaa9f5", + "installed_by": ["modules"] + }, "samtools/faidx": { "branch": "master", "git_sha": "04fbbc7c43cebc0b95d5b126f6d9fe4effa33519", @@ -65,6 +70,11 @@ "git_sha": "46eca555142d6e597729fcb682adcc791796f514", "installed_by": ["modules"] }, + "samtools/index": { + "branch": "master", + "git_sha": "46eca555142d6e597729fcb682adcc791796f514", + "installed_by": ["modules"] + }, "samtools/merge": { "branch": "master", "git_sha": "04fbbc7c43cebc0b95d5b126f6d9fe4effa33519", diff --git a/modules/local/cram_filter_align_bwamem2_fixmate_sort.nf b/modules/local/cram_filter_align_bwamem2_fixmate_sort.nf new file mode 100644 index 0000000..667e49b --- /dev/null +++ b/modules/local/cram_filter_align_bwamem2_fixmate_sort.nf @@ -0,0 +1,60 @@ +process CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT { + tag "$meta.chunk_id" + label "process_high" + + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/mulled-v2-1a6fe65bd6674daba65066aa796ed8f5e8b4687b:688e175eb0db54de17822ba7810cc9e20fa06dd5-0' : + 'biocontainers/mulled-v2-1a6fe65bd6674daba65066aa796ed8f5e8b4687b:688e175eb0db54de17822ba7810cc9e20fa06dd5-0' }" + + input: + tuple val(meta), path(cramfile), path(cramindex), val(from), val(to), val(base), val(chunkid), val(rglines) + tuple val(meta2), path(fasta) + tuple val(meta3), path(index) + + output: + tuple val(meta), path("*.bam"), emit: mappedbam + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def args1 = task.ext.args1 ?: '' + def args2 = task.ext.args2 ?: '' + def args3 = task.ext.args3 ?: '' + def args4 = task.ext.args4 ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + // Please be aware one of the tools here required mem = 28 * reference size!!! + + """ + INDEX=`find -L ./ -name "*.amb" | sed 's/\\.amb\$//'` + + cram_filter -n ${from}-${to} ${cramfile} - | \\ + samtools fastq ${args1} | \\ + bwa-mem2 mem \$INDEX -t ${task.cpus} ${args2} - | \\ + samtools fixmate ${args3} - - | \\ + samtools sort ${args4} -@${task.cpus} -T ${base}_${chunkid}_sort_tmp -o ${prefix}_${base}_${chunkid}_mem.bam - + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + bwa-mem2: \$(bwa-mem2 --version | sed 's/bwa-mem2 //g') + END_VERSIONS + """ + // temp removal staden_io_lib: \$(echo \$(staden_io_lib --version 2>&1) | sed 's/^.*staden_io_lib //; s/Using.*\$//') CAUSES ERROR + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + def base = "45022_3#2" + def chunkid = "1" + """ + touch ${prefix}_${base}_${chunkid}_mem.bam + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + bwamem2: \$(echo \$(bwa-mem2 version 2>&1) | sed 's/.* //') + END_VERSIONS + """ +} diff --git a/modules/local/cram_filter_minimap2_filter5end_fixmate_sort.nf b/modules/local/cram_filter_minimap2_filter5end_fixmate_sort.nf new file mode 100644 index 0000000..ceafb86 --- /dev/null +++ b/modules/local/cram_filter_minimap2_filter5end_fixmate_sort.nf @@ -0,0 +1,59 @@ +process CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT { + tag "$meta.chunk_id" + label "process_high" + + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/mulled-v2-1a6fe65bd6674daba65066aa796ed8f5e8b4687b:688e175eb0db54de17822ba7810cc9e20fa06dd5-0' : + 'biocontainers/mulled-v2-1a6fe65bd6674daba65066aa796ed8f5e8b4687b:688e175eb0db54de17822ba7810cc9e20fa06dd5-0' }" + + input: + tuple val(meta), path(cramfile), path(cramindex), val(from), val(to), val(base), val(chunkid), val(rglines), val(ref), path(reference) + + output: + tuple val(meta), path("*.bam"), emit: mappedbam + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + def args1 = task.ext.args1 ?: '' + def args2 = task.ext.args2 ?: '' + def args3 = task.ext.args3 ?: '' + def args4 = task.ext.args4 ?: '' + def prefix = task.ext.prefix ?: "${meta.id}" + def VERSION = "1.15" // Staden_io versions break the pipeline + """ + cram_filter -n ${from}-${to} ${cramfile} - | \\ + samtools fastq ${args1} - | \\ + minimap2 -t${task.cpus} -R '${rglines}' ${args2} ${ref} - | \\ + ${projectDir}/bin/grep_pg.sh | \\ + perl ${projectDir}/bin/filter_five_end.pl | \\ + ${projectDir}/bin/awk_filter_reads.sh | \\ + samtools fixmate ${args3} - - | \\ + samtools sort ${args4} -@${task.cpus} -T ${base}_${chunkid}_sort_tmp -o ${prefix}_${base}_${chunkid}_mm.bam - + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + minimap2: \$(minimap2 --version | sed 's/minimap2 //g') + staden_io: $VERSION + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "${meta.id}" + def base = "45022_3#2" + def chunkid = "1" + """ + touch ${prefix}_${base}_${chunkid}_mm.bam + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + minimap2: \$(echo \$(minimap2 version 2>&1) | sed 's/.* //') + staden_io: $VERSION + END_VERSIONS + """ +} diff --git a/modules/local/generate_cram_csv.nf b/modules/local/generate_cram_csv.nf new file mode 100644 index 0000000..f3cfe96 --- /dev/null +++ b/modules/local/generate_cram_csv.nf @@ -0,0 +1,37 @@ +process GENERATE_CRAM_CSV { + tag "${meta.id}" + label 'process_single' + + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://quay.io/sanger-tol/cramfilter_bwamem2_minimap2_samtools_perl:0.001-c1' : + 'sanger-tol/cramfilter_bwamem2_minimap2_samtools_perl:0.001-c1' }" + + input: + tuple val(meta), path(crampath), path(crai) + + + output: + tuple val(meta), path('*.csv'), emit: csv + path "versions.yml", emit: versions + + script: + def prefix = task.ext.prefix ?: "${meta.id}" + """ + generate_cram_csv.sh $crampath ${prefix}_cram.csv $crai ${params.chunk_size} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + END_VERSIONS + """ + + stub: + """ + touch ${meta.id}.csv + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//' ) + END_VERSIONS + """ +} diff --git a/modules/local/samtools_addreplacerg.nf b/modules/local/samtools_addreplacerg.nf new file mode 100644 index 0000000..6515c3e --- /dev/null +++ b/modules/local/samtools_addreplacerg.nf @@ -0,0 +1,56 @@ +process SAMTOOLS_ADDREPLACERG { + tag "$meta.id" + label 'process_low' + + conda "bioconda::samtools=1.20" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/samtools:1.20--h50ea8bc_0' : + 'biocontainers/samtools:1.20--h50ea8bc_0' }" + + input: + tuple val(meta), path(input) + + output: + tuple val(meta), path("${prefix}.bam"), emit: bam, optional: true + tuple val(meta), path("${prefix}.cram"), emit: cram, optional: true + tuple val(meta), path("${prefix}.sam"), emit: sam, optional: true + path "versions.yml", emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: "${meta.id}" + file_type = input.getExtension() + if ("$input" == "${prefix}.${file_type}") error "Input and output names are the same, use \"task.ext.prefix\" to disambiguate!" + """ + samtools \\ + addreplacerg \\ + --threads $task.cpus \\ + $args \\ + -o ${prefix}.${file_type} \\ + $input + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//') + END_VERSIONS + """ + + stub: + def args = task.ext.args ?: '' + prefix = task.ext.prefix ?: "${meta.id}" + file_type = input.getExtension() + if ("$input" == "${prefix}.${file_type}") error "Input and output names are the same, use \"task.ext.prefix\" to disambiguate!" + + """ + touch ${prefix}.${file_type} + ${index} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//') + END_VERSIONS + """ +} diff --git a/modules/nf-core/minimap2/index/environment.yml b/modules/nf-core/minimap2/index/environment.yml new file mode 100644 index 0000000..8a912a1 --- /dev/null +++ b/modules/nf-core/minimap2/index/environment.yml @@ -0,0 +1,7 @@ +name: minimap2_index +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - bioconda::minimap2=2.28 diff --git a/modules/nf-core/minimap2/index/main.nf b/modules/nf-core/minimap2/index/main.nf new file mode 100644 index 0000000..3832021 --- /dev/null +++ b/modules/nf-core/minimap2/index/main.nf @@ -0,0 +1,44 @@ +process MINIMAP2_INDEX { + label 'process_low' + + // Note: the versions here need to match the versions used in minimap2/align + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/minimap2:2.28--he4a0461_0' : + 'biocontainers/minimap2:2.28--he4a0461_0' }" + + input: + tuple val(meta), path(fasta) + + output: + tuple val(meta), path("*.mmi"), emit: index + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + """ + minimap2 \\ + -t $task.cpus \\ + -d ${fasta.baseName}.mmi \\ + $args \\ + $fasta + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + minimap2: \$(minimap2 --version 2>&1) + END_VERSIONS + """ + + stub: + """ + touch ${fasta.baseName}.mmi + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + minimap2: \$(minimap2 --version 2>&1) + END_VERSIONS + """ +} diff --git a/modules/nf-core/minimap2/index/meta.yml b/modules/nf-core/minimap2/index/meta.yml new file mode 100644 index 0000000..1d29e3f --- /dev/null +++ b/modules/nf-core/minimap2/index/meta.yml @@ -0,0 +1,43 @@ +name: minimap2_index +description: Provides fasta index required by minimap2 alignment. +keywords: + - index + - fasta + - reference +tools: + - minimap2: + description: | + A versatile pairwise aligner for genomic and spliced nucleotide sequences. + homepage: https://github.com/lh3/minimap2 + documentation: https://github.com/lh3/minimap2#uguide + licence: ["MIT"] +input: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - fasta: + type: file + description: | + Reference database in FASTA format. +output: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - index: + type: file + description: Minimap2 fasta index. + pattern: "*.mmi" + - versions: + type: file + description: File containing software versions + pattern: "versions.yml" +authors: + - "@yuukiiwa" + - "@drpatelh" +maintainers: + - "@yuukiiwa" + - "@drpatelh" diff --git a/modules/nf-core/minimap2/index/tests/main.nf.test b/modules/nf-core/minimap2/index/tests/main.nf.test new file mode 100644 index 0000000..97840ff --- /dev/null +++ b/modules/nf-core/minimap2/index/tests/main.nf.test @@ -0,0 +1,32 @@ +nextflow_process { + + name "Test Process MINIMAP2_INDEX" + script "../main.nf" + process "MINIMAP2_INDEX" + + tag "modules" + tag "modules_nfcore" + tag "minimap2" + tag "minimap2/index" + + test("minimap2 index") { + + when { + process { + """ + input[0] = [ + [ id:'test' ], + file(params.modules_testdata_base_path + 'genomics/sarscov2/genome/genome.fasta', checkIfExists: true) + ] + """ + } + } + + then { + assert process.success + assert snapshot(process.out).match() + } + + } + +} \ No newline at end of file diff --git a/modules/nf-core/minimap2/index/tests/main.nf.test.snap b/modules/nf-core/minimap2/index/tests/main.nf.test.snap new file mode 100644 index 0000000..0b09882 --- /dev/null +++ b/modules/nf-core/minimap2/index/tests/main.nf.test.snap @@ -0,0 +1,68 @@ +{ + "Should run without failures": { + "content": [ + { + "0": [ + [ + { + "id": "test_ref" + }, + "genome.mmi:md5,72e450f12dc691e763c697463bdb1571" + ] + ], + "1": [ + "versions.yml:md5,0fced0ee8015e7f50b82566e3db8f7b0" + ], + "index": [ + [ + { + "id": "test_ref" + }, + "genome.mmi:md5,72e450f12dc691e763c697463bdb1571" + ] + ], + "versions": [ + "versions.yml:md5,0fced0ee8015e7f50b82566e3db8f7b0" + ] + } + ], + "meta": { + "nf-test": "0.8.4", + "nextflow": "23.10.1" + }, + "timestamp": "2024-03-18T11:46:30.000058092" + }, + "minimap2 index": { + "content": [ + { + "0": [ + [ + { + "id": "test" + }, + "genome.mmi:md5,72e450f12dc691e763c697463bdb1571" + ] + ], + "1": [ + "versions.yml:md5,2f8340380c6741e9261a284262a90bde" + ], + "index": [ + [ + { + "id": "test" + }, + "genome.mmi:md5,72e450f12dc691e763c697463bdb1571" + ] + ], + "versions": [ + "versions.yml:md5,2f8340380c6741e9261a284262a90bde" + ] + } + ], + "meta": { + "nf-test": "0.8.4", + "nextflow": "23.10.1" + }, + "timestamp": "2024-04-05T10:58:29.828187662" + } +} \ No newline at end of file diff --git a/modules/nf-core/minimap2/index/tests/tags.yml b/modules/nf-core/minimap2/index/tests/tags.yml new file mode 100644 index 0000000..e5ef8e1 --- /dev/null +++ b/modules/nf-core/minimap2/index/tests/tags.yml @@ -0,0 +1,2 @@ +minimap2/index: + - modules/nf-core/minimap2/index/** diff --git a/modules/nf-core/samtools/index/environment.yml b/modules/nf-core/samtools/index/environment.yml new file mode 100644 index 0000000..260d516 --- /dev/null +++ b/modules/nf-core/samtools/index/environment.yml @@ -0,0 +1,8 @@ +name: samtools_index +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - bioconda::samtools=1.20 + - bioconda::htslib=1.20 diff --git a/modules/nf-core/samtools/index/main.nf b/modules/nf-core/samtools/index/main.nf new file mode 100644 index 0000000..e002585 --- /dev/null +++ b/modules/nf-core/samtools/index/main.nf @@ -0,0 +1,49 @@ +process SAMTOOLS_INDEX { + tag "$meta.id" + label 'process_low' + + conda "${moduleDir}/environment.yml" + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/samtools:1.20--h50ea8bc_0' : + 'biocontainers/samtools:1.20--h50ea8bc_0' }" + + input: + tuple val(meta), path(input) + + output: + tuple val(meta), path("*.bai") , optional:true, emit: bai + tuple val(meta), path("*.csi") , optional:true, emit: csi + tuple val(meta), path("*.crai"), optional:true, emit: crai + path "versions.yml" , emit: versions + + when: + task.ext.when == null || task.ext.when + + script: + def args = task.ext.args ?: '' + """ + samtools \\ + index \\ + -@ ${task.cpus-1} \\ + $args \\ + $input + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//') + END_VERSIONS + """ + + stub: + def args = task.ext.args ?: '' + def extension = file(input).getExtension() == 'cram' ? + "crai" : args.contains("-c") ? "csi" : "bai" + """ + touch ${input}.${extension} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + samtools: \$(echo \$(samtools --version 2>&1) | sed 's/^.*samtools //; s/Using.*\$//') + END_VERSIONS + """ +} diff --git a/modules/nf-core/samtools/index/meta.yml b/modules/nf-core/samtools/index/meta.yml new file mode 100644 index 0000000..01a4ee0 --- /dev/null +++ b/modules/nf-core/samtools/index/meta.yml @@ -0,0 +1,57 @@ +name: samtools_index +description: Index SAM/BAM/CRAM file +keywords: + - index + - bam + - sam + - cram +tools: + - samtools: + description: | + SAMtools is a set of utilities for interacting with and post-processing + short DNA sequence read alignments in the SAM, BAM and CRAM formats, written by Heng Li. + These files are generated as output by short read aligners like BWA. + homepage: http://www.htslib.org/ + documentation: http://www.htslib.org/doc/samtools.html + doi: 10.1093/bioinformatics/btp352 + licence: ["MIT"] +input: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - bam: + type: file + description: BAM/CRAM/SAM file + pattern: "*.{bam,cram,sam}" +output: + - meta: + type: map + description: | + Groovy Map containing sample information + e.g. [ id:'test', single_end:false ] + - bai: + type: file + description: BAM/CRAM/SAM index file + pattern: "*.{bai,crai,sai}" + - crai: + type: file + description: BAM/CRAM/SAM index file + pattern: "*.{bai,crai,sai}" + - csi: + type: file + description: CSI index file + pattern: "*.{csi}" + - versions: + type: file + description: File containing software versions + pattern: "versions.yml" +authors: + - "@drpatelh" + - "@ewels" + - "@maxulysse" +maintainers: + - "@drpatelh" + - "@ewels" + - "@maxulysse" diff --git a/modules/nf-core/samtools/index/tests/csi.nextflow.config b/modules/nf-core/samtools/index/tests/csi.nextflow.config new file mode 100644 index 0000000..0ed260e --- /dev/null +++ b/modules/nf-core/samtools/index/tests/csi.nextflow.config @@ -0,0 +1,7 @@ +process { + + withName: SAMTOOLS_INDEX { + ext.args = '-c' + } + +} diff --git a/modules/nf-core/samtools/index/tests/main.nf.test b/modules/nf-core/samtools/index/tests/main.nf.test new file mode 100644 index 0000000..ca34fb5 --- /dev/null +++ b/modules/nf-core/samtools/index/tests/main.nf.test @@ -0,0 +1,140 @@ +nextflow_process { + + name "Test Process SAMTOOLS_INDEX" + script "../main.nf" + process "SAMTOOLS_INDEX" + tag "modules" + tag "modules_nfcore" + tag "samtools" + tag "samtools/index" + + test("bai") { + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("crai") { + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.recalibrated.sorted.cram', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("csi") { + config "./csi.nextflow.config" + + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot( + file(process.out.csi[0][1]).name, + process.out.versions + ).match() } + ) + } + } + + test("bai - stub") { + options "-stub" + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("crai - stub") { + options "-stub" + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/homo_sapiens/illumina/cram/test.paired_end.recalibrated.sorted.cram', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } + + test("csi - stub") { + options "-stub" + config "./csi.nextflow.config" + + when { + process { + """ + input[0] = Channel.of([ + [ id:'test', single_end:false ], // meta map + file(params.modules_testdata_base_path + 'genomics/sarscov2/illumina/bam/test.paired_end.sorted.bam', checkIfExists: true) + ]) + """ + } + } + + then { + assertAll ( + { assert process.success }, + { assert snapshot(process.out).match() } + ) + } + } +} diff --git a/modules/nf-core/samtools/index/tests/main.nf.test.snap b/modules/nf-core/samtools/index/tests/main.nf.test.snap new file mode 100644 index 0000000..799d199 --- /dev/null +++ b/modules/nf-core/samtools/index/tests/main.nf.test.snap @@ -0,0 +1,250 @@ +{ + "csi - stub": { + "content": [ + { + "0": [ + + ], + "1": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.csi:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "2": [ + + ], + "3": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ], + "bai": [ + + ], + "crai": [ + + ], + "csi": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.csi:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "versions": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:51:53.9057" + }, + "crai - stub": { + "content": [ + { + "0": [ + + ], + "1": [ + + ], + "2": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.recalibrated.sorted.cram.crai:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "3": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ], + "bai": [ + + ], + "crai": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.recalibrated.sorted.cram.crai:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "csi": [ + + ], + "versions": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:51:45.931558" + }, + "bai - stub": { + "content": [ + { + "0": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.bai:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "1": [ + + ], + "2": [ + + ], + "3": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ], + "bai": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.bai:md5,d41d8cd98f00b204e9800998ecf8427e" + ] + ], + "crai": [ + + ], + "csi": [ + + ], + "versions": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:51:34.807525" + }, + "csi": { + "content": [ + "test.paired_end.sorted.bam.csi", + [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:52:55.688799" + }, + "crai": { + "content": [ + { + "0": [ + + ], + "1": [ + + ], + "2": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.recalibrated.sorted.cram.crai:md5,14bc3bd5c89cacc8f4541f9062429029" + ] + ], + "3": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ], + "bai": [ + + ], + "crai": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.recalibrated.sorted.cram.crai:md5,14bc3bd5c89cacc8f4541f9062429029" + ] + ], + "csi": [ + + ], + "versions": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:51:17.609533" + }, + "bai": { + "content": [ + { + "0": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.bai:md5,704c10dd1326482448ca3073fdebc2f4" + ] + ], + "1": [ + + ], + "2": [ + + ], + "3": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ], + "bai": [ + [ + { + "id": "test", + "single_end": false + }, + "test.paired_end.sorted.bam.bai:md5,704c10dd1326482448ca3073fdebc2f4" + ] + ], + "crai": [ + + ], + "csi": [ + + ], + "versions": [ + "versions.yml:md5,802c9776d9c5e95314e888cf18e96d77" + ] + } + ], + "meta": { + "nf-test": "0.9.0", + "nextflow": "24.04.3" + }, + "timestamp": "2024-07-22T16:51:04.16585" + } +} \ No newline at end of file diff --git a/modules/nf-core/samtools/index/tests/tags.yml b/modules/nf-core/samtools/index/tests/tags.yml new file mode 100644 index 0000000..e0f58a7 --- /dev/null +++ b/modules/nf-core/samtools/index/tests/tags.yml @@ -0,0 +1,2 @@ +samtools/index: + - modules/nf-core/samtools/index/** diff --git a/nextflow.config b/nextflow.config index 35f7911..536987e 100644 --- a/nextflow.config +++ b/nextflow.config @@ -16,6 +16,10 @@ params { fasta = null header = null + // Aligner option + short_aligner = "minimap2" // Can choose minimap2 and bwamem2 + chunk_size = 10000 + // Output format options outfmt = "cram" compression = "crumble" diff --git a/nextflow_schema.json b/nextflow_schema.json index 2209d6b..91eb9f6 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -86,6 +86,26 @@ } } }, + "aligner_options": { + "title": "Aligner options", + "type": "object", + "description": "", + "default": "", + "properties": { + "short_aligner": { + "type": "string", + "default": "minimap2", + "description": "Alignment method to use for short-reads (options are: \"minimap2\" or \"bwamem2\"", + "enum": ["minimap2", "bwamem2"] + }, + "chunk_size": { + "type": "integer", + "default": 10000, + "description": "Chunk size for parallel processing of CRAM files, in number of containers", + "minimum": 1 + } + } + }, "execution": { "title": "Execution", "type": "object", @@ -285,6 +305,9 @@ { "$ref": "#/definitions/reference_genome_options" }, + { + "$ref": "#/definitions/aligner_options" + }, { "$ref": "#/definitions/execution" }, diff --git a/seq_cache_populate.pl b/seq_cache_populate.pl new file mode 100644 index 0000000..e69de29 diff --git a/subworkflows/local/align_short.nf b/subworkflows/local/align_short.nf index 6b58e4e..4f49cdc 100644 --- a/subworkflows/local/align_short.nf +++ b/subworkflows/local/align_short.nf @@ -2,21 +2,26 @@ // Align short read (HiC and Illumina) data against the genome // -include { SAMTOOLS_FASTQ } from '../../modules/nf-core/samtools/fastq/main' -include { BWAMEM2_MEM } from '../../modules/nf-core/bwamem2/mem/main' -include { SAMTOOLS_MERGE } from '../../modules/nf-core/samtools/merge/main' -include { SAMTOOLS_SORMADUP } from '../../modules/local/samtools_sormadup' - +include { SAMTOOLS_FASTQ } from '../../modules/nf-core/samtools/fastq/main' +include { SAMTOOLS_MERGE } from '../../modules/nf-core/samtools/merge/main' +include { SAMTOOLS_SORMADUP } from '../../modules/local/samtools_sormadup' +include { SAMTOOLS_INDEX } from '../../modules/nf-core/samtools/index/main' +include { GENERATE_CRAM_CSV } from '../../modules/local/generate_cram_csv' +include { SAMTOOLS_SORMADUP as CONVERT_CRAM } from '../../modules/local/samtools_sormadup' +include { SAMTOOLS_ADDREPLACERG } from '../../modules/local/samtools_addreplacerg' +include { MINIMAP2_MAPREDUCE } from '../../subworkflows/local/minimap2_mapreduce' +include { BWAMEM2_MAPREDUCE } from '../../subworkflows/local/bwamem2_mapreduce' workflow ALIGN_SHORT { take: - fasta // channel: [ val(meta), /path/to/fasta ] + fasta // channel: [ val(meta), /path/to/fasta ] reference_tuple index // channel: [ val(meta), /path/to/bwamem2/ ] - reads // channel: [ val(meta), /path/to/datafile ] + reads // channel: [ val(meta), /path/to/datafile ] hic_reads_path main: ch_versions = Channel.empty() + ch_merged_bam = Channel.empty() // Check file types and branch reads @@ -28,23 +33,60 @@ workflow ALIGN_SHORT { | set { ch_reads } - // Convert from CRAM to FASTQ only if CRAM files were provided as input - SAMTOOLS_FASTQ ( ch_reads.cram, true ) - ch_versions = ch_versions.mix ( SAMTOOLS_FASTQ.out.versions.first() ) - - - SAMTOOLS_FASTQ.out.interleaved - | mix ( ch_reads.fastq ) - | set { ch_reads_fastq } - + // Convert FASTQ to CRAM only if FASTQ were provided as input + CONVERT_CRAM ( ch_reads.fastq, fasta ) + ch_versions = ch_versions.mix ( CONVERT_CRAM.out.versions ) + + SAMTOOLS_ADDREPLACERG ( CONVERT_CRAM.out.bam ) + ch_versions = ch_versions.mix ( SAMTOOLS_ADDREPLACERG.out.versions ) + + SAMTOOLS_ADDREPLACERG.out.cram + | mix ( ch_reads.cram ) + | set { ch_reads_cram } + + // Index the CRAM file + SAMTOOLS_INDEX ( ch_reads_cram ) + ch_versions = ch_versions.mix( SAMTOOLS_INDEX.out.versions ) + + ch_reads_cram + | join ( SAMTOOLS_INDEX.out.crai ) + | set { ch_reads_cram_crai } + + + // + // MODULE: generate a CRAM CSV file containing the required parametres for CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT + // + GENERATE_CRAM_CSV( ch_reads_cram_crai ) + ch_versions = ch_versions.mix( GENERATE_CRAM_CSV.out.versions ) + + // + // SUBWORKFLOW: mapping hic reads using minimap2 or bwamem2 + // + if (params.short_aligner.startsWith("minimap")) { + MINIMAP2_MAPREDUCE ( + fasta, + GENERATE_CRAM_CSV.out.csv + ) + ch_versions = ch_versions.mix( MINIMAP2_MAPREDUCE.out.versions ) + ch_merged_bam = ch_merged_bam.mix(MINIMAP2_MAPREDUCE.out.mergedbam) + } else { + BWAMEM2_MAPREDUCE ( + fasta, + GENERATE_CRAM_CSV.out.csv, + index + ) + ch_versions = ch_versions.mix( BWAMEM2_MAPREDUCE.out.versions ) + ch_merged_bam = ch_merged_bam.mix(BWAMEM2_MAPREDUCE.out.mergedbam) + } - // Align Fastq to Genome and output sorted BAM - BWAMEM2_MEM ( ch_reads_fastq, index, fasta, true ) - ch_versions = ch_versions.mix ( BWAMEM2_MEM.out.versions.first() ) + ch_merged_bam + | combine( ch_reads_cram_crai ) + | map { meta_bam, bam, meta_cram, cram, crai -> [ meta_cram, bam ] } + | set { ch_merged_bam } - // Collect all BWAMEM2 output by sample name - BWAMEM2_MEM.out.bam + // Collect all BAM output by sample name + ch_merged_bam | map { meta, bam -> [['id': meta.id.split('_')[0..-2].join('_'), 'datatype': meta.datatype], meta.read_count, bam] } | groupTuple( by: [0] ) | map { meta, read_counts, bams -> [meta + [read_count: read_counts.sum()], bams] } @@ -58,7 +100,7 @@ workflow ALIGN_SHORT { // Merge, but only if there is more than 1 file SAMTOOLS_MERGE ( ch_bams.multi_bams, [ [], [] ], [ [], [] ] ) - ch_versions = ch_versions.mix ( SAMTOOLS_MERGE.out.versions.first() ) + ch_versions = ch_versions.mix ( SAMTOOLS_MERGE.out.versions ) SAMTOOLS_MERGE.out.bam diff --git a/subworkflows/local/bwamem2_mapreduce.nf b/subworkflows/local/bwamem2_mapreduce.nf new file mode 100644 index 0000000..13711fb --- /dev/null +++ b/subworkflows/local/bwamem2_mapreduce.nf @@ -0,0 +1,85 @@ +#!/usr/bin/env nextflow + +// +// MODULE IMPORT BLOCK +// +include { CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT } from '../../modules/local/cram_filter_align_bwamem2_fixmate_sort' +include { SAMTOOLS_MERGE } from '../../modules/nf-core/samtools/merge/main' + +workflow BWAMEM2_MAPREDUCE { + take: + fasta // Channel: tuple [ val(meta), path( file ) ] + csv_ch + index + + + main: + ch_versions = Channel.empty() + mappedbam_ch = Channel.empty() + + + csv_ch + .splitCsv() + .map{ cram_id, cram_info -> + tuple([ + id: cram_id.id, + chunk_id: cram_id.id + "_" + cram_info[5] + ], + file(cram_info[0]), + cram_info[1], + cram_info[2], + cram_info[3], + cram_info[4], + cram_info[5], + cram_info[6] + ) + } + .set { ch_filtering_input } + + + // + // MODULE: map hic reads in each chunk using bwamem2 + // + CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT ( + ch_filtering_input, + fasta, + index + + ) + ch_versions = ch_versions.mix( CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT.out.versions ) + mappedbam_ch = CRAM_FILTER_ALIGN_BWAMEM2_FIXMATE_SORT.out.mappedbam + + // + // LOGIC: PREPARING BAMS FOR MERGE + // + mappedbam_ch + .map{ meta, file -> + tuple( file ) + } + .collect() + .map { file -> + tuple ( + [ + id: file[0].toString().split('/')[-1].split('_')[0] + '_' + file[0].toString().split('/')[-1].split('_')[1] + ], + file + ) + } + .set { collected_files_for_merge } + + + // + // MODULE: MERGE POSITION SORTED BAM FILES AND MARK DUPLICATES + // + SAMTOOLS_MERGE ( + collected_files_for_merge, + fasta, + [ [], [] ] + ) + ch_versions = ch_versions.mix ( SAMTOOLS_MERGE.out.versions.first() ) + + + emit: + mergedbam = SAMTOOLS_MERGE.out.bam + versions = ch_versions.ifEmpty(null) +} diff --git a/subworkflows/local/minimap2_mapreduce.nf b/subworkflows/local/minimap2_mapreduce.nf new file mode 100644 index 0000000..35b5aae --- /dev/null +++ b/subworkflows/local/minimap2_mapreduce.nf @@ -0,0 +1,101 @@ +#!/usr/bin/env nextflow + +// +// MODULE IMPORT BLOCK +// +include { CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT } from '../../modules/local/cram_filter_minimap2_filter5end_fixmate_sort' +include { SAMTOOLS_MERGE } from '../../modules/nf-core/samtools/merge/main' +include { MINIMAP2_INDEX } from '../../modules/nf-core/minimap2/index/main' + + +workflow MINIMAP2_MAPREDUCE { + take: + fasta // Channel: tuple [ val(meta), path( file ) ] + csv_ch + + + main: + ch_versions = Channel.empty() + mappedbam_ch = Channel.empty() + + + // + // MODULE: generate minimap2 mmi file + // + MINIMAP2_INDEX ( + fasta + ) + ch_versions = ch_versions.mix( MINIMAP2_INDEX.out.versions ) + + + // + // LOGIC: generate input channel for mapping + // + csv_ch + .splitCsv() + .combine ( fasta ) + .combine ( MINIMAP2_INDEX.out.index ) + .map{ cram_id, cram_info, ref_id, ref_dir, mmi_id, mmi_path-> + tuple([ + id: cram_id.id, + chunk_id: cram_id.id + "_" + cram_info[5] + ], + file(cram_info[0]), + cram_info[1], + cram_info[2], + cram_info[3], + cram_info[4], + cram_info[5], + cram_info[6], + mmi_path.toString(), + ref_dir + ) + } + .set { ch_filtering_input } + + + // + // MODULE: map hic reads by 10,000 container per time + // + CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT ( + ch_filtering_input + + ) + ch_versions = ch_versions.mix( CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT.out.versions ) + mappedbam_ch = CRAM_FILTER_MINIMAP2_FILTER5END_FIXMATE_SORT.out.mappedbam + + + // + // LOGIC: PREPARING BAMS FOR MERGE + // + mappedbam_ch + .map{ meta, file -> + tuple( file ) + } + .collect() + .map { file -> + tuple ( + [ + id: file[0].toString().split('/')[-1].split('_')[0] + '_' + file[0].toString().split('/')[-1].split('_')[1] + ], + file + ) + } + .set { collected_files_for_merge } + + + // + // MODULE: MERGE POSITION SORTED BAM FILES AND MARK DUPLICATES + // + SAMTOOLS_MERGE ( + collected_files_for_merge, + fasta, + [ [], [] ] + ) + ch_versions = ch_versions.mix ( SAMTOOLS_MERGE.out.versions.first() ) + + + emit: + mergedbam = SAMTOOLS_MERGE.out.bam + versions = ch_versions.ifEmpty(null) +} diff --git a/subworkflows/local/prepare_genome.nf b/subworkflows/local/prepare_genome.nf index aa951b6..f7a4754 100644 --- a/subworkflows/local/prepare_genome.nf +++ b/subworkflows/local/prepare_genome.nf @@ -30,7 +30,7 @@ workflow PREPARE_GENOME { // Unmask genome fasta UNMASK ( ch_fasta ) - ch_versions = ch_versions.mix ( UNMASK.out.versions.first() ) + ch_versions = ch_versions.mix ( UNMASK.out.versions ) // Generate BWA index if ( checkShortReads( params.input ) ) { @@ -42,14 +42,14 @@ workflow PREPARE_GENOME { if ( params.bwamem2_index.endsWith('.tar.gz') ) { ch_bwamem2_index = UNTAR ( ch_bwamem ).untar - ch_versions = ch_versions.mix ( UNTAR.out.versions.first() ) + ch_versions = ch_versions.mix ( UNTAR.out.versions ) } else { ch_bwamem2_index = ch_bwamem } } else { ch_bwamem2_index = BWAMEM2_INDEX ( UNMASK.out.fasta ).index - ch_versions = ch_versions.mix ( BWAMEM2_INDEX.out.versions.first() ) + ch_versions = ch_versions.mix ( BWAMEM2_INDEX.out.versions ) } } else { ch_bwamem2_index = Channel.empty() diff --git a/workflows/readmapping.nf b/workflows/readmapping.nf index c70a406..a71c77e 100644 --- a/workflows/readmapping.nf +++ b/workflows/readmapping.nf @@ -104,7 +104,7 @@ workflow READMAPPING { // ALIGN_HIC ( PREPARE_GENOME.out.fasta, PREPARE_GENOME.out.bwaidx, ch_reads.hic ) ch_versions = ch_versions.mix ( ALIGN_HIC.out.versions ) - + ALIGN_ILLUMINA ( PREPARE_GENOME.out.fasta, PREPARE_GENOME.out.bwaidx, ch_reads.illumina ) ch_versions = ch_versions.mix ( ALIGN_ILLUMINA.out.versions )