Hi I am opening a new discussion because this time I have do nf-core #sarek

Hi, I am opening a new discussion because this ti...

RATTINA Vimel

02/23/2026, 3:46 PM

Hi, I am opening a new discussion because this time I have do not have any error however: -the output bam file is still located in work directory and not in my result directory -no vcf / annotated vcf are generated. I tried to run my sample with GRCh37 and GRCh38, with the correct dbsnp, known_indels and vep_genome. Both analyses have no error but no vcf output also. Here the command I used:

nextflow run nf-core/sarek -profile singularity --input /mnt/HDD2/sarek_inhouse/samplesheet_mapping_NGShem.csv --outdir /mnt/HDD2/sarek_inhouse/2026-02-04_SRP_NGShem/ -r 3.8.1 -with-report -with-trace -with-timeline -with-dag -params-file ipa_params.json

Copy code

[...]
Pulling Singularity image <https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/34/34e733a9ae16a27e80fe00f863ea1479c96416017f24a907996126283e7ecd4d/data> [cache /mnt/HDD2/sarek_inhouse/work/singularity/community-cr-prod.seqera.io-docker-registry-v2-blobs-sha256-34-34e733a9ae16a27e80fe00f863ea1479c96416017f24a907996126283e7ecd4d-data.img]
-[nf-core/sarek] Pipeline completed successfully-
WARN: Singularity cache directory has not been defined -- Remote image will be stored in the path: /mnt/HDD2/sarek_inhouse/work/singularity -- Use the environment variable NXF_SINGULARITY_CACHEDIR to specify a different location
Completed at: 19-Feb-2026 19:21:11
Duration    : 1h 51m 56s
CPU hours   : 4.3
Succeeded   : 9

Here the sample sheet:

cat samplesheet_mapping_NGShem.csv

Copy code

patient,sex,status,sample,lane,fastq_1,fastq_2
C2508328-star,NA,1,C2508328-star,L1,/mnt/SHARE_NGSDATA/IPA_NHE/260204_VH01811_146_AAHJLJCM5/Analysis/1/Data/fastq/C2508328-star_S1_R1_001.fastq.gz,/mnt/SHARE_NGSDATA/IPA_NHE/260204_VH01811_146_AAHJLJCM5/Analysis/1/Data/fastq/C2508328-star_S1_R2_001.fastq.gz

Here the parameters:

cat ipa_params.json

Copy code

{
  "igenomes_ignore": true,
  "genome": "GATK.GRCh38",
  "markdup_java_options": "-Xmx20g",
  "umi_read_structure": "5M2S+T 5M2S+T",
  "save_mapped": true,
  "tools": "freebayes,manta,mpileup,vep",
  "vep_cache": "/mnt/ARC/IPA/CBIOMOL_ARCH/References/vep/",
  "vep_cache_version": "114",
  "vep_genome": "GRCh38",
  "vep_species": "homo_sapiens_refseq",
  "dbsnp": "/mnt/HDD2/sarek_inhouse_hg38/bgzip_input/Homo_sapiens_assembly38.dbsnp138.vcf.gz",
  "intervals": "/mnt/ARC/IPA/CBIOMOL_ARCH/MISEQ_bioinfo/Heme.panel/Annotation/BED_hg38/HemePanel.padded.Nochr.hg38.v4.bed",
  "known_indels": "/mnt/HDD2/sarek_inhouse_hg38/bgzip_input/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz",
  "pon": "/mnt/HDD2/sarek_inhouse_hg38/bgzip_input/MuTect2_PON.NEXTSEQ.SNP.hg38.vcf.gz",
  "max_cpus": 24,
  "max_memory": "128 GB"
}

And the log file attached. Do you know how to have a bam file in the result directory and launch the sarek analysis from fastq input to annotated vcf output? I am also trying with Mutect2, I removed it because it was not able to find the tumour sample. Best, Vimel

👀 2

Maxime Garcia

02/23/2026, 5:15 PM

something weird is definitively happening

Maxime Garcia

02/23/2026, 5:18 PM

As fat as I can see, it uses FASTQTOBAM to create some bam, and then BAM2FASTQ to create some fastq. But it seems to stop there and won't go into the next step and align the resulting FASTQ files. Could you check the output from either of these processes so that we can figure if the reads are correct. It seems to me, we still have some UMI related issues

