Skip to content

Commit f66b26c

Browse files
DavideBrexDavideBrex
committed
feat: fixed bug in ChIP-rx formula
1 parent 7dd5548 commit f66b26c

File tree

5 files changed

+94
-68
lines changed

5 files changed

+94
-68
lines changed

‎workflow/rules/align.smk

Lines changed: 16 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -9,21 +9,27 @@ if config["aligner"] == "bowtie2":
99
index=temp("{}results/bam/{{id}}.tmp.bam.bai".format(outdir)),
1010
threads: config["threads"]["bowtie2"]
1111
params:
12-
index=config["resources"]["ref"]["index"]
13-
if config["resources"]["ref"]["index"] != ""
14-
else "resources/reference_genome/index/index_ref",
12+
index=(
13+
config["resources"]["ref"]["index"]
14+
if config["resources"]["ref"]["index"] != ""
15+
else "resources/reference_genome/index/index_ref"
16+
),
1517
bowtie2=config["params"]["bowtie2"]["global"],
1618
samtools_mem=config["params"]["samtools"]["memory"],
1719
inputsel=(
18-
lambda wildcards, input: " -U {0} ".format(input.reads)
19-
if len(input.reads) == 1
20-
else config["params"]["bowtie2"]["pe"]
21-
+ " -1 {0} -2 {1}".format(*input.reads)
20+
lambda wildcards, input: (
21+
" -U {0} ".format(input.reads)
22+
if len(input.reads) == 1
23+
else config["params"]["bowtie2"]["pe"]
24+
+ " -1 {0} -2 {1}".format(*input.reads)
25+
)
2226
),
2327
forSamblaster=(
24-
lambda wildcards, input: " --ignoreUnmated ".format(input.reads)
25-
if len(input.reads) == 1
26-
else ""
28+
lambda wildcards, input: (
29+
" --ignoreUnmated ".format(input.reads)
30+
if len(input.reads) == 1
31+
else ""
32+
)
2733
),
2834
message:
2935
"Aligning {input} with parameters {params.bowtie2}"

‎workflow/rules/callPeaks.smk

