https://nf-co.re logo
Join Slack
Powered by
Hi again! Regarding the new version of Sarek (v. ...
# sarek
m
Hi again! Regarding the new version of Sarek (v. 3.7.1), is this the normal behavior when using bbsplit? I have used an standard script that I have for variant calling and annotation: nextflow run $pipeline \ -profile $profile \ -r $version \ --input $samples \ --outdir $outdir \ --tools $tools \ --bbsplit_index /path/bbsplit \ --bbsplit_fasta_list ref_bbsplit.csv \ --genome $genome \ --pon $pon \ --multiqc_title "report" \ -c $config \ --wes \ --pon_tbi $pon_index \ --intervals "HyperExomeV2_capture_targets.bed" \ and I have realised that after bbsplit is applied, the process is finished without going to the next steps (see image). Is there any argument that I am not including to allow Sarek to continue the normal workflow?
r
def not. are you able to reproduce it with a small test example you could share?
what $genome are you using?
what is in your 
config
?
markdownhelp
m
Hi Rike, I am using for the genome --> GATK.GRCh38 parameter. On the other hand, on config file I got this option for bbsplit command --> withName: 'BBMAP_BBSPLIT' { cpus = 24 memory = 90.GB time = 60.h } as I got some issues in the past for RNAseq pipeline. The rest of the config is quite standard, I could share if you need it
r
and what is in 
tools
m
version="3.7.1" pipeline="nf-core/sarek" profile="singularity" genome="GATK.GRCh38" config="pdx11.config" tools="mutect2,snpeff,bbsplit"
r
thank you. I am a bit baffeled. @Maxime Garcia do you spot anything? It should be using the normal bwa indices. @Manuel Ramos can you check out these tests here: https://github.com/nf-core/sarek/blob/master/tests/bbsplit.nf.test and try to run them on the command line, slowly swapping out reference parameters (like the genome, intervals, etc) ? I think we need a minimal example that we can run to debug this. Also if you in the meantime have the nextflow.log, it would help investigate more
m
sure! I'll check it out later thanks. I was wondering if not using the argument _`--save_bbsplit_reads` could make the workflow unable to retreive the full human-specific FastQ to keep processing ? (as it is splitted by the workflow to speed up)
m
that's weird indeed
r
@Manuel Ramos would you mind testing this? it should not matter, because that parameter just copies the reads to the publishDir but nextflow internally is only using data from the work directory. But you never know
á
Hello! Same here, I put bbsplit, strelka, mutect and snpeff but the process is "successfully finished" right after the bbsplit step is done...And the only thing I got were some stats files...
r
Do youhave a tiny example we can test with ?
jsut tested on dev and it works for me 🫠
á
I'm sorry, I don't think I can share my data here to test.. Is it possible to know which data do you get in return? Our problem seems to be that the process works but just in an informative way, we are not getting the human/mouse fastqs splitted... Maybe we lost the point and the bbsplit step is only informative, but I was expecting a process similar to the one in the RNAseq pipeline
r
I am not sure I am following. Can you describe a bit better what you get in the RNASeq pipeline vs what you are missing here. The implementation is basically a lift from the rnaseq pipeline
m
Hi, what Angel is mentioning is that, in this case, once bbsplit is finished multiQC is run directly. In that MultiQC you only obtain info from FastQC, FastP + BBsplit. In addition, the only output you obtain from this run (in folder preprocessing) is a subfolder called 
bbsplit
, with the following: 
ls
0001.PDX11_tumor-1.stats.txt  0002.PDX11_tumor-1.stats.txt  0003.PDX11_tumor-1.stats.txt  0004.PDX11_tumor-1.stats.txt
Those are the stats from each splitted chunk of the original FASTQ file (that's why i suggested if it could be that the merging was not being applied properly). When you run it on RNAseq pipeline, the process goes on with the mapping, etc.
r
yep so this is definitely not the correct or expected behavior. What should happen that after bbsplit is run bwamem is run on the primary reads for each split and sample and markdup then merges the individual bams. However, I cannot reproduce the behavior you are observing neither with the small tests nor with my example above.
t
Is there any chance that BBSplit removes basically all the reads, so that there is nothing left to map?
á
I’m not sure. Is there a threshold number of reads that stops the process?
m
Can BBmap run and not output anything?
r
can you check the work directory of the BBmap process to see if any 
*primary*
reads were produced?
➕ 1
m
HI all, for my case, bbsplit only removes a small amount of the reads (attached image). If you look into the work directory of bbsplit, there is the following: 
0004.PDX11_tumor-1_1.fastp.fastq.gz
0004.PDX11_tumor-1_hs_1.fastq.gz
0004.PDX11_tumor-1.log
0004.PDX11_tumor-1_mm_2.fastq.gz
index_writable  versions.yml
0004.PDX11_tumor-1_2.fastp.fastq.gz
0004.PDX11_tumor-1_hs_2.fastq.gz
0004.PDX11_tumor-1_mm_1.fastq.gz
0004.PDX11_tumor-1.stats.txt      input_index
r
can you share the .command.sh?
m
Copy code
#!/usr/bin/env bash -C -e -u -o pipefail
# If using a pre-built index, create writable structure: symlink all files except
# summary.txt (which we copy to modify). When we stage in the index files the time
# stamps get disturbed, which bbsplit doesn't like. Fix the time stamps in summaries.
if [ "true" == "true" ]; then
    find input_index/ref -type f | while read -r f; do
        target="index_writable/${f#input_index/}"
        mkdir -p "$(dirname "$target")"
        [[ $(basename "$f") == "summary.txt" ]] && cp "$f" "$target" || ln -s "$(realpath "$f")" "$target"
    done
    find index_writable/ref/genome -name summary.txt | while read -r summary_file; do
        src=$(grep '^source' "$summary_file" | cut -f2- -d$'\t' | sed 's|.*/ref/|index_writable/ref/|')
        mod=$(echo "System.out.println(java.nio.file.Files.getLastModifiedTime(java.nio.file.Paths.get(\"$src\")).toMillis());" | jshell -J-Djdk.lang.Process.launchMechanism=vfork - 2>
/dev/null | grep -oE '^[0-9]{12,14}$')
        sed -e 's|bbsplit_index/ref|index_writable/ref|' -e "s|^last modified.*|last modified\t$mod|" "$summary_file" > ${summary_file}.tmp && mv ${summary_file}.tmp ${summary_file}
    done
fi

# Run BBSplit

bbsplit.sh \
    -Xmx73728M \
    path=index_writable \
    threads=24 \
    in=0004.PDX11_tumor-1_1.fastp.fastq.gz in2=0004.PDX11_tumor-1_2.fastp.fastq.gz \
    basename=0004.PDX11_tumor-1_%_#.fastq.gz \
    refstats=0004.PDX11_tumor-1.stats.txt \
    build=1 ambiguous2=best maxindel=150000 2>| >(tee 0004.PDX11_tumor-1.log >&2)

# Summary files will have an absolute path that will make the index
# impossible to use in other processes - fix paths and rename atomically
if [ -d bbsplit_build/ref/genome ]; then
    find bbsplit_build/ref/genome -name summary.txt | while read -r summary_file; do
        sed "s|^source.*|source\t$(grep '^source' "$summary_file" | cut -f2- -d$'\t' | sed 's|.*/bbsplit_build|bbsplit_index|')|" "$summary_file" > ${summary_file}.tmp && mv ${summary_
file}.tmp ${summary_file}
    done
    mv bbsplit_build bbsplit_index
fi

cat <<-END_VERSIONS > versions.yml
"NFCORE_SAREK:sarek:FASTQ_PREPROCESS_GATK:BBMAP_BBSPLIT":
    bbmap: $(bbversion.sh | grep -v "Duplicate cpuset")
END_VERSIONS
r
thank you, and the .command.log as well?
m
Copy code
java -ea -Xmx73728M -Xms73728M -cp /opt/conda/opt/bbmap-39.18-0/current/ align2.BBSplitter ow=t fastareadlen=500 minhits=1 minratio=0.56 maxindel=20 qtrim=rl untrim=t trimq=6 -Xmx73728
M path=index_writable threads=24 in=0004.PDX11_tumor-1_1.fastp.fastq.gz in2=0004.PDX11_tumor-1_2.fastp.fastq.gz basename=0004.PDX11_tumor-1_%_#.fastq.gz refstats=0004.PDX11_tumor-1.sta
ts.txt build=1 ambiguous2=best maxindel=150000
Executing align2.BBSplitter [ow=t, fastareadlen=500, minhits=1, minratio=0.56, maxindel=20, qtrim=rl, untrim=t, trimq=6, -Xmx73728M, path=index_writable, threads=24, in=0004.PDX11_tumo
r-1_1.fastp.fastq.gz, in2=0004.PDX11_tumor-1_2.fastp.fastq.gz, basename=0004.PDX11_tumor-1_%_#.fastq.gz, refstats=0004.PDX11_tumor-1.stats.txt, build=1, ambiguous2=best, maxindel=15000
0]

Merged reference file index_writable/ref/genome/1/merged_ref_7002394016.fa.gz already exists; skipping merge.
Ref merge time:     	0.026 seconds.
Executing align2.BBMap [ow=t, fastareadlen=500, minhits=1, minratio=0.56, maxindel=20, qtrim=rl, untrim=t, trimq=6, threads=24, in=0004.PDX11_tumor-1_1.fastp.fastq.gz, in2=0004.PDX11_t
umor-1_2.fastp.fastq.gz, refstats=0004.PDX11_tumor-1.stats.txt, build=1, ambiguous2=best, maxindel=150000, ref=index_writable/ref/genome/1/merged_ref_7002394016.fa.gz, out_hs=0004.PDX1
1_tumor-1_hs_#.fastq.gz, out_mm=0004.PDX11_tumor-1_mm_#.fastq.gz]
Version 39.18

Set MINIMUM_ALIGNMENT_SCORE_RATIO to 0.560
Reference set statistics will be written to 0004.PDX11_tumor-1.stats.txt
Retaining first best site only for ambiguous mappings.
NOTE:	Deleting contents of index_writable/ref/genome/1 because reference is specified and overwrite=true
NOTE:	Deleting contents of index_writable/ref/index/1 because reference is specified and overwrite=true
Writing reference.
Executing dna.FastaToChromArrays2 [index_writable/ref/genome/1/merged_ref_7002394016.fa.gz, 1, writeinthread=false, genscaffoldinfo=true, retain, waitforwriting=false, gz=true, maxlen=
536670912, writechroms=true, minscaf=1, midpad=300, startpad=8000, stoppad=8000, nodisk=false]

Set genScaffoldInfo=true
Writing chunk 1
Writing chunk 2
Writing chunk 3
Writing chunk 4
Writing chunk 5
Writing chunk 6
Writing chunk 7
Writing chunk 8
Writing chunk 9
Writing chunk 10
Writing chunk 11
Writing chunk 12
Writing chunk 13
Set genome to 1

Loaded Reference:	0.003 seconds.
Loading index for chunk 1-13, build 1
No index available; generating from reference genome: /path/work/52/db5e32f3da0e03670cd515f8461bfb/index
_writable/ref/index/1/chr1-3_index_k13_c2_b1.block
No index available; generating from reference genome:/path/work/52/db5e32f3da0e03670cd515f8461bfb/index
_writable/ref/index/1/chr8-11_index_k13_c2_b1.block
No index available; generating from reference genome: /path/work/52/db5e32f3da0e03670cd515f8461bfb/index
_writable/ref/index/1/chr4-7_index_k13_c2_b1.block
No index available; generating from reference genome: /path/work/52/db5e32f3da0e03670cd515f8461bfb/index
_writable/ref/index/1/chr12-13_index_k13_c2_b1.block
Indexing threads started for block 8-11
Indexing threads started for block 12-13
Indexing threads started for block 0-3
Indexing threads started for block 4-7
Indexing threads finished for block 12-13
Indexing threads finished for block 0-3
Indexing threads finished for block 8-11
Indexing threads finished for block 4-7
Generated Index:	141.646 seconds.
Finished Writing:	7.525 seconds.
Analyzed Index:   	11.573 seconds.
Reads that map to multiple references will be written to the first reference's stream only.
Creating ref-set statistics table:    	0.184 seconds.
Cleared Memory:    	0.228 seconds.
Processing reads in paired-ended mode.
Started read stream.
Started 24 mapping threads.
Detecting finished threads: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23

   ------------------   Results   ------------------   

Genome:                	1
Key Length:            	13
Max Indel:             	150000
Minimum Score Ratio:  	0.56
Mapping Mode:         	normal
Reads Used:           	4868422	(730254101 bases)

Mapping:          	848.687 seconds.
Reads/sec:       	5736.42
kBases/sec:      	860.45


Pairing data:   	pct pairs	num pairs 	pct bases	   num bases

mated pairs:     	 91.8946% 	  2236909 	 91.8952% 	   671068428
bad pairs:       	  1.7892% 	    43552 	  1.7892% 	    13065492
insert size avg: 	  225.47


Read 1 data:      	pct reads	num reads 	pct bases	   num bases

mapped:          	 94.8616% 	  2309132 	 94.8616% 	   346367574
unambiguous:     	 89.5408% 	  2179612 	 89.5408% 	   326939722
ambiguous:       	  5.3208% 	   129520 	  5.3208% 	    19427852
low-Q discards:  	  0.0000% 	        0 	  0.0000% 	           0

perfect best site:	 53.4527% 	  1301151 	 53.4527% 	   195171430
semiperfect site:	 53.4724% 	  1301631 	 53.4724% 	   195243428
rescued:         	  5.6037% 	   136405

Match Rate:      	      NA 	       NA 	 89.9172% 	   330723284
Error Rate:      	 43.5534% 	  1005729 	 10.0815% 	    37080676
Sub Rate:        	 42.5138% 	   981723 	  2.5864% 	     9512889
Del Rate:        	  7.7681% 	   179380 	  5.8295% 	    21441237
Ins Rate:        	 18.5228% 	   427725 	  1.6657% 	     6126550
N Rate:          	  0.1717% 	     3966 	  0.0013% 	        4851


Read 2 data:      	pct reads	num reads 	pct bases	   num bases

mapped:          	 95.0594% 	  2313947 	 95.0594% 	   347085573
unambiguous:     	 89.7244% 	  2184080 	 89.7244% 	   327605898
ambiguous:       	  5.3351% 	   129867 	  5.3351% 	    19479675
low-Q discards:  	  0.0000% 	        0 	  0.0000% 	           0

perfect best site:	 52.6249% 	  1281001 	 52.6252% 	   192147618
semiperfect site:	 52.6428% 	  1281437 	 52.6431% 	   192213018
rescued:         	  6.2850% 	   152991

Match Rate:      	      NA 	       NA 	 89.7773% 	   330694213
Error Rate:      	 44.6128% 	  1032360 	 10.2220% 	    37652502
Sub Rate:        	 43.6285% 	  1009582 	  2.7416% 	    10098493
Del Rate:        	  7.1506% 	   165467 	  5.7727% 	    21263732
Ins Rate:        	 18.0051% 	   416646 	  1.7077% 	     6290277
N Rate:          	  0.0494% 	     1143 	  0.0007% 	        2590

Total time:     	1055.667 seconds.
r
For some reason there is no primary reference specified in the bbsplit call, but this is what the module is looking for to stage data into the output: https://github.com/nf-core/sarek/blob/20f41d1ce8b7ba296ee22adc71fe2da2ebcae93f/modules/nf-core/bbmap/bbsplit/main.nf#L20 to know which reads to continue processing. Normally the primary reference is specified during the index building https://github.com/nf-core/sarek/blob/20f41d1ce8b7ba296ee22adc71fe2da2ebcae93f/modules/nf-core/bbmap/bbsplit/main.nf#L53 if you run the pipeline ahead of time to specify the index. HOw did you build the bbsplit index
a
Hello
did you manage to fix this? I'm having a similar issue with v3.8.1. With tools=bbsplit,mutect2,snpeff the pipeline only runs fastqc, fastp and multiqc, skipping everything else.
r
hm yes this was fixed. Do you have a minimal reproducible example?
a
I can paste the command and the log
The command is the following: nextflow run nfcore/sarek -r 3.8.1 -profile singularity --input samplefile_test.csv -c nextflow.conf --genome GATK.GRCh38 --igenomes_base /mnt/petasan_general_bioinformatics_RW/refs/igenomes --msisensorpro_scan false --msisensor2_models false --wes --tools bbsplit --outdir results_sarek_bbsplit Nextflow version: 25.10.2
Log:
image.png
3 Views
Open in Slack
PreviousNext
Powered by