Skip to content

First pass at adding bamslice #62

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mattsoup wants to merge 5 commits into
base: master
Choose a base branch
from
Open

First pass at adding bamslice #62

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3303f07
First pass at adding bamslice
mattsoup Jun 29, 2026
b007616
Updating bamslice to include adapter dimers, plus corresponding fastp…
Jul 2, 2026
9a0b448
More useful comment
Jul 2, 2026
3701b9c
Updating snapshot to account for new fastp version
mattsoup Jul 2, 2026
2517bed
Updating snapshot from the cluster rather than my laptop
mattsoup Jul 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 0 additions & 2 deletions conf/test.config
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,6 @@ params {
max_time = '6.h'

outputDir = "test_output"
fastp_path = null

}

trace.overwrite = true
Expand Down
49 changes: 24 additions & 25 deletions main.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ try {
include { createVersionsFile } from './lib/versions.nf'
include { format_ngs_agg_opts } from './modules/aggregate_results'
include { fastp } from './modules/fastp'
include { mergeFastpJson } from './modules/merge_fastp_json'
include { alignReads } from './modules/align_reads'
include { mergeAndMarkDuplicates } from './modules/merge_and_mark_duplicates'
include { methylDackel_mbias } from './modules/methyldackel_mbias'
Expand Down Expand Up @@ -83,32 +84,30 @@ workflow {
}

///////// Trim, align and mark duplicates //////////
fastp( passed_bams )
if (params.single_end) {
fastq_chunks = fastp.out.trimmed_fastq
.flatMap { library, fq_files ->
def fq_list = fq_files instanceof List ? fq_files : [fq_files]

fq_list.findAll { it.baseName.contains('.1.trimmed.fastq') }.collect { fq1 ->
def chunk_name = fq1.baseName.split(".1.trimmed.fastq")[0]
[library, chunk_name, fq1]
}
// Split each uBAM into byte-range chunks so trimming runs in parallel per chunk.
chunk_size = params.bamslice_chunk_size as long
bam_chunks = passed_bams.flatMap { library, bam ->
def file_size = bam.size()
def chunks = []
for (long start = 0; start < file_size; start += chunk_size) {
long end = Math.min(start + chunk_size, file_size)
chunks << tuple(library, start, end)
}
chunks
}
Comment on lines +87 to 97
else {
fastq_chunks = fastp.out.trimmed_fastq
.flatMap { library, fq_files ->
def fq_list = fq_files instanceof List ? fq_files : [fq_files]
def chunk_groups = fq_list.groupBy {
it.baseName.replaceAll(/\.[12]\.trimmed\.fastq$/, '')
}

chunk_groups.collect { chunk_prefix, files ->
def r1 = files.find { it.baseName.contains('.1.trimmed.fastq') }
def r2 = files.find { it.baseName.contains('.2.trimmed.fastq') }

[library, chunk_prefix, [r1, r2]]
}

fastp( passed_bams.combine(bam_chunks, by:0) )
mergeFastpJson( fastp.out.fastp_json.groupTuple() )

fastq_chunks = fastp.out.trimmed_fastq.map { library, chunk_name, fq_files ->
def fq_list = fq_files instanceof List ? fq_files : [fq_files]
if (params.single_end) {
tuple(library, chunk_name, fq_list[0])
}
else {
def r1 = fq_list.find { it.name.contains('.1.trimmed.fastq') }
def r2 = fq_list.find { it.name.contains('.2.trimmed.fastq') }
tuple(library, chunk_name, [r1, r2])
}
}
alignReads( passed_bams.combine(fastq_chunks, by:0), params.reference_list.bwa_index )
Expand Down Expand Up @@ -163,7 +162,7 @@ workflow {
['--bam', mergeAndMarkDuplicates.out.md_bams.map{ tuple(it[0], it[1]) }],
['--bai', mergeAndMarkDuplicates.out.md_bams.map{ tuple(it[0], it[2]) }],
['--metadata_bam_file', bams],
['--fastp', fastp.out.fastp_json],
['--fastp', mergeFastpJson.out.merged_json],
['--aln', picard_metrics.out.for_agg ],
['--gc', gc_bias.out.for_agg ],
['--dup', mergeAndMarkDuplicates.out.log],
Expand Down
26 changes: 13 additions & 13 deletions modules/fastp.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,36 +1,36 @@
process fastp {
tag "${library}"
label 'low_cpu'
conda "bioconda::samtools=1.21 bioconda::fastp=1.0.1"
conda "bioconda::samtools=1.21 bioconda::fastp=1.3.6 bioconda::bamslice=0.2.1"
publishDir "${params.outputDir}/fastp"

input:
tuple val(library), path(bam)
tuple val(library), path(bam), val(start_offset), val(end_offset)

output:
tuple val(library), path("*.trimmed.fastq.gz"), emit: trimmed_fastq
tuple val(library), path("${library}.fastp.json"), emit: fastp_json
tuple val(library), val("${library}_${start_offset}_${end_offset}"), path("${library}_${start_offset}_${end_offset}*.trimmed.fastq.gz"), emit: trimmed_fastq
tuple val(library), val("${library}_${start_offset}_${end_offset}"), path("${library}_${start_offset}_${end_offset}.fastp.json"), emit: fastp_json
tuple val("${task.process}"), val('samtools'), eval('samtools --version | head -n 1 | sed \'s/^samtools //\''), topic: versions
tuple val("${task.process}"), val('fastp'), eval('fastp --version 2>&1 | cut -f 2 -d " "'), topic: versions

tuple val("${task.process}"), val('bamslice'), eval('bamslice --version | cut -f 2 -d " "'), topic: versions

script:
def fastp_path = params.fastp_path ? params.fastp_path : ''
def fastp_args = params.single_end ? "--out1 ${library}.1.trimmed.fastq.gz" : "--interleaved_in --out1 ${library}.1.trimmed.fastq.gz --out2 ${library}.2.trimmed.fastq.gz"
def chunk = "${library}_${start_offset}_${end_offset}"
def fastp_args = params.single_end ? "--out1 ${chunk}.1.trimmed.fastq.gz" : "--interleaved_in --out1 ${chunk}.1.trimmed.fastq.gz --out2 ${chunk}.2.trimmed.fastq.gz"
"""
set +o pipefail
inst_name=\$(samtools view ${bam} | head -n 1 | cut -d ":" -f 1)
set -o pipefail

trim_polyg=\$(echo "\${inst_name}" | awk '{if (\$1~/^A0|^NB|^NS|^VH/) {print "--trim_poly_g"} else {print ""}}')
echo \${trim_polyg} | awk '{ if (length(\$1)>0) { print "2-color instrument: poly-g trim mode on" } }'
samtools fastq -n ${bam} | \\
${fastp_path}fastp --stdin \\

bamslice --input ${bam} --start-offset ${start_offset} --end-offset ${end_offset} | \\
fastp --stdin \\
-l 2 -Q \${trim_polyg} \\
--thread 1 \\
--overrepresentation_analysis \\
-j "${library}.fastp.json" \\
--split_by_lines ${params.fastq_split_lines} \\
-j "${chunk}.fastp.json" \\
${fastp_args}
"""
}
18 changes: 18 additions & 0 deletions modules/merge_fastp_json.nf
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
process mergeFastpJson {
tag "${library}"
label 'low_cpu'
conda 'bioconda::bamslice=0.2.1'
publishDir "${params.outputDir}/fastp"

input:
tuple val(library), val(_slices), path(json_files)

output:
tuple val(library), path("${library}.fastp.json"), emit: merged_json
tuple val("${task.process}"), val('bamslice'), eval('bamslice --version | cut -f 2 -d " "'), topic: versions

script:
"""
fastp-merge ${json_files.join(' ')} -o ${library}.fastp.json
"""
}
2 changes: 1 addition & 1 deletion nextflow.config
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ params {
workflow = 'EM-seq'
outputDir = "em-seq_output"
tmp_dir = '/tmp'
fastq_split_lines = 8000000 // approx 4M reads
bamslice_chunk_size = 212121330 // For 2x100 reads, ~4.8M reads
enable_neb_agg = false
skip_target_bed = false
target_bed_slop = 50
Expand Down
20 changes: 10 additions & 10 deletions tests/main.nf.test
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,14 +11,14 @@ nextflow_pipeline {
params.genome = 'test'
params.email = 'foo@example.com'
params.flowcell = "test_pipeline"
params.fastq_split_lines = 4000
params.bamslice_chunk_size = 250000
params.enable_neb_agg = true

}

then {
def bwameth_file_1 = bam("${launchDir}/test_output/bwameth_align/0001.emseq-test1.aln.bam").getStatistics()
def bwameth_file_2 = bam("${launchDir}/test_output/bwameth_align/0002.emseq-test1.aln.bam").getStatistics()
def bwameth_file_1 = bam("${launchDir}/test_output/bwameth_align/emseq-test1_0_250000.aln.bam").getStatistics()
def bwameth_file_2 = bam("${launchDir}/test_output/bwameth_align/emseq-test1_250000_421605.aln.bam").getStatistics()
def fastqc = path("${launchDir}/test_output/stats/fastqc/emseq-test1.md_fastqc.html").exists()
def fastp = path("${launchDir}/test_output/fastp/emseq-test1.fastp.json").text.tokenize('\n')[1..33]
def markdup = bam("${launchDir}/test_output/markduped_bams/emseq-test1.md.bam").getStatistics()
Expand All @@ -30,7 +30,7 @@ nextflow_pipeline {
def alignment_metrics = path("${launchDir}/test_output/stats/picard_alignment_metrics/emseq-test1.alignment_summary_metrics.txt").text.tokenize('\n')[5..8]
def methyldackel_extract = path("${launchDir}/test_output/methylDackelExtracts/emseq-test1_CpG.methylKit.gz").md5
def mbias = path("${launchDir}/test_output/methylDackelExtracts/mbias/emseq-test1.combined_mbias.tsv").md5
def nonconverted = path("${launchDir}/test_output/bwameth_align/0001.emseq-test1.nonconverted_counts.tsv").text.tokenize('\n')
def nonconverted = path("${launchDir}/test_output/bwameth_align/emseq-test1_0_250000.nonconverted_counts.tsv").text.tokenize('\n')
Comment on lines 20 to +33
def nonconverted_combined = path("${launchDir}/test_output/stats/nonconverted_counts/emseq-test1.nonconverted_counts.for_agg.tsv").text.tokenize('\n')
def combined_summaries = path("${launchDir}/test_output/methylKit_intersections/positional_summaries_combined.tsv").exists()
def positional_summary = path("${launchDir}/test_output/methylKit_intersections/emseq-test1_vs_emseq_test_regions_positional_summary.tsv").text.tokenize('\n')
Expand Down Expand Up @@ -77,15 +77,15 @@ nextflow_pipeline {
params.genome = 'test'
params.email = 'foo@example.com'
params.flowcell = "test_pipeline"
params.fastq_split_lines = 4000
params.bamslice_chunk_size = 250000
params.enable_neb_agg = true
params.single_end = true

}

then {
def bwameth_file_1 = bam("${launchDir}/test_output/bwameth_align/0001.emseq-test1.aln.bam").getStatistics()
def bwameth_file_2 = bam("${launchDir}/test_output/bwameth_align/0002.emseq-test1.aln.bam").getStatistics()
def bwameth_file_1 = bam("${launchDir}/test_output/bwameth_align/emseq-test1_0_250000.aln.bam").getStatistics()
def bwameth_file_2 = bam("${launchDir}/test_output/bwameth_align/emseq-test1_250000_421605.aln.bam").getStatistics()
def fastqc = path("${launchDir}/test_output/stats/fastqc/emseq-test1.md_fastqc.html").exists()
def fastp = path("${launchDir}/test_output/fastp/emseq-test1.fastp.json").text.tokenize('\n')[1..33]
def markdup = bam("${launchDir}/test_output/markduped_bams/emseq-test1.md.bam").getStatistics()
Expand All @@ -96,7 +96,7 @@ nextflow_pipeline {
def alignment_metrics = path("${launchDir}/test_output/stats/picard_alignment_metrics/emseq-test1.alignment_summary_metrics.txt").text.tokenize('\n')[5..8]
def methyldackel_extract = path("${launchDir}/test_output/methylDackelExtracts/emseq-test1_CpG.methylKit.gz").md5
def mbias = path("${launchDir}/test_output/methylDackelExtracts/mbias/emseq-test1.combined_mbias.tsv").md5
def nonconverted = path("${launchDir}/test_output/bwameth_align/0001.emseq-test1.nonconverted_counts.tsv").text.tokenize('\n')
def nonconverted = path("${launchDir}/test_output/bwameth_align/emseq-test1_0_250000.nonconverted_counts.tsv").text.tokenize('\n')
Comment on lines 87 to +99
def nonconverted_combined = path("${launchDir}/test_output/stats/nonconverted_counts/emseq-test1.nonconverted_counts.for_agg.tsv").text.tokenize('\n')
def combined_summaries = path("${launchDir}/test_output/methylKit_intersections/positional_summaries_combined.tsv").exists()
def positional_summary = path("${launchDir}/test_output/methylKit_intersections/emseq-test1_vs_emseq_test_regions_positional_summary.tsv").text.tokenize('\n')
Expand Down
Loading
Loading