😬 1

George Seed

02/24/2026, 9:07 AM

I might be encountering this too... the pipeline doesn't disastrously crash, but it doesn't... produce output either - v3.8.1

RATTINA Vimel

02/24/2026, 10:44 AM

Hi, Thanks for the debugging, good to know, I did not go further as I was expected a success and found a bam in work/ directory. My bad! I decided to focus on few regions of interests. As I already have a bam from in-house pipeline, I created a sub-bam thanks to

samtools view -L .bed

and made R1/R2 fastq files with

samtools fastq

. Now the fastq files are 1ko instead of 5Go, the same analysis went from ~2H to 1min which is better for debugging. Here the commands:

export JAVA_CMD=/usr/lib/jvm/java-25-openjdk-amd64/bin/java

export JAVA_HOME=/usr/lib/jvm/java-25-openjdk-amd64

nextflow run nf-core/sarek -profile singularity --input /mnt/HDD2/sarek_inhouse_hg19/sub_input/sub_samplesheet_mapping_NGShem.csv --outdir /mnt/HDD2/sarek_inhouse_hg19/2026-02-04_SRP_NGShem/ -r 3.8.1 -with-report -with-trace -with-timeline -with-dag -params-file /mnt/HDD2/sarek_inhouse_hg19/sub_input/sub_ipa_params.json

Copy code

[...]
WARN: FASTQ file(/mnt/HDD2/sarek_inhouse_hg19/sub_input/subC2508328-star_R1.fastq.gz): Cannot extract flowcell ID from @lib2:3347859
WARN: Singularity cache directory has not been defined -- Remote image will be stored in the path: /mnt/HDD2/sarek_inhouse_hg19/work/singularity -- Use the environment variable NXF_SINGULARITY_CACHEDIR to specify a different location
Completed at: 24-Feb-2026 10:49:46
Duration    : 1m
CPU hours   : (a few seconds)
Succeeded   : 11

The sample sheet:

cat sub_samplesheet_mapping_NGShem.csv

Copy code

patient,sex,status,sample,lane,fastq_1,fastq_2
C2508328-star,NA,1,C2508328-star,L1,/mnt/HDD2/sarek_inhouse_hg19/sub_input/subC2508328-star_R1.fastq.gz,/mnt/HDD2/sarek_inhouse_hg19/sub_input/subC2508328-star_R2.fastq.gz

The params.json, I modified it thanks to

~/.nextflow/assets/nf-core/sarek/nextflow.config

inspiration:

cat sub_input/sub_ipa_params.json

Copy code

{
  "igenomes_ignore": true,
  "genome": "GATK.GRCh37",
  "markdup_java_options": "-Xmx20g",
  "umi_read_structure": "5M2S+T 5M2S+T",
  "save_mapped": true,
  "save_output_as_bam": true,
  "tools": "strelka,mutect2,freebayes,manta,mpileup,vep",
  "skip_tools": "markduplicates",
  "wes": true,
  "vep_cache": "/mnt/ARC/IPA/CBIOMOL_ARCH/References/vep/",
  "vep_cache_version": "114",
  "vep_genome": "GRCh37",
  "vep_species": "homo_sapiens_refseq",
  "dbsnp": "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/dbsnp_138.b37.vcf.gz",
  "intervals": "/mnt/HDD2/sarek_inhouse_hg19/sub_input/sub_HemePanel.padded.Nochr.v4.bed",
  "known_indels": "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/1000G_phase1.indels.b37.vcf.gz",
  "known_snps": "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/1000G_phase1.snps.high_confidence.b37.vcf.gz",
  "pon": "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/MuTect2_PON.vcf.gz",
  "max_cpus": 90,
  "max_memory": "120 GB"
}

Regarding the output, thanks to the execution trace I know which outputs are in which directories and can check them:

cut -f2,4,5 /mnt/HDD2/sarek_inhouse_hg19/2026-02-04_SRP_NGShem/pipeline_info/execution_trace_2026-02-24_10-48-43.txt

Copy code

