Skip to content

Commit 5ae22e9

Browse files
guillaumecharbonnierguillaumecharbonnier
committed
Misc update prior to fishpond script for Salva
2 parents 6953ead + ed2e7a4 commit 5ae22e9

14 files changed

Lines changed: 424 additions & 32 deletions

‎src/snakemake/envs/arriba.yaml‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,4 +5,6 @@ channels:
55
- bioconda
66
- nodefaults
77
dependencies:
8-
- arriba
8+
# Pinning version to get consistent annotations between hosts for TALL
9+
# but more recent version should work too for independent projects
10+
- arriba = 2.4.0

‎src/snakemake/envs/cytonaut.yaml‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
# Cytonaut requirements
2+
# Will likely be splitted into star and cite-seq-count envs later
3+
channels:
4+
- conda-forge
5+
- bioconda
6+
- nodefaults
7+
dependencies:
8+
- star == 2.7.10b
9+
- cite-seq-count

‎src/snakemake/envs/jbrowser.yaml‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
name:
2+
jbrowser
3+
channels:
4+
- conda-forge
5+
- bioconda
6+
- nodefaults
7+
dependencies:
8+
# From https://jbrowse.org/jb2/docs/quickstart_web/#prerequisites
9+
- nodejs
10+
- samtools
11+
- tabix

‎src/snakemake/envs/mixcr.yaml‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
channels:
2+
- milaboratories
3+
- bioconda
4+
- conda-forge
5+
- nodefaults
6+
dependencies:
7+
- mixcr

‎src/snakemake/envs/star.yaml‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,5 +3,7 @@ channels:
33
- bioconda
44
- nodefaults
55
dependencies:
6-
- star < 2.7.9 # 2.7.10 seems to lead to 'kernel too old' (untested)
6+
# - star < 2.7.9 # 2.7.10 seems to lead to 'kernel too old' (untested)
7+
# On 2025-02-19, I am trying to get a more recent version of STAR to match Cytonaut requirement
8+
- star
79
- samtools

‎src/snakemake/rules/agent_trimmer.smk‎

Lines changed: 25 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -22,8 +22,9 @@ rule agent_trim_pe:
2222
output:
2323
fq1 = "out/{tool}{extra}/{filler}_1.fastq.gz",
2424
fq2 = "out/{tool}{extra}/{filler}_2.fastq.gz",
25-
mbc = "out/{tool}{extra}/{filler}.txt.gz",
26-
properties = "out/{tool}{extra}/{filler}.properties",
25+
mbc = "out/{tool}{extra}/{filler}.txt.gz"#,
26+
# properties seems to be not generated anymore in latest versions of AGeNT
27+
#properties = "out/{tool}{extra}/{filler}.properties",
2728
params:
2829
outprefix = "out/{tool}{extra}/{filler}",
2930
extra = params_extra
@@ -35,13 +36,32 @@ rule agent_trim_pe:
3536
"../envs/agent.yaml"
3637
shell:
3738
"""
38-
{input.agent} trim -fq1 {input.fq1} -fq2 {input.fq2} -out {params.outprefix} {params.extra}
39-
39+
bash {input.agent} trim -fq1 {input.fq1} -fq2 {input.fq2} -out {params.outprefix} {params.extra}
40+
# Try to add a sleep here to avoid some issue with rclone VFS
41+
#timeout=120
42+
#interval=2
43+
#elapsed=0
44+
#
45+
#while [ $elapsed -lt $timeout ]; do
46+
# if [ -f "{params.outprefix}_R1.fastq.gz" ] && [ -f "{params.outprefix}_R2.fastq.gz" ] && \
47+
# [ -f "{params.outprefix}_MBC.txt.gz" ] && [ -f "{params.outprefix}_STATS.properties" ]; then
48+
# break
49+
# fi
50+
# sleep $interval
51+
# elapsed=$((elapsed + interval))
52+
#done
53+
#
54+
#if [ $elapsed -ge $timeout ]; then
55+
# echo "Timeout: Output files did not appear within $timeout seconds." >&2
56+
# exit 1
57+
#fi
58+
59+
# Then move files if they are available:
4060
mv {params.outprefix}_R1.fastq.gz {output.fq1}
4161
mv {params.outprefix}_R2.fastq.gz {output.fq2}
4262
mv {params.outprefix}_MBC.txt.gz {output.mbc}
43-
mv {params.outprefix}_STATS.properties {output.properties}
4463
"""
64+
# mv {params.outprefix}_STATS.properties {output.properties}
4565