Lines changed: 17 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,9 @@ rule macs2_callNarrowPeak:
1313
log:
1414
"{}results/logs/peakCalling/macs2/{{sample}}_callpeak.log".format(outdir),
1515
params:
16-
isPairedEnd=lambda w: "--format BAM --nomodel"
17-
if is_single_end(w.sample)
18-
else "--format BAMPE",
16+
isPairedEnd=lambda w: (
17+
"--format BAM --nomodel" if is_single_end(w.sample) else "--format BAMPE"
18+
),
1919
gsize=config["params"]["deeptools"]["effective_genome_length"],
2020
qval=config["params"]["peakCalling"]["macs2"]["qvalue"],
2121
otherParams=config["params"]["peakCalling"]["macs2"]["extraOptions"],
@@ -49,11 +49,13 @@ rule macs2_callNormPeaks_narrow:
4949
outdir, sample_to_input[w.sample]
5050
),
5151
logFile="{}results/logs/spike/{{sample}}.normFactor".format(outdir),
52-
logFileInput=lambda wildcards: "{}results/logs/spike/{}.normFactor".format(
53-
outdir, sample_to_input[wildcards.sample]
54-
)
55-
if not pd.isna(sample_to_input[wildcards.sample])
56-
else "{}results/logs/spike/{{sample}}.normFactor".format(outdir),
52+
logFileInput=lambda wildcards: (
53+
"{}results/logs/spike/{}.normFactor".format(
54+
outdir, sample_to_input[wildcards.sample]
55+
)
56+
if not pd.isna(sample_to_input[wildcards.sample])
57+
else "{}results/logs/spike/{{sample}}.normFactor".format(outdir)
58+
),
5759
output:
5860
multiext(
5961
"{}results/peakCallingNorm/{{sample}}".format(outdir),
@@ -101,11 +103,13 @@ rule macs2_callNormPeaks_broad:
101103
outdir, sample_to_input[w.sample]
102104
),
103105
logFile="{}results/logs/spike/{{sample}}.normFactor".format(outdir),
104-
logFileInput=lambda wildcards: "{}results/logs/spike/{}.normFactor".format(
105-
outdir, sample_to_input[wildcards.sample]
106-
)
107-
if not pd.isna(sample_to_input[wildcards.sample])
108-
else "{}results/logs/spike/{{sample}}.normFactor".format(outdir),
106+
logFileInput=lambda wildcards: (
107+
"{}results/logs/spike/{}.normFactor".format(
108+
outdir, sample_to_input[wildcards.sample]
109+
)
110+
if not pd.isna(sample_to_input[wildcards.sample])
111+
else "{}results/logs/spike/{{sample}}.normFactor".format(outdir)
112+
),
109113
output:
110114
multiext(
111115
"{}results/peakCallingNorm/{{sample}}".format(outdir),

‎workflow/rules/differentialPeaks.smk

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,25 @@
11

22
rule differential_peaks:
33
input:
4-
rawReadsOnPeaks="{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
5-
outdir=outdir
6-
)
7-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
8-
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
9-
outdir=outdir
4+
rawReadsOnPeaks=(
5+
"{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
6+
outdir=outdir
7+
)
8+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
9+
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
10+
outdir=outdir
11+
)
1012
),
1113
logFile=get_normFactor_by_antibody,
1214
output:
13-
diffTab="{outdir}results/differentialAnalysis/NormalisedPeaks/{{antibody}}/{{antibody}}_{{contrast}}_diffPeaks.tsv".format(
14-
outdir=outdir
15-
)
16-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
17-
else "{outdir}results/differentialAnalysis/{{antibody}}/{{antibody}}_{{contrast}}_diffPeaks.tsv".format(
18-
outdir=outdir
15+
diffTab=(
16+
"{outdir}results/differentialAnalysis/NormalisedPeaks/{{antibody}}/{{antibody}}_{{contrast}}_diffPeaks.tsv".format(
17+
outdir=outdir
18+
)
19+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
20+
else "{outdir}results/differentialAnalysis/{{antibody}}/{{antibody}}_{{contrast}}_diffPeaks.tsv".format(
21+
outdir=outdir
22+
)
1923
),
2024
params:
2125
padjCutoff=config["diffPeakAnalysis"]["padjust"],

‎workflow/rules/mergeAndAnnotatePeaks.smk

Lines changed: 43 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,33 @@
11

22
rule consensus_peaks:
33
input:
4-
narrowCalled=expand(
5-
"{}results/peakCallingNorm/{{sample}}_narrowPeaks.narrowPeak".format(
6-
outdir
7-
),
8-
sample=narrowSamples,
9-
)
10-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
11-
else expand(
12-
"{}results/peakCalling/macs2/{{sample}}_peaks.narrowPeak".format(outdir),
13-
sample=narrowSamples,
14-
),
15-
broadCalled=expand(
16-
"{}results/peakCallingNorm/{{sample}}_broadPeaks.broadPeak".format(outdir),
17-
sample=broadSamples,
18-
)
19-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
20-
else expand(
21-
"{}results/peakCalling/epic2/{{sample}}_broadPeaks.bed".format(outdir),
22-
sample=broadSamples,
4+
narrowCalled=(
5+
expand(
6+
"{}results/peakCallingNorm/{{sample}}_narrowPeaks.narrowPeak".format(
7+
outdir
8+
),
9+
sample=narrowSamples,
10+
)
11+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
12+
else expand(
13+
"{}results/peakCalling/macs2/{{sample}}_peaks.narrowPeak".format(
14+
outdir
15+
),
16+
sample=narrowSamples,
17+
)
18+
),
19+
broadCalled=(
20+
expand(
21+
"{}results/peakCallingNorm/{{sample}}_broadPeaks.broadPeak".format(
22+
outdir
23+
),
24+
sample=broadSamples,
25+
)
26+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
27+
else expand(
28+
"{}results/peakCalling/epic2/{{sample}}_broadPeaks.bed".format(outdir),
29+
sample=broadSamples,
30+
)
2331
),
2432
output:
2533
"{}results/peakCallingNorm/mergedPeaks/{{antibody}}_consensusPeaks.bed".format(
@@ -53,20 +61,24 @@ rule consensus_peaks:
5361
rule count_reads_on_peaks:
5462
input:
5563
bamFiles=get_bams_by_antibody,
56-
consensus_peaks="{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_consensusPeaks.bed".format(
57-
outdir=outdir
58-
)
59-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
60-
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_consensusPeaks.bed".format(
61-
outdir=outdir
64+
consensus_peaks=(
65+
"{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_consensusPeaks.bed".format(
66+
outdir=outdir
67+
)
68+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
69+
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_consensusPeaks.bed".format(
70+
outdir=outdir
71+
)
6272
),
6373
output:
64-
output_tsv="{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
65-
outdir=outdir
66-
)
67-
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
68-
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
69-
outdir=outdir
74+
output_tsv=(
75+
"{outdir}results/peakCallingNorm/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
76+
outdir=outdir
77+
)
78+
if config["diffPeakAnalysis"]["useSpikeinCalledPeaks"]
79+
else "{outdir}results/peakCalling/mergedPeaks/{{antibody}}_readsOnConsensusPeaks.tsv".format(
80+
outdir=outdir
81+
)
7082
),
7183
params:
7284
joined_bams=lambda w, input: ",".join(input.bamFiles),

‎workflow/scripts/computeNormFactors.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,8 @@
5151
) as input_file:
5252
info_input = input_file.read().strip().split("\n")
5353

54-
gamma = int(info_input[2].split(":")[-1]) / int(
55-
info_input[1].split(":")[-1]
54+
gamma = int(info_input[1].split(":")[-1]) / int(
55+
info_input[0].split(":")[-1]
5656
) # ratio spike/samples in input
5757
alpha = gamma / Nspike * 1000000 # normalization factor
5858

0 commit comments

Comments
 (0)