|
10
|
1 """
|
|
|
2 Parameters:
|
|
16
|
3 - fastq_dump_args: arguments to pass to fastq dumping tool
|
|
|
4 - fastq_dump_retries: number of retry attempts before dying
|
|
|
5 - fastq_dump_tool: name of tool to run, e.g. 'fastq-dump' or 'fasterq-dump'
|
|
10
|
6 - star: arguments to pass to STAR aligner
|
|
16
|
7 - salmon_args: arguments to pass to salmon quant
|
|
10
|
8 - unique_qual: minimum MAPQ needed to be counted in unique BW [default: 10]
|
|
16
|
9 - bw_bed: name of BED file to use with bwtool
|
|
|
10 - max_unalign: maximum number of unaligned reads to save per run accession
|
|
10
|
11 """
|
|
|
12
|
|
16
|
13 STEPS = ['download', 'fastq_check', 'align', 'sort',
|
|
|
14 'bamcount',
|
|
|
15 'salmon',
|
|
|
16 'align_unmapped',
|
|
|
17 'extract_jx']
|
|
10
|
18
|
|
16
|
19 FILES = ['sjout.zst', 'fastq_check.tsv.zst',
|
|
|
20 'unmapped.bam', 'unmapped.fastq.zst',
|
|
10
|
21 'bamcount_nonref.csv.zst',
|
|
|
22 'bamcount_auc.tsv',
|
|
|
23 'bamcount_frag.tsv',
|
|
|
24 'Chimeric.out.junction.zst',
|
|
|
25 'all.exon_bw_count.zst', 'unique.exon_bw_count.zst',
|
|
16
|
26 'all.bw.zst', 'unique.bw.zst',
|
|
|
27 'salmon.tsv.zst',
|
|
|
28 'jx_bed.zst',
|
|
|
29 'manifest'] + list(map(lambda x: x + '.log', STEPS))
|
|
10
|
30
|
|
14
|
31 #get rid of the pesky warnings about current paths
|
|
|
32 #and make sure we're not assuming any wrong relative paths
|
|
12
|
33 import os
|
|
|
34 config['output']=os.path.abspath(config['output'])
|
|
13
|
35 config['temp']=os.path.abspath(config['temp'])
|
|
14
|
36 config['exon_bed']=os.path.abspath(config['exon_bed'])
|
|
|
37 config['ref']=os.path.abspath(config['ref'])
|
|
12
|
38
|
|
10
|
39 import subprocess
|
|
|
40 def run_command(cmd_args):
|
|
|
41 cmd_args = ' '.join(cmd_args)
|
|
|
42 try:
|
|
|
43 subprocess.check_call(args=cmd_args, shell=True)
|
|
|
44 except subprocess.CalledProcessError as cpe:
|
|
|
45 sys.stderr.write("error in run_command for command: %s\n" % cmd_args)
|
|
|
46 raise cpe
|
|
|
47
|
|
16
|
48 import re
|
|
32
|
49 STEPS_FILES_FILTER=re.compile(r'(unmapped)|(download)|(salmon)|(extract_jx)|(jx_bed)|(manifest)|(nonref)|(Chimeric)')
|
|
16
|
50 def remove_steps_files():
|
|
|
51 #modify STEP and FILES
|
|
|
52 #so we don't run download or unmapped steps
|
|
|
53 global FILES
|
|
|
54 global STEPS
|
|
|
55 for i in reversed(range(0,len(STEPS))):
|
|
|
56 f = STEPS[i]
|
|
|
57 if STEPS_FILES_FILTER.search(f):
|
|
|
58 del(STEPS[i])
|
|
|
59 for i in reversed(range(0,len(FILES))):
|
|
|
60 f = FILES[i]
|
|
|
61 if STEPS_FILES_FILTER.search(f):
|
|
|
62 del(FILES[i])
|
|
10
|
63
|
|
16
|
64
|
|
11
|
65 FASTQ_PATT=re.compile(r'([^_\.]+)(_(\d+))?\.((fastq)|fq)(\.gz)?$')
|
|
10
|
66 def prep_for_galaxy_run():
|
|
16
|
67 remove_steps_files()
|
|
10
|
68 try:
|
|
|
69 os.mkdir(config['temp'])
|
|
|
70 except OSError as ose:
|
|
|
71 pass
|
|
|
72 fastqs = config['inputs'].split(',')
|
|
|
73 m = FASTQ_PATT.search(fastqs[0])
|
|
11
|
74 run_acc = 'sample'
|
|
|
75 if m is not None:
|
|
|
76 run_acc = m.group(1)
|
|
12
|
77 study_acc = 'study'
|
|
10
|
78 if 'study' in config:
|
|
|
79 study_acc = config['study']
|
|
|
80 genome = 'hg38'
|
|
|
81 if 'genome' in config:
|
|
|
82 genome = config['genome']
|
|
|
83 method = 'sra'
|
|
|
84 # SRR,SRP,genome
|
|
|
85 # e.g. SRR1557855,SRP045778,ce10
|
|
|
86 #create links which will be used in the align step
|
|
|
87 i = 1
|
|
|
88 for f in fastqs:
|
|
|
89 newf = '%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, i)
|
|
15
|
90 if 'compressed' in config:
|
|
|
91 run_command(['zcat',f,'>',newf])
|
|
|
92 else:
|
|
|
93 os.symlink(os.path.abspath(f), newf)
|
|
10
|
94 #create fastq 0
|
|
|
95 if i == 1:
|
|
15
|
96 os.symlink(os.path.abspath(newf), '%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, 0))
|
|
10
|
97 i += 1
|
|
|
98 #create fastq 2 if not paired
|
|
|
99 if i == 2:
|
|
|
100 open('%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, 2), "w").close()
|
|
|
101 return([run_acc, study_acc, genome, method])
|
|
|
102
|
|
|
103
|
|
|
104 def get_accessions(wildcards):
|
|
|
105 """
|
|
|
106 Grouping of SRRs with the same SRP could happen here
|
|
|
107 """
|
|
|
108 #if running under galaxy where the user will input the
|
|
|
109 #FASTQs, study, and genome directly
|
|
|
110 if 'inputs' in config:
|
|
|
111 (run_acc, study_acc, genome, method) = prep_for_galaxy_run()
|
|
|
112 for ext in FILES:
|
|
|
113 yield os.path.join(config['output'], '%s_%s_%s_%s.%s' % (run_acc, study_acc, genome, method, ext))
|
|
|
114 #here to get the make_galaxy_links rule to fire
|
|
|
115 yield os.path.join(config['output'], '%s_%s_%s_%s.done' % (run_acc, study_acc, genome, method))
|
|
|
116 #not running under galaxy, takes a list of accessions
|
|
|
117 else:
|
|
|
118 for fn in config['input'].split():
|
|
|
119 with open(fn, 'r') as fh:
|
|
|
120 for ln in fh:
|
|
|
121 if ln.count(',') < 2:
|
|
|
122 continue
|
|
|
123 toks = ln.rstrip().split(',')
|
|
|
124 assert 3 <= len(toks) <= 4
|
|
|
125 method = 'sra'
|
|
|
126 if len(toks) == 4:
|
|
|
127 method = toks[3]
|
|
|
128 # SRR,SRP,genome
|
|
|
129 # e.g. SRR1557855,SRP045778,ce10
|
|
|
130 for ext in FILES:
|
|
|
131 yield os.path.join(config['output'], '%s_%s_%s_%s.%s' % (toks[0], toks[1], toks[2], method, ext))
|
|
|
132
|
|
|
133 rule all:
|
|
|
134 input:
|
|
|
135 get_accessions
|
|
|
136
|
|
|
137 rule make_manifest:
|
|
|
138 input:
|
|
16
|
139 config['output'] + '/{quad}.salmon.tsv.zst',
|
|
10
|
140 config['output'] + '/{quad}.sjout.zst',
|
|
16
|
141 config['output'] + '/{quad}.jx_bed.zst',
|
|
10
|
142 config['output'] + '/{quad}.Chimeric.out.junction.zst',
|
|
16
|
143 config['output'] + '/{quad}.unmapped.bam',
|
|
|
144 config['output'] + '/{quad}.unmapped.fastq.zst',
|
|
10
|
145 config['output'] + '/{quad}.bamcount_nonref.csv.zst',
|
|
|
146 config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
147 config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
16
|
148 config['output'] + '/{quad}.fastq_check.tsv.zst',
|
|
10
|
149 config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
150 config['output'] + '/{quad}.unique.exon_bw_count.zst',
|
|
16
|
151 config['output'] + '/{quad}.all.bw.zst',
|
|
|
152 config['output'] + '/{quad}.unique.bw.zst',
|
|
|
153 config['output'] + '/{quad}.align.log',
|
|
|
154 config['output'] + '/{quad}.extract_jx.log',
|
|
|
155 config['output'] + '/{quad}.bamcount.log',
|
|
|
156 config['output'] + '/{quad}.align_unmapped.log',
|
|
|
157 config['output'] + '/{quad}.download.log',
|
|
|
158 config['output'] + '/{quad}.fastq_check.log',
|
|
|
159 config['output'] + '/{quad}.sort.log',
|
|
|
160 config['output'] + '/{quad}.salmon.log'
|
|
10
|
161 output:
|
|
|
162 config['output'] + '/{quad}.manifest'
|
|
|
163 params:
|
|
|
164 quad=lambda wildcards: wildcards.quad
|
|
|
165 run:
|
|
|
166 with open(output[0], 'wt') as fh:
|
|
|
167 for fn in FILES:
|
|
|
168 fh.write(params.quad + "." + fn + '\n')
|
|
|
169
|
|
|
170 def galaxy_link_files(op):
|
|
16
|
171 remove_steps_files()
|
|
10
|
172 a = [op + '/' + f for f in FILES]
|
|
|
173 return a
|
|
|
174
|
|
|
175 rule make_galaxy_links:
|
|
|
176 input:
|
|
|
177 config['output'] + '/{quad}.sjout.zst',
|
|
16
|
178 config['output'] + '/{quad}.fastq_check.tsv.zst',
|
|
10
|
179 config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
180 config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
|
181 config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
182 config['output'] + '/{quad}.unique.exon_bw_count.zst',
|
|
16
|
183 config['output'] + '/{quad}.all.bw.zst',
|
|
|
184 config['output'] + '/{quad}.unique.bw.zst',
|
|
|
185 config['output'] + '/{quad}.fastq_check.log',
|
|
|
186 config['output'] + '/{quad}.align.log',
|
|
|
187 config['output'] + '/{quad}.sort.log',
|
|
|
188 config['output'] + '/{quad}.bamcount.log',
|
|
10
|
189 output:
|
|
|
190 config['output'] + '/{quad}.done'
|
|
|
191 params:
|
|
|
192 quad=lambda wildcards: wildcards.quad,
|
|
|
193 out=galaxy_link_files(config['output'])
|
|
|
194 run:
|
|
|
195 for (i,fn) in enumerate(input):
|
|
|
196 os.symlink(os.path.abspath(fn), params.out[i])
|
|
16
|
197 os.symlink(os.path.abspath(input[-1]), output[0])
|
|
10
|
198
|
|
|
199
|
|
|
200 rule bamcount:
|
|
|
201 input:
|
|
|
202 bam=config['temp'] + '/{quad}-sorted.bam',
|
|
|
203 bamidx=config['temp'] + '/{quad}-sorted.bam.bai',
|
|
|
204 bed=lambda wildcards: '%s/%s/gtf/%s' % (config['ref'], wildcards.quad.split('_')[2], config.get('bw_bed', 'exons.bed'))
|
|
|
205 output:
|
|
|
206 nonref=config['output'] + '/{quad}.bamcount_nonref.csv.zst',
|
|
|
207 auc=config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
208 frag=config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
|
209 all_bw=config['output'] + '/{quad}.all.bw.zst',
|
|
|
210 unique_bw=config['output'] + '/{quad}.unique.bw.zst',
|
|
|
211 all_bw_count=config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
212 unique_bw_count=config['output'] + '/{quad}.unique.exon_bw_count.zst'
|
|
|
213 log:
|
|
|
214 config['output'] + '/{quad}.bamcount.log'
|
|
|
215 params:
|
|
|
216 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
217 uniq_qual=config.get('unique_qual', 10)
|
|
|
218 threads: 4
|
|
|
219 shell:
|
|
|
220 """
|
|
|
221 TMP={config[temp]}/{params.srr}_bamcount
|
|
14
|
222 bamcount {input.bam} \
|
|
10
|
223 --threads {threads} \
|
|
|
224 --coverage \
|
|
|
225 --no-head \
|
|
|
226 --require-mdz \
|
|
|
227 --min-unique-qual {params.uniq_qual} \
|
|
|
228 --frag-dist ${{TMP}} \
|
|
|
229 --bigwig ${{TMP}} \
|
|
|
230 --annotation {input.bed} ${{TMP}} \
|
|
|
231 --auc ${{TMP}} \
|
|
|
232 --alts ${{TMP}} \
|
|
|
233 2>&1 | tee -a {log}
|
|
|
234
|
|
|
235 #
|
|
|
236 # --alts
|
|
|
237 #
|
|
|
238
|
|
|
239 (time zstd ${{TMP}}.alts.tsv -o {output.nonref}) 2>&1 | tee -a {log}
|
|
|
240 size=$(wc -c < {output.nonref})
|
|
|
241 echo "COUNT_NonrefSize ${{size}}"
|
|
|
242 rm -f ${{TMP}}.alts.tsv
|
|
|
243
|
|
|
244 #
|
|
|
245 # --auc
|
|
|
246 #
|
|
|
247 mv ${{TMP}}.auc.tsv {output.auc}
|
|
|
248 size=$(wc -c < {output.auc})
|
|
|
249 echo "COUNT_AucSize ${{size}}"
|
|
|
250 rm -f ${{TMP}}.auc.tsv
|
|
|
251
|
|
|
252 #
|
|
|
253 # --frag-dist
|
|
|
254 #
|
|
|
255 mv ${{TMP}}.frags.tsv {output.frag}
|
|
|
256 size=$(wc -c < {output.frag})
|
|
|
257 echo "COUNT_FragDistSize ${{size}}"
|
|
|
258 rm -f ${{TMP}}.frags.tsv
|
|
|
259
|
|
|
260 #
|
|
|
261 # --bigwig
|
|
|
262 #
|
|
|
263
|
|
|
264 (time zstd ${{TMP}}.all.bw -o {output.all_bw}) 2>&1 | tee -a {log}
|
|
|
265 size=$(wc -c < {output.all_bw})
|
|
|
266 echo "COUNT_BwSize ${{size}}"
|
|
|
267 rm -f ${{TMP}}.all.bw
|
|
|
268
|
|
|
269 (time zstd ${{TMP}}.unique.bw -o {output.unique_bw}) 2>&1 | tee -a {log}
|
|
|
270 size=$(wc -c < {output.unique_bw})
|
|
|
271 echo "COUNT_BwSize ${{size}}"
|
|
|
272 rm -f ${{TMP}}.unique.bw
|
|
|
273
|
|
|
274 #
|
|
|
275 # --annotation
|
|
|
276 #
|
|
|
277
|
|
|
278 (time zstd ${{TMP}}.all.tsv -o {output.all_bw_count}) 2>&1 | tee -a {log}
|
|
|
279 size=$(wc -c < {output.all_bw_count})
|
|
|
280 echo "COUNT_BwQuantSize ${{size}}"
|
|
|
281 rm -f ${{TMP}}.all.tsv
|
|
|
282
|
|
|
283 (time zstd ${{TMP}}.unique.tsv -o {output.unique_bw_count}) 2>&1 | tee -a {log}
|
|
|
284 size=$(wc -c < {output.unique_bw_count})
|
|
|
285 echo "COUNT_BwQuantSize ${{size}}"
|
|
|
286 rm -f ${{TMP}}.unique.tsv
|
|
|
287
|
|
|
288 # Check that all temporaries were properly purged
|
|
|
289 set +o pipefail ; num_files=$(ls -d ${{TMP}}* 2>/dev/null | wc -l)
|
|
|
290 if (( $num_files > 0 )) ; then
|
|
|
291 echo "Failed to purge files (ignore . and ..): $(ls -ad ${{TMP}}*)"
|
|
|
292 exit 1
|
|
|
293 fi
|
|
|
294
|
|
|
295 echo "COUNT_BamcountComplete 1"
|
|
|
296 """
|
|
|
297
|
|
16
|
298 rule bw_zstd:
|
|
|
299 input:
|
|
|
300 config['temp'] + '/{prefix}.bw'
|
|
|
301 output:
|
|
|
302 config['output'] + '/{prefix}.bw.zst'
|
|
|
303 shell:
|
|
|
304 """
|
|
|
305 zstd {input} -o {output}
|
|
|
306 size=$(wc -c < {output})
|
|
|
307 echo "COUNT_BwBytes ${{size}}"
|
|
|
308 echo "COUNT_BwZstdComplete 1"
|
|
|
309 """
|
|
|
310
|
|
10
|
311 rule sort:
|
|
|
312 input:
|
|
|
313 config['temp'] + '/{quad}.bam'
|
|
|
314 wildcard_constraints:
|
|
|
315 quad="[^-]+"
|
|
|
316 output:
|
|
|
317 bam=temp(config['temp'] + '/{quad}-sorted.bam'),
|
|
|
318 bai=temp(config['temp'] + '/{quad}-sorted.bam.bai')
|
|
|
319 log:
|
|
|
320 config['output'] + '/{quad}.sort.log'
|
|
|
321 params:
|
|
|
322 srr=lambda wildcards: wildcards.quad.split('_')[0]
|
|
|
323 threads: 8
|
|
|
324 shell:
|
|
|
325 """
|
|
|
326 TMP="{config[temp]}/sort_temp.{params.srr}"
|
|
|
327 mkdir -p ${{TMP}}
|
|
|
328 time samtools sort \
|
|
|
329 -T ${{TMP}}/samtools_temp \
|
|
|
330 -@ {threads} \
|
|
|
331 -m 64M \
|
|
|
332 -o {output.bam} {input} 2>&1 | tee -a {log}
|
|
|
333 rm -rf ${{TMP}}
|
|
|
334 size=$(wc -c < {output.bam})
|
|
|
335 echo "COUNT_SortedBAMBytes ${{size}}"
|
|
|
336
|
|
|
337 time samtools index -@ {threads} {output.bam} 2>&1 | tee -a {log}
|
|
|
338 echo "COUNT_SortComplete 1"
|
|
|
339 """
|
|
|
340
|
|
16
|
341 rule salmon:
|
|
|
342 input:
|
|
|
343 reads0=config['temp'] + '/{quad}_0.fastq',
|
|
|
344 reads1=config['temp'] + '/{quad}_1.fastq',
|
|
|
345 reads2=config['temp'] + '/{quad}_2.fastq',
|
|
|
346 index1=lambda wildcards: '%s/%s/salmon_index/hash.bin' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
347 index2=lambda wildcards: '%s/%s/salmon_index/sa.bin' % (config['ref'], wildcards.quad.split('_')[2])
|
|
|
348 output:
|
|
|
349 config['output'] + '/{quad}.salmon.tsv.zst'
|
|
|
350 log:
|
|
|
351 config['output'] + '/{quad}.salmon.log'
|
|
|
352 params:
|
|
|
353 index_base=lambda wildcards: '%s/%s/salmon_index' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
354 salmon_args=config.get('salmon_args', '')
|
|
|
355 threads: 8
|
|
|
356 shell:
|
|
|
357 """
|
|
|
358 READ_FILES="-r {input.reads0}"
|
|
|
359 if [[ -s {input.reads2} ]] ; then
|
|
|
360 READ_FILES="-1 {input.reads1} -2 {input.reads2}"
|
|
|
361 fi
|
|
|
362 if set -o pipefail && time salmon quant \
|
|
|
363 --libType U \
|
|
|
364 --quiet \
|
|
|
365 --validateMappings \
|
|
|
366 -i {params.index_base} \
|
|
|
367 -p {threads} \
|
|
|
368 {params.salmon_args} \
|
|
|
369 ${{READ_FILES}} \
|
|
|
370 --output salmon_quant \
|
|
|
371 --minAssignedFrags 1 \
|
|
|
372 2>&1 | tee -a {log}
|
|
|
373 then
|
|
|
374 echo "COUNT_SalmonSuccess 1"
|
|
|
375
|
|
|
376 time zstd salmon_quant/quant.sf -o {output} 2>&1 | tee -a {log}
|
|
|
377 size=$(wc -c < {output})
|
|
|
378 echo "COUNT_SalmonQuantBytes ${{size}}"
|
|
|
379 else
|
|
|
380 touch {output}
|
|
|
381 echo "COUNT_SalmonFailure 1"
|
|
|
382 fi
|
|
|
383
|
|
|
384 rm -rf salmon_quant
|
|
|
385 echo "COUNT_SalmonComplete 1"
|
|
|
386 """
|
|
|
387
|
|
|
388 rule align_unmapped:
|
|
|
389 input:
|
|
|
390 unmapped1=config['temp'] + '/{quad}_1.unmappedfastq',
|
|
|
391 unmapped2=config['temp'] + '/{quad}_2.unmappedfastq',
|
|
|
392 index=lambda wildcards: '%s/%s/unmapped_hisat2_idx/genome.1.ht2' % (config['ref'], wildcards.quad.split('_')[2])
|
|
|
393 output:
|
|
|
394 bam=config['output'] + '/{quad}.unmapped.bam',
|
|
|
395 sample=config['output'] + '/{quad}.unmapped.fastq.zst'
|
|
|
396 log:
|
|
|
397 config['output'] + '/{quad}.align_unmapped.log'
|
|
|
398 params:
|
|
|
399 index_base=lambda wildcards: '%s/%s/unmapped_hisat2_idx/genome' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
400 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
401 hisat2_params=config.get('hisat2', ''),
|
|
|
402 max_unalign=config.get('max_unalign', 100000)
|
|
|
403 threads: 16
|
|
|
404 shell:
|
|
|
405 """
|
|
|
406 TMP="{config[temp]}/align_unmapped_temp.{params.srr}"
|
|
|
407 READ_FILES="-1 {input.unmapped1} -2 {input.unmapped2}"
|
|
|
408 if [[ ! -s {input.unmapped2} ]] ; then
|
|
|
409 READ_FILES="-U {input.unmapped1}"
|
|
|
410 fi
|
|
|
411 time hisat2 \
|
|
|
412 $READ_FILES \
|
|
|
413 -t --mm \
|
|
|
414 -x {params.index_base} \
|
|
|
415 --threads {threads} \
|
|
|
416 {params.hisat2_params} \
|
|
|
417 --un ${{TMP}}.un \
|
|
|
418 --un-conc ${{TMP}}.un_conc \
|
|
|
419 -S ${{TMP}}.sam \
|
|
|
420 2>&1 | tee -a {log}
|
|
|
421
|
|
|
422 #
|
|
|
423 # Make BAM file out of aligned reads
|
|
|
424 #
|
|
|
425 time samtools view \
|
|
|
426 -b -F 4 \
|
|
|
427 -o {output.bam} \
|
|
|
428 ${{TMP}}.sam \
|
|
|
429 2>&1 | tee -a {log}
|
|
|
430 rm -f ${{TMP}}.sam
|
|
|
431
|
|
|
432 #
|
|
|
433 # Save a subset of the doubly unaligned reads
|
|
|
434 #
|
|
|
435 if [[ ! -s {input.unmapped2} ]] ; then
|
|
|
436 # unpaired
|
|
|
437 fastq-sample \
|
|
|
438 -n {params.max_unalign} \
|
|
|
439 -o ${{TMP}}.samp \
|
|
|
440 ${{TMP}}.un
|
|
|
441 rm -f ${{TMP}}.un
|
|
|
442 else
|
|
|
443 # paired-end
|
|
|
444 fastq-sample \
|
|
|
445 -n {params.max_unalign} \
|
|
|
446 -o ${{TMP}}.samp \
|
|
|
447 ${{TMP}}.1.un_conc ${{TMP}}.2.un_conc
|
|
|
448 rm -f ${{TMP}}.1.un_conc ${{TMP}}.2.un_conc
|
|
|
449
|
|
|
450 # interleave the two output fastqs into single file
|
|
|
451 paste ${{TMP}}.samp.1.fastq ${{TMP}}.samp.2.fastq \
|
|
|
452 | paste - - - - \
|
|
|
453 | awk -v OFS="\n" -v FS="\t" '{{print($1,$3,$5,$7,$2,$4,$6,$8)}}' \
|
|
|
454 > ${{TMP}}.samp.fastq
|
|
|
455 rm -f ${{TMP}}.samp.1.fastq ${{TMP}}.samp.2.fastq
|
|
|
456 fi
|
|
|
457
|
|
|
458 test -f ${{TMP}}.samp.fastq
|
|
|
459 time zstd ${{TMP}}.samp.fastq -o {output.sample} 2>&1 | tee -a {log}
|
|
|
460 size=$(wc -c < {output.sample})
|
|
|
461 echo "COUNT_UnmappedSampleBytes ${{size}}"
|
|
|
462 rm -f ${{TMP}}.samp.fastq
|
|
|
463
|
|
|
464 size=$(wc -c < {output.bam})
|
|
|
465 echo "COUNT_UnmappedBamBytes ${{size}}"
|
|
|
466
|
|
|
467 echo "COUNT_AlignUnmappedComplete 1"
|
|
|
468 """
|
|
|
469
|
|
|
470 rule extract_jx:
|
|
|
471 input:
|
|
|
472 bam=config['temp'] + '/{quad}-sorted.bam',
|
|
|
473 bamidx=config['temp'] + '/{quad}-sorted.bam.bai',
|
|
|
474 fa=lambda wildcards: '%s/%s/fasta/genome.fa' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
475 gtf=lambda wildcards: '%s/%s/gtf/genes.gtf' % (config['ref'], wildcards.quad.split('_')[2])
|
|
|
476 output:
|
|
|
477 config['output'] + '/{quad}.jx_bed.zst'
|
|
|
478 params:
|
|
|
479 srr=lambda wildcards: wildcards.quad.split('_')[0]
|
|
|
480 log:
|
|
|
481 config['output'] + '/{quad}.extract_jx.log'
|
|
|
482 shell:
|
|
|
483 """
|
|
|
484 TMP="{config[temp]}/extract_jx.{params.srr}"
|
|
|
485 time regtools junctions extract \
|
|
|
486 -i 20 -a 1 \
|
|
|
487 -o ${{TMP}}.jx_tmp \
|
|
|
488 {input.bam} 2>&1 | tee -a {log}
|
|
|
489 time zstd ${{TMP}}.jx_tmp -o {output} 2>&1 | tee -a {log}
|
|
|
490 rm -f ${{TMP}}.jx_tmp
|
|
|
491
|
|
|
492 size=$(wc -c < {output})
|
|
|
493 echo "COUNT_ExtractJxBytes ${{size}}"
|
|
|
494
|
|
|
495 echo "COUNT_ExtractJxComplete 1"
|
|
|
496 """
|
|
|
497
|
|
10
|
498 rule align:
|
|
|
499 input:
|
|
|
500 reads0=config['temp'] + '/{quad}_0.fastq',
|
|
|
501 reads1=config['temp'] + '/{quad}_1.fastq',
|
|
|
502 reads2=config['temp'] + '/{quad}_2.fastq',
|
|
|
503 index1=lambda wildcards: '%s/%s/star_idx/SAindex' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
504 index2=lambda wildcards: '%s/%s/star_idx/SA' % (config['ref'], wildcards.quad.split('_')[2])
|
|
|
505 wildcard_constraints:
|
|
|
506 quad="[^-]+"
|
|
|
507 output:
|
|
|
508 bam=temp(config['temp'] + '/{quad}.bam'),
|
|
|
509 jxs=config['output'] + '/{quad}.sjout.zst',
|
|
|
510 chimeric=config['output'] + '/{quad}.Chimeric.out.junction.zst',
|
|
|
511 unmapped1=config['temp'] + '/{quad}_1.unmappedfastq',
|
|
|
512 unmapped2=config['temp'] + '/{quad}_2.unmappedfastq'
|
|
|
513 log:
|
|
|
514 config['output'] + '/{quad}.align.log'
|
|
|
515 params:
|
|
|
516 index_base=lambda wildcards: '%s/%s/star_idx' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
517 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
518 star_params="%s %s" % (config.get('star', ''), '--genomeLoad LoadAndRemove' if 'inputs' not in config else '')
|
|
|
519 threads: 16
|
|
|
520 shell:
|
|
|
521 """
|
|
|
522 READ_FILES="{input.reads0}"
|
|
|
523 if [[ -s {input.reads2} ]] ; then
|
|
|
524 READ_FILES="{input.reads1} {input.reads2}"
|
|
|
525 fi
|
|
|
526 TMP="{config[temp]}/align_temp.{params.srr}"
|
|
|
527 rm -rf ${{TMP}}
|
|
|
528 time STAR \
|
|
|
529 {params.star_params} \
|
|
|
530 --runMode alignReads \
|
|
|
531 --runThreadN {threads} \
|
|
|
532 --genomeDir {params.index_base} \
|
|
|
533 --readFilesIn ${{READ_FILES}} \
|
|
|
534 --twopassMode None \
|
|
|
535 --outTmpDir ${{TMP}} \
|
|
|
536 --outReadsUnmapped Fastx \
|
|
|
537 --outMultimapperOrder Random \
|
|
|
538 --outSAMreadID Number \
|
|
|
539 --outSAMtype BAM Unsorted \
|
|
|
540 --outSAMmode NoQS \
|
|
|
541 --outSAMattributes NH MD \
|
|
|
542 --chimOutType Junctions \
|
|
|
543 --chimOutJunctionFormat 1 \
|
|
|
544 --chimSegmentReadGapMax 3 \
|
|
|
545 --chimJunctionOverhangMin 12 \
|
|
|
546 --chimSegmentMin 12 2>&1 | tee -a {log}
|
|
|
547
|
|
|
548 # Full set of output files:
|
|
|
549 #
|
|
|
550 # Aligned.out.bam
|
|
|
551 # Chimeric.out.junction
|
|
|
552 # Log.final.out
|
|
|
553 # Log.out
|
|
|
554 # Log.progress.out
|
|
|
555 # SJ.out.tab
|
|
|
556 # Unmapped.out.mate1
|
|
|
557 # Unmapped.out.mate2 (if any reads were paired-end)
|
|
|
558
|
|
|
559 #
|
|
|
560 # Logs
|
|
|
561 #
|
|
|
562 rm -rf ${{TMP}}
|
|
|
563 cat Log.out >> {log}
|
|
|
564 cat Log.final.out >> {log}
|
|
|
565 rm -f Log*.out
|
|
|
566
|
|
|
567 #
|
|
|
568 # Junctions
|
|
|
569 #
|
|
|
570 test -f SJ.out.tab
|
|
|
571 time zstd SJ.out.tab -o {output.jxs} 2>&1 | tee -a {log}
|
|
|
572 rm -f SJ.out.tab
|
|
|
573 size=$(wc -c < {output.jxs})
|
|
|
574 echo "COUNT_CompressedJxBytes ${{size}}"
|
|
|
575
|
|
|
576 #
|
|
|
577 # Chimerics
|
|
|
578 #
|
|
|
579 test -f Chimeric.out.junction
|
|
|
580 test -s Chimeric.out.junction
|
|
|
581 sort -k1,1 -n -k2,2 Chimeric.out.junction > Chimeric.out.junction.sorted
|
|
|
582 time zstd Chimeric.out.junction.sorted -o {output.chimeric} 2>&1 | tee -a {log}
|
|
|
583 rm -f Chimeric.out.junction Chimeric.out.junction.sorted
|
|
|
584 size=$(wc -c < {output.chimeric})
|
|
|
585 echo "COUNT_ChimericBytes ${{size}}"
|
|
|
586
|
|
|
587 #
|
|
|
588 # Unmapped
|
|
|
589 #
|
|
|
590 touch {output.unmapped2}
|
|
|
591 test -f Unmapped.out.mate1
|
|
|
592 mv Unmapped.out.mate1 {output.unmapped1}
|
|
|
593 if [[ -f Unmapped.out.mate2 ]] ; then
|
|
|
594 mv Unmapped.out.mate2 {output.unmapped2}
|
|
|
595 fi
|
|
|
596
|
|
|
597 #
|
|
|
598 # Alignments
|
|
|
599 #
|
|
|
600 size=$(wc -c < Aligned.out.bam)
|
|
|
601 echo "COUNT_BAMBytes ${{size}}"
|
|
|
602 mv Aligned.out.bam {output.bam}
|
|
16
|
603
|
|
10
|
604 echo "COUNT_AlignComplete 1"
|
|
|
605 """
|
|
16
|
606
|
|
|
607 rule fastq_check:
|
|
|
608 input:
|
|
|
609 reads0=config['temp'] + '/{quad}_0.fastq',
|
|
|
610 reads1=config['temp'] + '/{quad}_1.fastq',
|
|
|
611 reads2=config['temp'] + '/{quad}_2.fastq'
|
|
|
612 output:
|
|
|
613 config['output'] + '/{quad}.fastq_check.tsv.zst'
|
|
|
614 log:
|
|
|
615 config['output'] + '/{quad}.fastq_check.log'
|
|
|
616 params:
|
|
|
617 srr=lambda wildcards: wildcards.quad.split('_')[0]
|
|
|
618 shell:
|
|
|
619 """
|
|
|
620 TMP="{config[temp]}/fastq_check-{params.srr}.tsv"
|
|
|
621 touch ${{TMP}}
|
|
|
622 if [[ -s {input.reads0} ]] ; then
|
|
|
623 time seqtk fqchk -q0 {input.reads0} >>${{TMP}} 2>>{log}
|
|
|
624 fi
|
|
|
625 if [[ -s {input.reads1} ]] ; then
|
|
|
626 time seqtk fqchk -q0 {input.reads1} >>${{TMP}} 2>>{log}
|
|
|
627 fi
|
|
|
628 if [[ -s {input.reads2} ]] ; then
|
|
|
629 time seqtk fqchk -q0 {input.reads2} >>${{TMP}} 2>>{log}
|
|
|
630 fi
|
|
|
631 time zstd ${{TMP}} -o {output} 2>&1 | tee -a {log}
|
|
|
632 size=$(wc -c < {output})
|
|
|
633 echo "COUNT_FastqCheckBytes ${{size}}"
|
|
|
634
|
|
|
635 echo "COUNT_FastqCheckComplete 1"
|
|
|
636 """
|
|
|
637
|
|
|
638 rule download:
|
|
|
639 output:
|
|
|
640 temp(config['temp'] + '/{quad}_0.fastq'),
|
|
|
641 temp(config['temp'] + '/{quad}_1.fastq'),
|
|
|
642 temp(config['temp'] + '/{quad}_2.fastq')
|
|
|
643 log:
|
|
|
644 config['output'] + '/{quad}.download.log'
|
|
|
645 params:
|
|
|
646 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
647 method=lambda wildcards: wildcards.quad.split('_')[3],
|
|
|
648 fd_params=config.get('fastq_dump_args', ''),
|
|
|
649 retries=config.get('fastq_dump_retries', '5')
|
|
|
650 threads: 4
|
|
|
651 shell:
|
|
|
652 """
|
|
|
653 set -x
|
|
|
654 SUCCESS=0
|
|
|
655 TIMEOUT=10
|
|
|
656 PARAMS=""
|
|
|
657 TMP="{config[temp]}/dl-{params.srr}"
|
|
|
658 ! test -d ${{TMP}}
|
|
|
659
|
|
|
660 if [[ {params.method} == "sra" ]] ; then
|
|
|
661 USE_FASTERQ=1
|
|
|
662 for i in {{ 1..{params.retries} }} ; do
|
|
|
663 if time prefetch -t fasp -O ${{TMP}} -L info {params.srr} 2>&1 >> {log} ; then
|
|
|
664 SUCCESS=1
|
|
|
665 break
|
|
|
666 else
|
|
|
667 echo "COUNT_SraRetries 1"
|
|
|
668 TIMEOUT=$((${{TIMEOUT}} * 2))
|
|
|
669 sleep ${{TIMEOUT}}
|
|
|
670 fi
|
|
|
671 done
|
|
|
672 if (( $SUCCESS == 0 )) ; then
|
|
|
673 echo "COUNT_SraFailures 1"
|
|
|
674 exit 1
|
|
|
675 fi
|
|
|
676 test -f ${{TMP}}/*.sra
|
|
|
677 size=$(cat ${{TMP}}/*.sra | wc -c)
|
|
|
678 echo "COUNT_SraBytesDownloaded ${{size}}"
|
|
|
679 if (( ${{USE_FASTERQ}} == 1 )) ; then
|
|
|
680 time fasterq-dump ${{TMP}}/*.sra \
|
|
|
681 -e {threads} \
|
|
|
682 -t ${{TMP}} \
|
|
|
683 -L info \
|
|
|
684 --split-files \
|
|
|
685 --skip-technical \
|
|
|
686 -o {params.srr}.fastq \
|
|
|
687 2>&1 >> {log}
|
|
|
688 test -f {params.srr}_2.fastq || mv {params.srr}.fastq {params.srr}_0.fastq
|
|
|
689 else
|
|
|
690 time fastq-dump ${{TMP}}/*.sra \
|
|
|
691 -L info \
|
|
|
692 --split-files \
|
|
|
693 --skip-technical \
|
|
|
694 -O . \
|
|
|
695 2>&1 >> {log}
|
|
|
696 test -f {params.srr}_2.fastq || mv {params.srr}_1.fastq {params.srr}_0.fastq
|
|
|
697 fi
|
|
|
698 rm -rf ${{TMP}}
|
|
|
699 elif [[ {params.method} == "gdc" ]] ; then
|
|
|
700 TOKEN=~/gdc/latest.txt
|
|
|
701 if [[ ! -f ${{TOKEN}} ]] ; then
|
|
|
702 echo "ERROR: no GDC token file found at ${{TOKEN}}"
|
|
|
703 exit 1
|
|
|
704 fi
|
|
|
705 for i in {{ 1..{params.retries} }} ; do
|
|
|
706 if time gdc-client download \
|
|
|
707 -t ${{TOKEN}} \
|
|
|
708 --log-file {{TMP}}/log.txt \
|
|
|
709 -d {{TEMP}} \
|
|
|
710 {params.srr} 2>&1 >> {log}
|
|
|
711 then
|
|
|
712 SUCCESS=1
|
|
|
713 break
|
|
|
714 else
|
|
|
715 echo "COUNT_GdcRetries 1"
|
|
|
716 TIMEOUT=$((${{TIMEOUT}} * 2))
|
|
|
717 sleep ${{TIMEOUT}}
|
|
|
718 fi
|
|
|
719 done
|
|
|
720 if (( $SUCCESS == 0 )) ; then
|
|
|
721 echo "COUNT_GdcFailures 1"
|
|
|
722 exit 1
|
|
|
723 fi
|
|
|
724 test -d {{TEMP}}/{params.srr}
|
|
|
725 test -f {{TEMP}}/{params.srr}/*.tar.gz
|
|
|
726
|
|
|
727 echo "=== gdc-client log.txt begin ===" >> {log}
|
|
|
728 cat {{TMP}}/log.txt >> {log}
|
|
|
729 echo "=== gdc-client log.txt end===" >> {log}
|
|
|
730
|
|
|
731 size=$(cat {{TEMP}}/{params.srr}/*.tar.gz | wc -c)
|
|
|
732 echo "COUNT_GdcBytesDownloaded ${{size}}"
|
|
|
733
|
|
|
734 tar zxvf {{TEMP}}/{params.srr}/*.tar.gz
|
|
|
735 rm -rf {{TEMP}}
|
|
|
736
|
|
|
737 num_1s=$(ls -1 *_1.fastq | wc -l)
|
|
|
738 num_2s=$(ls -1 *_2.fastq | wc -l)
|
|
|
739 if (( ${{num_1s}} == 0 )) ; then
|
|
|
740 echo "ERROR: No _1.fastq files output"
|
|
|
741 exit 1
|
|
|
742 fi
|
|
|
743 if (( ${{num_1s}} > 1 )) ; then
|
|
|
744 echo "ERROR: More than one _1.fastq file found"
|
|
|
745 exit 1
|
|
|
746 fi
|
|
|
747 if (( ${{num_2s}} == 0 )) ; then
|
|
|
748 # unpaired
|
|
|
749 mv *_1.fastq {params.srr}_0.fastq
|
|
|
750 else
|
|
|
751 # paired-end
|
|
|
752 mv *_1.fastq {params.srr}_1.fastq
|
|
|
753 mv *_2.fastq {params.srr}_2.fastq
|
|
|
754 fi
|
|
|
755 fi
|
|
|
756
|
|
|
757 # Next chunk expects the FASTQ files to exist in the current directory
|
|
|
758 # named {{params.srr}}_{{0,1,2}}.fastq
|
|
|
759 size=0
|
|
|
760 for i in {{0..2}} ; do
|
|
|
761 fn={params.srr}_${{i}}.fastq
|
|
|
762 if [[ -f ${{fn}} ]] ; then
|
|
|
763 echo "COUNT_FilesDownloaded 1"
|
|
|
764 else
|
|
|
765 touch ${{fn}}
|
|
|
766 fi
|
|
|
767 size=$((${{size}} + $(wc -c < ${{fn}})))
|
|
|
768 mv ${{fn}} {config[temp]}/{wildcards.quad}_${{i}}.fastq
|
|
|
769 done
|
|
|
770 echo "COUNT_BytesDownloaded ${{size}}"
|
|
|
771 echo "COUNT_DownloadComplete 1"
|
|
|
772 """
|