4666
rule agent_trim_2_lanes_pe:
4767
"""

‎src/snakemake/rules/agent_wget.smk‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,12 +6,27 @@ rule agent_wget:
66
"""
77
output:
88
agent = "out/agent/agent/agent.sh",
9-
locatit = "out/agent/agent/lib/locatit-2.0.5.jar"
9+
# Note locatit was removed in later versions of AGeNT. Currently commented but it would require adjustment to downstream RNA-Seq analyses.
10+
# locatit = "out/agent/agent/lib/locatit-2.0.5.jar"
1011
shell:
1112
"""
1213
OUTDIR=out/agent
1314
mkdir -p $OUTDIR
1415
cd $OUTDIR
16+
# TODO: version is hardcoded in name but not the target url.
17+
# Need to find a specific URL to get the needed version.
18+
# This is actually not the 2.0.5 version.
19+
# The downstream rules expect these libs:
20+
# -rw-rw-r-- 1 gcharbonnier shinyusers 20K Aug 26 2022 Agent-3.0.6.jar
21+
# -rw-rw-r-- 1 gcharbonnier shinyusers 4,5M Mar 9 2022 creak-1.0.5.jar
22+
# -rw-rw-r-- 1 gcharbonnier shinyusers 119K Nov 18 2021 junit-3.8.1.jar
23+
# -rw-rw-r-- 1 gcharbonnier shinyusers 6,3M Sep 25 2020 locatit-2.0.5.jar
24+
# -rw-rw-r-- 1 gcharbonnier shinyusers 22M Aug 26 2022 trimmer-3.0.5.jar
25+
# But the current download leads to these libs:
26+
# -rw-rw-r-- 1 gcharbonnier shinyusers 16K Jun 28 2024 Agent-3.1.2.jar
27+
# -rw-rw-r-- 1 gcharbonnier shinyusers 4,3M Apr 24 2024 creak-1.1.1.jar
28+
# -rw-rw-r-- 1 gcharbonnier shinyusers 119K Nov 18 2021 junit-3.8.1.jar
29+
# -rw-rw-r-- 1 gcharbonnier shinyusers 22M Jun 28 2024 trimmer-3.1.2.jar
1530
wget 'https://dt4ei3l3hxs7z.cloudfront.net/?elqTrackId=30b3c5b8c3bd44f7b3a01b66ab2a30a5&elqaid=3928&elqat=2' --output-document 'AGeNT_2.0.5.zip'
1631
unzip AGeNT_2.0.5.zip
1732
chmod +x {output.agent}

‎src/snakemake/rules/arriba.smk‎

Lines changed: 125 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,8 @@ rule arriba_extra:
22
"""
33
Created:
44
2024-03-29 20:01:56
5-
Note:
6-
The run_arriba.sh script is available here:
7-
https://github.com/suhrig/arriba/blob/master/run_arriba.sh
8-
and I could improve my snakemake workflow splitting it into the STAR rule, the arriba call and the samtools index.
5+
Modified:
6+
2025-10-28 - Replaced run_arriba.sh with direct STAR/arriba calls to support --outTmpDir
97
Test:
108
out/arriba/pe_fq_GRCh38viral_ENSEMBL104/agent/trim_-v2/ln/alias/sst/all_samples/fastq/113281_PICCL/fusions.tsv
119
"""
@@ -19,7 +17,8 @@ rule arriba_extra:
1917
bam = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}/Aligned.sortedByCoord.out.bam",
2018
bai = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}/Aligned.sortedByCoord.out.bam.bai"
2119
params:
22-
outdir = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}"
20+
outdir = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}",
21+
tmpdir = "/tmp/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}"
2322
wildcard_constraints:
2423
tool="arriba/pe_fq",
2524
arriba_fa_id = "GRCh38viral",
@@ -30,34 +29,139 @@ rule arriba_extra:
3029
8
3130
shell:
3231
"""
32+
WDIR=`pwd`
3333
ARRIBA_FILES=$CONDA_PREFIX/var/lib/arriba
3434
35-
# Using arriba_fa_id we can decide which suffix to use for blacklist, known_fusions and protein_domains:
35+
# Clean tmp dir if previous execution was interrupted
36+
rm -rf {params.tmpdir}
37+
mkdir -p `dirname {params.tmpdir}`
3638
39+
# Using arriba_fa_id we can decide which suffix to use for blacklist, known_fusions and protein_domains:
3740
if [[ {wildcards.arriba_fa_id} == *"GRCh38"* ]] || [[ {wildcards.arriba_fa_id} == *"hg38"* ]]; then
38-
STEM="_hg38_GRCh38_v2.4.0"
41+
BASE="_hg38_GRCh38"
3942
elif [[ {wildcards.arriba_fa_id} == *"GRCh37"* ]] || [[ {wildcards.arriba_fa_id} == *"hg19"* ]]; then
40-
STEM="_hg19_hs37d5_GRCh37_v2.4.0"
43+
BASE="_hg19_hs37d5_GRCh37"
4144
else
4245
echo "Invalid arriba_fa_id"
4346
exit 1
4447
fi
48+
49+
# Auto-detect the version from available files in ARRIBA_FILES
50+
BLACKLIST_FILE=$(ls $ARRIBA_FILES/blacklist${{BASE}}_v*.tsv.gz 2>/dev/null | head -n1)
51+
if [[ -z "$BLACKLIST_FILE" ]]; then
52+
echo "Error: No blacklist file found matching pattern: $ARRIBA_FILES/blacklist${{BASE}}_v*.tsv.gz"
53+
exit 1
54+
fi
55+
56+
# Extract the full stem (including version) from the detected file
57+
STEM=$(basename "$BLACKLIST_FILE" .tsv.gz | sed "s/blacklist//")
4558
4659
cd {params.outdir}
4760
48-
run_arriba.sh \
49-
{WDIR}/out/arriba/download_references_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id}/STAR_index_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id} \
50-
{WDIR}/{input.gtf} \
51-
{WDIR}/{input.fa} \
52-
$ARRIBA_FILES/blacklist$STEM.tsv.gz \
53-
$ARRIBA_FILES/known_fusions$STEM.tsv.gz \
54-
$ARRIBA_FILES/protein_domains$STEM.gff3 \
55-
{threads} \
56-
{WDIR}/{input.fq1} \
57-
{WDIR}/{input.fq2}
58-
61+
# Run STAR with arriba settings and output to tmpdir for FIFO operations
62+
STAR \
63+
--runThreadN {threads} \
64+
--genomeDir $WDIR/out/arriba/download_references_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id}/STAR_index_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id} \
65+
--genomeLoad NoSharedMemory \
66+
--readFilesIn $WDIR/{input.fq1} $WDIR/{input.fq2} \
67+
--readFilesCommand zcat \
68+
--outStd BAM_Unsorted \
69+
--outSAMtype BAM Unsorted \
70+
--outSAMunmapped Within \
71+
--outBAMcompression 0 \
72+
--outTmpDir {params.tmpdir} \
73+
--outFilterMultimapNmax 50 \
74+
--peOverlapNbasesMin 10 \
75+
--alignSplicedMateMapLminOverLmate 0.5 \
76+
--alignSJstitchMismatchNmax 5 -1 5 5 \
77+
--chimSegmentMin 10 \
78+
--chimOutType WithinBAM HardClip \
79+
--chimJunctionOverhangMin 10 \
80+
--chimScoreDropMax 30 \
81+
--chimScoreJunctionNonGTAG 0 \
82+
--chimScoreSeparation 1 \
83+
--chimSegmentReadGapMax 3 \
84+
--chimMultimapNmax 50 | \
85+
tee Aligned.out.bam | \
86+
arriba \
87+
-x /dev/stdin \
88+
-o fusions.tsv \
89+
-O fusions.discarded.tsv \
90+
-a $WDIR/{input.fa} \
91+
-g $WDIR/{input.gtf} \
92+
-b $ARRIBA_FILES/blacklist$STEM.tsv.gz \
93+
-k $ARRIBA_FILES/known_fusions$STEM.tsv.gz \
94+
-t $ARRIBA_FILES/known_fusions$STEM.tsv.gz \
95+
-p $ARRIBA_FILES/protein_domains$STEM.gff3
96+
97+
# Sort and index BAM
98+
samtools sort -@ {threads} -m $((40000/{threads}))M -T tmp -O bam Aligned.out.bam > Aligned.sortedByCoord.out.bam
99+
rm -f Aligned.out.bam
100+
samtools index Aligned.sortedByCoord.out.bam
59101
"""
60102

103+
# Commenting original rule for reference.
104+
# I had to replicate run_arriba.sh content here to be able to use --outTmpDir in STAR call to avoid FIFO issue
105+
# rule arriba_extra:
106+
# """
107+
# Created:
108+
# 2024-03-29 20:01:56
109+
# Note:
110+
# The run_arriba.sh script is available here:
111+
# https://github.com/suhrig/arriba/blob/master/run_arriba.sh
112+
# and I could improve my snakemake workflow splitting it into the STAR rule, the arriba call and the samtools index.
113+
# Test:
114+
# out/arriba/pe_fq_GRCh38viral_ENSEMBL104/agent/trim_-v2/ln/alias/sst/all_samples/fastq/113281_PICCL/fusions.tsv
115+
# """
116+
# input:
117+
# fq1 = "out/{filler}_1.fastq.gz",
118+
# fq2 = "out/{filler}_2.fastq.gz",
119+
# gtf = "out/arriba/download_references_{arriba_fa_id}_{arriba_gtf_id}/{arriba_gtf_id}.gtf",
120+
# fa = "out/arriba/download_references_{arriba_fa_id}_{arriba_gtf_id}/{arriba_fa_id}.fa"
121+
# output:
122+
# tsv = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}/fusions.tsv",
123+
# bam = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}/Aligned.sortedByCoord.out.bam",
124+
# bai = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}/Aligned.sortedByCoord.out.bam.bai"
125+
# params:
126+
# outdir = "out/{tool}_{arriba_fa_id}_{arriba_gtf_id}/{filler}"
127+
# wildcard_constraints:
128+
# tool="arriba/pe_fq",
129+
# arriba_fa_id = "GRCh38viral",
130+
# arriba_gtf_id = "ENSEMBL104"
131+
# conda:
132+
# "../envs/arriba.yaml"
133+
# threads:
134+
# 8
135+
# shell:
136+
# """
137+
# ARRIBA_FILES=$CONDA_PREFIX/var/lib/arriba
138+
139+
# # Using arriba_fa_id we can decide which suffix to use for blacklist, known_fusions and protein_domains:
140+
141+
# if [[ {wildcards.arriba_fa_id} == *"GRCh38"* ]] || [[ {wildcards.arriba_fa_id} == *"hg38"* ]]; then
142+
# STEM="_hg38_GRCh38_v2.5.0"
143+
# elif [[ {wildcards.arriba_fa_id} == *"GRCh37"* ]] || [[ {wildcards.arriba_fa_id} == *"hg19"* ]]; then
144+
# STEM="_hg19_hs37d5_GRCh37_v2.5.0"
145+
# else
146+
# echo "Invalid arriba_fa_id"
147+
# exit 1
148+
# fi
149+
150+
# cd {params.outdir}
151+
152+
# run_arriba.sh \
153+
# {WDIR}/out/arriba/download_references_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id}/STAR_index_{wildcards.arriba_fa_id}_{wildcards.arriba_gtf_id} \
154+
# {WDIR}/{input.gtf} \
155+
# {WDIR}/{input.fa} \
156+
# $ARRIBA_FILES/blacklist$STEM.tsv.gz \
157+
# $ARRIBA_FILES/known_fusions$STEM.tsv.gz \
158+
# $ARRIBA_FILES/protein_domains$STEM.gff3 \
159+
# {threads} \
160+
# {WDIR}/{input.fq1} \
161+
# {WDIR}/{input.fq2}
162+
163+
# """
164+
61165
rule arriba_draw_fusions:
62166
"""
63167
Created:
@@ -91,9 +195,9 @@ rule arriba_draw_fusions:
91195
# Using arriba_fa_id we can decide which suffix to use for cytobands and protein_domains:
92196
93197
if [[ {wildcards.arriba_fa_id} == *"GRCh38"* ]] || [[ {wildcards.arriba_fa_id} == *"hg38"* ]]; then
94-
STEM="_hg38_GRCh38_v2.4.0"
198+
STEM="_hg38_GRCh38_v2.5.0"
95199
elif [[ {wildcards.arriba_fa_id} == *"GRCh37"* ]] || [[ {wildcards.arriba_fa_id} == *"hg19"* ]]; then
96-
STEM="_hg19_hs37d5_GRCh37_v2.4.0"
200+
STEM="_hg19_hs37d5_GRCh37_v2.5.0"
97201
else
98202
echo "Invalid arriba_fa_id"
99203
exit 1

‎src/snakemake/rules/mixcr.smk‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
rule mixcr_analyze_pe:
2+
"""
3+
Created:
4+
2025-05-26
5+
Doc:
6+
https://mixcr.readthedocs.io/en/master/
7+
Aim:
8+
Aligns sequencing reads to V(D)J reference using MiXCR.
9+
Test:
10+
out/mixcr/analyze_rna-seq_--species_hsa/agent/trim_-v2/ln/alias/sst/all_samples/fastq/1379_BOUNAD_RNA_XT_HS2.clna
11+
"""
12+
input:
13+
fq1 = "out/{filler}_1.fastq.gz",
14+
fq2 = "out/{filler}_2.fastq.gz"
15+
# ref = lambda wildcards: eval(mwconf['ids'][wildcards.mixcr_ref_id])
16+
output:
17+
# clna = "out/{tool}{extra}_{mixcr_ref_id}/{filler}.clna"
18+
clna = "out/{tool}{extra}/{filler}.clna"
19+
# log:
20+
# "out/{tool}{extra}/{filler}.log"
21+
# "out/{tool}{extra}_{mixcr_ref_id}/{filler}.log"
22+
benchmark:
23+
# "out/{tool}{extra}_{mixcr_ref_id}/{filler}.benchmark.tsv"
24+
"out/{tool}{extra}/{filler}.benchmark.tsv"
25+
params:
26+
extra = params_extra,
27+
outprefix = "out/{tool}{extra}/{filler}"
28+
wildcard_constraints:
29+
tool = "mixcr/analyze_pe"
30+
conda:
31+
"../envs/mixcr.yaml"
32+
priority: 10
33+
threads: 2
34+
shell:
35+
"""
36+
# mixcr analyze {params.extra} -f {input.fq1} {input.fq2} {params.outprefix}
37+
# Replacing the above command with the following to try to avoid these kinds of issues:
38+
# https://github.com/guillaumecharbonnier/mw-lib/issues/38
39+
tmpdir=$(mktemp -d /tmp/mixcr_analyze_pe.XXXXXX)
40+
tmpprefix=$(basename {params.outprefix})
41+
mkdir -p "$tmpdir"/"$tmpprefix"
42+
mixcr analyze {params.extra} -f {input.fq1} {input.fq2} "$tmpdir"/"$tmpprefix"
43+
mv "$tmpdir"/* $(dirname {params.outprefix})/
44+
rmdir "$tmpdir"
45+
"""
46+
47+
rule mixcr_analyze_se:
48+
"""
49+
Created:
50+
2025-05-26
51+
Doc:
52+
https://mixcr.readthedocs.io/en/master/
53+
Aim:
54+
Aligns sequencing reads to V(D)J reference using MiXCR.
55+
Test:
56+
Never tried
57+
"""
58+
input:
59+
fq = "out/{filler}.fastq.gz",
60+
output:
61+
clna = "out/{tool}{extra}/{filler}.clna"
62+
# log:
63+
# "out/{tool}{extra}/{filler}.log"
64+
# "out/{tool}{extra}_{mixcr_ref_id}/{filler}.log"
65+
benchmark:
66+
# "out/{tool}{extra}_{mixcr_ref_id}/{filler}.benchmark.tsv"
67+
"out/{tool}{extra}/{filler}.benchmark.tsv"
68+
params:
69+
extra = params_extra,
70+
outdir = "out/{tool}{extra}/{filler}"
71+
wildcard_constraints:
72+
tool = "mixcr/analyze_se"
73+
conda:
74+
"../envs/mixcr.yaml"
75+
priority: 10
76+
threads: 2
77+
shell:
78+
"""
79+
mixcr analyze {params.extra} -f {input.fq} {params.outdir}
80+
"""

0 commit comments

Comments
 (0)