hash	name	status
70/3c4727	NFCORE_SAREK:PREPARE_INTERVALS:CREATE_INTERVALS_BED (sub_HemePanel.padded.Nochr.v4.bed)	COMPLETED
ff/b630f0	NFCORE_SAREK:PREPARE_GENOME:TABIX_PON (MuTect2_PON.vcf)	COMPLETED
3e/1f88fe	NFCORE_SAREK:PREPARE_INTERVALS:TABIX_BGZIPTABIX_INTERVAL_COMBINED (sub_HemePanel.padded.Nochr.v4)	COMPLETED
da/913111	NFCORE_SAREK:PREPARE_INTERVALS:TABIX_BGZIPTABIX_INTERVAL_SPLIT (2_25470896-25471131)	COMPLETED
1c/701730	NFCORE_SAREK:PREPARE_GENOME:TABIX_KNOWN_INDELS (1000G_phase1.indels.b37.vcf)	COMPLETED
d1/d0b645	NFCORE_SAREK:sarek:FASTQC (C2508328-star-L1)	COMPLETED
05/9f23ab	NFCORE_SAREK:sarek:FASTQ_PREPROCESS_GATK:FASTQ_CREATE_UMI_CONSENSUS_FGBIO:FASTQTOBAM (C2508328-star-L1)	COMPLETED
21/2b3eda	NFCORE_SAREK:sarek:FASTQ_PREPROCESS_GATK:FASTQ_CREATE_UMI_CONSENSUS_FGBIO:BAM2FASTQ (C2508328-star-L1)	COMPLETED
c0/5b0885	NFCORE_SAREK:PREPARE_GENOME:TABIX_KNOWN_SNPS (1000G_phase1.snps.high_confidence.b37.vcf)	COMPLETED
b9/4a20f9	NFCORE_SAREK:PREPARE_GENOME:TABIX_DBSNP (dbsnp_138.b37.vcf)	COMPLETED
b1/d5e080	NFCORE_SAREK:sarek:MULTIQC	COMPLETED

I can find these files and they are not empty (16 alignments):

samtools view -h work/05/*/C2508328-star-L1.bam > C2508328-star-L1.sam

grep -Pvc "^@" C2508328-star-L1.sam

Copy code

zcat work/21/*/C2508328-star-L1_interleaved.fq.gz | grep -Pc "^\+$"

Copy code

However, there are not other files (expect the input files with symbolic links). In our in-house pipeline, after

fgbio FastqToBam

, we are using

samtools fastq

then

bwa mem

. I am not familiar with

fgbio bam2fastq

and the

interleaved.fq.gz

file. Do you see something anormal? I also joined the log file. Best, Vimel

Thomas Danhorn

02/25/2026, 5:56 AM

What you are doing with

samtools view

and

grep

can ne more succinctly done with

samtools view -c

(you can also add other option like

-f

-F

-q

, etc., for filtering, e.g. to only count R1, or only uniquely mapping reads). For counting reads in FASTQs, I use

wc -l

and divide by

. If I am interpreting your finding correctly, you only have 16 reads in your BAM and FASTQ -- even for a test that would be way too few -- unless I am missing something, something seems to be very wrong here, and I would fix that first, before trying to debug why your output is missing (not sure if sarek requires a minimum number of aligned reads like #CE8SSJV3N, but it may well be related). Once you have figured that out and still have issues, you could try to see if there is a regression by trying with a few older versions of sarek (you can use

-r

to pick the version, e.g.

-r 3.3.2

RATTINA Vimel

02/25/2026, 8:00 AM

Hi Thomas, Thanks for the

samtools view -c

tip, I also like to look at the sam file this is why I converted it first. For counting reads, there are many ways to do, if we found the correct number, all methods are ok. I have the same error with my entire fastq and bed with GRCh38 so it is not linked to a minimum number of aligned reads.

Copy code

samtools view -c work/7b/*/C2508328-star-L1.bam 
146208632

zcat work/10/*/C2508328-star-L1_interleaved.fq.gz  | grep -Pc "^\+$"
146208632

I reduce the input data at maximum to quick up the analysis and go up to the end. However the subset data are coming from a deduplicated bam file so you may have right on this point. I will give a try with an old version. Best, Vimel

🤞 1

👍 1

RATTINA Vimel

03/05/2026, 8:43 AM

Hi everybody, I tried with an old version

-r 3.3.2

but I still do not have a proper bam in output directory. It seems there is something wrong with the UMI step with UMI inside the sequences at least. What I tried this time was to first put the UMI on the header with umi_tools (take ~3H):

Copy code

umi_tools extract --extract-method=regex --bc-pattern='^(?P<umi_1>.{5})(?P<discard_1>.{2}).*' --bc-pattern2='^(?P<umi_1>.{5})(?P<discard_1>.{2}).*' -I before_umitools/C2508328-star_S1_R1_001.fastq.gz --read2-in=before_umitools/C2508328-star_S1_R2_001.fastq.gz --stdout=header/C2508328-star_R1.fastq.gz --read2-out=header/C2508328-star_R2.fastq.gz

-I

and

--reads2

are the input files and

--stdout

and

--read2-out

are the output files. Followed by a syntax editing for both ends:

Copy code

zcat header/C2508328-star_R1.fastq.gz | sed -E 's/ .*//' | sed -E 's/_/:5,5:/' | sed -E 's/(^@.*)([ATGCN]{5})$/\1,\2/' | gzip > header/renamed_C2508328-star_R1.fastq.gz&
mv header/renamed_C2508328-star_R1.fastq.gz header/C2508328-star_R1.fastq.gz

According to https://nf-co.re/fastquorum/1.2.0/docs/usage/, we need to use

umi_in_read_header true

when UMI are in the header and we should follow the fastp syntax: https://github.com/OpenGene/fastp?tab=readme-ov-file#unique-molecular-identifier-umi-processing On fastp page example, the UMI is in R1 only. In my case I have UMI on R1 and R2, I did this kind of reheader with UMI for umiVar variant caller (https://github.com/imgag/umiVar) and I did the same for sarek. Here what it looked like before

umi_tools

(*5M*~2S~_+T_ with the 5M in bold and 2S crossed out):

zcat before_umitools/C2508328-star_S1_R1_001.fastq.gz | head

@VH01811146AAHJLJCM511101184591000 1N0:CGTTCTCTTG+CCACTTGGTG *GTGAG*~AT~_ATGGTGAGCTCTTGGAGGGCGGGGCACAGAGGCCTGCCCCCCTTCATGTATCTCATACC[...]_

zcat before_umitools/C2508328-star_S1_R2_001.fastq.gz | head

@VH01811146AAHJLJCM511101184591000 2N0:CGTTCTCTTG+CCACTTGGTG *CTTGG*~AT~_CTCAAATCTGTCTCTAGCCTGCATGCGAATCATCGACCAGAAACTAATCTGTCCAACTG[...]_ Here after umi_tools and syntax editing:

zcat header/C2508328-star_R1.fastq.gz | head

@VH01811146AAHJLJCM511101184591000:*5,5*:*GTGAG*,CTTGG ATGGTGAGCTCTTGGAGGGCGGGGCACAGAGGCCTGCCCCCCTTCATGTATCTCATACC[...]

zcat header/C2508328-star_R2.fastq.gz | head

@VH01811146AAHJLJCM511101184591000:*5,5*:*GTGAG*,CTTGG CTCAAATCTGTCTCTAGCCTGCATGCGAATCATCGACCAGAAACTAATCTGTCCAACTG[...] Please let me know if it is not correct and this is why sarek cannot go further on analysis. Here the sarek command used:

Copy code

nextflow run nf-core/sarek -profile singularity \
--input /mnt/HDD2/sarek_inhouse_hg19/header/samplesheet_mapping_NGShem_header.csv \
--outdir /mnt/HDD2/sarek_inhouse_hg19/2026-02-04_SRP_NGShem/ \
-r 3.8.1 -with-report -with-trace -with-timeline -with-dag \
-params-file /mnt/HDD2/sarek_inhouse_hg19/header/ipa_params_header.json

The sample sheet:

Copy code

cat /mnt/HDD2/sarek_inhouse_hg19/header/samplesheet_mapping_NGShem_header.csv
patient,sex,status,sample,lane,fastq_1,fastq_2
C2508328-star,NA,1,C2508328-star,L1,/mnt/HDD2/sarek_inhouse_hg19/header/C2508328-star_R1.fastq.gz,/mnt/HDD2/sarek_inhouse_hg19/header/C2508328-star_R2.fastq.gz

And the params file:

Copy code

cat /mnt/HDD2/sarek_inhouse_hg19/header/ipa_params_header.json
{
    "igenomes_ignore":           true,
    "genome":                    "GATK.GRCh37",
    "markdup_java_options":      "-Xmx20g",
    "umi_in_read_header":        true,
    "save_mapped":               true,
    "save_output_as_bam":        true,
    "groupreadsbyumi_strategy":  "Adjacency",
    "groupreadsbyumi_edits":     1,
    "call_min_reads":            1,
    "call_min_baseq":            20,
    "filter_min_reads":          1,
    "filter_min_baseq":          20,
    "tools":                     "strelka,mutect2,freebayes,manta,mpileup,vep",
    "skip_tools":                "markduplicates",
    "wes":                       true,
    "vep_cache":                 "/mnt/ARC/IPA/CBIOMOL_ARCH/References/vep/",
    "vep_cache_version":         "114",
    "vep_genome":                "GRCh37",
    "vep_species":               "homo_sapiens_refseq",
    "dbsnp":                     "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/dbsnp_138.b37.vcf.gz",
    "intervals":                 "/mnt/ARC/IPA/CBIOMOL_ARCH/MISEQ_bioinfo/Heme.panel/Annotation/HemePanel.padded.Nochr.v4.bed",
    "known_indels":              "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/1000G_phase1.indels.b37.vcf.gz",
    "known_snps":                "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/1000G_phase1.snps.high_confidence.b37.vcf.gz",
    "pon":                       "/mnt/HDD2/sarek_inhouse_hg19/bgzip_input/MuTect2_PON.vcf.gz",
    "max_cpus":                  90,
    "max_memory":                "120 GB"
}

I also added some options specific to fgbio to precise I expect an UMI treatment (

groupreadsbyumi_

_min_reads

and

_min_baseq

). To let you know I also saw we can put

read_structure

in the sample sheet. I tried all combination with and without read structure information in samplesheet and/or params file. No success for all of them. Here the execution trace:

Copy code

cut -f2,4,5 2026-02-04_SRP_NGShem/pipeline_info/execution_trace_2026-03-04_12-44-08.txt
hash    name    status
2f/ea23de       NFCORE_SAREK:PREPARE_INTERVALS:CREATE_INTERVALS_BED (HemePanel.padded.Nochr.v4.bed)     COMPLETED
ff/aebe92       NFCORE_SAREK:PREPARE_GENOME:TABIX_PON (MuTect2_PON.vcf) COMPLETED
d5/3badea       NFCORE_SAREK:PREPARE_INTERVALS:TABIX_BGZIPTABIX_INTERVAL_COMBINED (HemePanel.padded.Nochr.v4)   COMPLETED
4f/ce894f       NFCORE_SAREK:PREPARE_GENOME:TABIX_KNOWN_INDELS (1000G_phase1.indels.b37.vcf)    COMPLETED
2b/7da49d       NFCORE_SAREK:PREPARE_INTERVALS:TABIX_BGZIPTABIX_INTERVAL_SPLIT (1_1718760-1718886)      COMPLETED
65/1406a7       NFCORE_SAREK:PREPARE_GENOME:TABIX_KNOWN_SNPS (1000G_phase1.snps.high_confidence.b37.vcf)        COMPLETED
86/8fe209       NFCORE_SAREK:PREPARE_GENOME:TABIX_DBSNP (dbsnp_138.b37.vcf)     COMPLETED
89/7e0f4e       NFCORE_SAREK:sarek:FASTQ_PREPROCESS_GATK:FASTP (C2508328-star-L1)       COMPLETED
60/9d90d1       NFCORE_SAREK:sarek:FASTQC (C2508328-star-L1)    COMPLETED
c1/b6476b       NFCORE_SAREK:sarek:MULTIQC      COMPLETED

Do you provide an sarek dataset with UMI (inside header, inside sequence) to test the pipeline? It could help to see what intermerdiar files are created to have a better understanding of what is wrong with my in-house data. Best, Vimel

2 Views

Open in Slack

Previous Next