|
10
|
1 """
|
|
|
2 Parameters:
|
|
|
3 - star: arguments to pass to STAR aligner
|
|
|
4 - unique_qual: minimum MAPQ needed to be counted in unique BW [default: 10]
|
|
|
5 """
|
|
|
6
|
|
|
7 STEPS = ['align', 'sort', 'bamcount']
|
|
|
8
|
|
|
9 FILES = ['sjout.zst',
|
|
|
10 'bamcount_nonref.csv.zst',
|
|
|
11 'bamcount_auc.tsv',
|
|
|
12 'bamcount_frag.tsv',
|
|
|
13 'Chimeric.out.junction.zst',
|
|
|
14 'all.exon_bw_count.zst', 'unique.exon_bw_count.zst',
|
|
|
15 'manifest']
|
|
|
16
|
|
12
|
17 import os
|
|
|
18 config['output']=os.path.abspath(config['output'])
|
|
|
19
|
|
10
|
20 import subprocess
|
|
|
21 def run_command(cmd_args):
|
|
|
22 cmd_args = ' '.join(cmd_args)
|
|
|
23 try:
|
|
|
24 subprocess.check_call(args=cmd_args, shell=True)
|
|
|
25 except subprocess.CalledProcessError as cpe:
|
|
|
26 sys.stderr.write("error in run_command for command: %s\n" % cmd_args)
|
|
|
27 raise cpe
|
|
|
28
|
|
|
29
|
|
|
30 import re
|
|
11
|
31 FASTQ_PATT=re.compile(r'([^_\.]+)(_(\d+))?\.((fastq)|fq)(\.gz)?$')
|
|
10
|
32 def prep_for_galaxy_run():
|
|
|
33 try:
|
|
|
34 os.mkdir(config['temp'])
|
|
|
35 except OSError as ose:
|
|
|
36 pass
|
|
|
37 fastqs = config['inputs'].split(',')
|
|
|
38 m = FASTQ_PATT.search(fastqs[0])
|
|
11
|
39 run_acc = 'sample'
|
|
|
40 if m is not None:
|
|
|
41 run_acc = m.group(1)
|
|
12
|
42 study_acc = 'study'
|
|
10
|
43 if 'study' in config:
|
|
|
44 study_acc = config['study']
|
|
|
45 genome = 'hg38'
|
|
|
46 if 'genome' in config:
|
|
|
47 genome = config['genome']
|
|
|
48 method = 'sra'
|
|
|
49 # SRR,SRP,genome
|
|
|
50 # e.g. SRR1557855,SRP045778,ce10
|
|
|
51 #create links which will be used in the align step
|
|
|
52 i = 1
|
|
|
53 for f in fastqs:
|
|
|
54 newf = '%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, i)
|
|
|
55 run_command(['zcat',f,'>',newf])
|
|
|
56 #create fastq 0
|
|
|
57 if i == 1:
|
|
|
58 try:
|
|
|
59 os.symlink(os.path.abspath(newf), '%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, 0))
|
|
|
60 except FileExistsError as fee:
|
|
|
61 pass
|
|
|
62 i += 1
|
|
|
63 #create fastq 2 if not paired
|
|
|
64 if i == 2:
|
|
|
65 open('%s/%s_%s_%s_%s_%d.fastq' % (config['temp'], run_acc, study_acc, genome, method, 2), "w").close()
|
|
|
66 #create expected file structure for annotated exon bed file & reference index
|
|
|
67 ref = config['ref']
|
|
|
68 config['ref'] = '.'
|
|
|
69 os.makedirs('%s/%s' % (config['ref'], genome))
|
|
|
70 os.symlink(ref, '%s/%s/star_idx' % (config['ref'], genome))
|
|
|
71 os.makedirs('%s/%s/gtf' % (config['ref'], genome))
|
|
|
72 os.symlink(config['exon_bed'], 'exons.tmp')
|
|
|
73 os.symlink('../../exons.tmp', '%s/%s/gtf/%s' % (config['ref'], genome, config.get('bw_bed', 'exons.bed')))
|
|
|
74 return([run_acc, study_acc, genome, method])
|
|
|
75
|
|
|
76
|
|
|
77 def get_accessions(wildcards):
|
|
|
78 """
|
|
|
79 Grouping of SRRs with the same SRP could happen here
|
|
|
80 """
|
|
|
81 #if running under galaxy where the user will input the
|
|
|
82 #FASTQs, study, and genome directly
|
|
|
83 if 'inputs' in config:
|
|
|
84 (run_acc, study_acc, genome, method) = prep_for_galaxy_run()
|
|
|
85 for ext in FILES:
|
|
|
86 yield os.path.join(config['output'], '%s_%s_%s_%s.%s' % (run_acc, study_acc, genome, method, ext))
|
|
|
87 #here to get the make_galaxy_links rule to fire
|
|
|
88 yield os.path.join(config['output'], '%s_%s_%s_%s.done' % (run_acc, study_acc, genome, method))
|
|
|
89 #not running under galaxy, takes a list of accessions
|
|
|
90 else:
|
|
|
91 for fn in config['input'].split():
|
|
|
92 with open(fn, 'r') as fh:
|
|
|
93 for ln in fh:
|
|
|
94 if ln.count(',') < 2:
|
|
|
95 continue
|
|
|
96 toks = ln.rstrip().split(',')
|
|
|
97 assert 3 <= len(toks) <= 4
|
|
|
98 method = 'sra'
|
|
|
99 if len(toks) == 4:
|
|
|
100 method = toks[3]
|
|
|
101 # SRR,SRP,genome
|
|
|
102 # e.g. SRR1557855,SRP045778,ce10
|
|
|
103 for ext in FILES:
|
|
|
104 yield os.path.join(config['output'], '%s_%s_%s_%s.%s' % (toks[0], toks[1], toks[2], method, ext))
|
|
|
105
|
|
|
106 rule all:
|
|
|
107 input:
|
|
|
108 get_accessions
|
|
|
109
|
|
|
110
|
|
|
111 rule make_manifest:
|
|
|
112 input:
|
|
|
113 config['output'] + '/{quad}.sjout.zst',
|
|
|
114 config['output'] + '/{quad}.Chimeric.out.junction.zst',
|
|
|
115 config['output'] + '/{quad}.bamcount_nonref.csv.zst',
|
|
|
116 config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
117 config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
|
118 config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
119 config['output'] + '/{quad}.unique.exon_bw_count.zst',
|
|
|
120 output:
|
|
|
121 config['output'] + '/{quad}.manifest'
|
|
|
122 params:
|
|
|
123 quad=lambda wildcards: wildcards.quad
|
|
|
124 run:
|
|
|
125 with open(output[0], 'wt') as fh:
|
|
|
126 for fn in FILES:
|
|
|
127 fh.write(params.quad + "." + fn + '\n')
|
|
|
128
|
|
|
129 def galaxy_link_files(op):
|
|
|
130 a = [op + '/' + f for f in FILES]
|
|
|
131 a.extend([op + '/align.log', op + '/bamcount.log'])
|
|
|
132 return a
|
|
|
133
|
|
|
134 rule make_galaxy_links:
|
|
|
135 input:
|
|
|
136 config['output'] + '/{quad}.sjout.zst',
|
|
|
137 config['output'] + '/{quad}.bamcount_nonref.csv.zst',
|
|
|
138 config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
139 config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
|
140 config['output'] + '/{quad}.Chimeric.out.junction.zst',
|
|
|
141 config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
142 config['output'] + '/{quad}.unique.exon_bw_count.zst',
|
|
|
143 config['output'] + '/{quad}.manifest'
|
|
|
144 output:
|
|
|
145 config['output'] + '/{quad}.done'
|
|
|
146 params:
|
|
|
147 quad=lambda wildcards: wildcards.quad,
|
|
|
148 out=galaxy_link_files(config['output'])
|
|
|
149 run:
|
|
|
150 inputs = input.extend([config['output'] + '/' + params.quad + '.align.log',
|
|
|
151 config['output'] + '/' + params.quad + '.bamcount.log'])
|
|
|
152 for (i,fn) in enumerate(input):
|
|
|
153 os.symlink(os.path.abspath(fn), params.out[i])
|
|
|
154 os.symlink(os.path.abspath(input[-3]), output[0])
|
|
|
155
|
|
|
156
|
|
|
157 rule bamcount:
|
|
|
158 input:
|
|
|
159 bam=config['temp'] + '/{quad}-sorted.bam',
|
|
|
160 bamidx=config['temp'] + '/{quad}-sorted.bam.bai',
|
|
|
161 #exe='/bamcount/bamcount',
|
|
|
162 exe='bamcount',
|
|
|
163 bed=lambda wildcards: '%s/%s/gtf/%s' % (config['ref'], wildcards.quad.split('_')[2], config.get('bw_bed', 'exons.bed'))
|
|
|
164 output:
|
|
|
165 nonref=config['output'] + '/{quad}.bamcount_nonref.csv.zst',
|
|
|
166 auc=config['output'] + '/{quad}.bamcount_auc.tsv',
|
|
|
167 frag=config['output'] + '/{quad}.bamcount_frag.tsv',
|
|
|
168 all_bw=config['output'] + '/{quad}.all.bw.zst',
|
|
|
169 unique_bw=config['output'] + '/{quad}.unique.bw.zst',
|
|
|
170 all_bw_count=config['output'] + '/{quad}.all.exon_bw_count.zst',
|
|
|
171 unique_bw_count=config['output'] + '/{quad}.unique.exon_bw_count.zst'
|
|
|
172 log:
|
|
|
173 config['output'] + '/{quad}.bamcount.log'
|
|
|
174 params:
|
|
|
175 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
176 uniq_qual=config.get('unique_qual', 10)
|
|
|
177 threads: 4
|
|
|
178 shell:
|
|
|
179 """
|
|
|
180 TMP={config[temp]}/{params.srr}_bamcount
|
|
|
181 {input.exe} {input.bam} \
|
|
|
182 --threads {threads} \
|
|
|
183 --coverage \
|
|
|
184 --no-head \
|
|
|
185 --require-mdz \
|
|
|
186 --min-unique-qual {params.uniq_qual} \
|
|
|
187 --frag-dist ${{TMP}} \
|
|
|
188 --bigwig ${{TMP}} \
|
|
|
189 --annotation {input.bed} ${{TMP}} \
|
|
|
190 --auc ${{TMP}} \
|
|
|
191 --alts ${{TMP}} \
|
|
|
192 2>&1 | tee -a {log}
|
|
|
193
|
|
|
194 #
|
|
|
195 # --alts
|
|
|
196 #
|
|
|
197
|
|
|
198 (time zstd ${{TMP}}.alts.tsv -o {output.nonref}) 2>&1 | tee -a {log}
|
|
|
199 size=$(wc -c < {output.nonref})
|
|
|
200 echo "COUNT_NonrefSize ${{size}}"
|
|
|
201 rm -f ${{TMP}}.alts.tsv
|
|
|
202
|
|
|
203 #
|
|
|
204 # --auc
|
|
|
205 #
|
|
|
206 mv ${{TMP}}.auc.tsv {output.auc}
|
|
|
207 size=$(wc -c < {output.auc})
|
|
|
208 echo "COUNT_AucSize ${{size}}"
|
|
|
209 rm -f ${{TMP}}.auc.tsv
|
|
|
210
|
|
|
211 #
|
|
|
212 # --frag-dist
|
|
|
213 #
|
|
|
214 mv ${{TMP}}.frags.tsv {output.frag}
|
|
|
215 size=$(wc -c < {output.frag})
|
|
|
216 echo "COUNT_FragDistSize ${{size}}"
|
|
|
217 rm -f ${{TMP}}.frags.tsv
|
|
|
218
|
|
|
219 #
|
|
|
220 # --bigwig
|
|
|
221 #
|
|
|
222
|
|
|
223 (time zstd ${{TMP}}.all.bw -o {output.all_bw}) 2>&1 | tee -a {log}
|
|
|
224 size=$(wc -c < {output.all_bw})
|
|
|
225 echo "COUNT_BwSize ${{size}}"
|
|
|
226 rm -f ${{TMP}}.all.bw
|
|
|
227
|
|
|
228 (time zstd ${{TMP}}.unique.bw -o {output.unique_bw}) 2>&1 | tee -a {log}
|
|
|
229 size=$(wc -c < {output.unique_bw})
|
|
|
230 echo "COUNT_BwSize ${{size}}"
|
|
|
231 rm -f ${{TMP}}.unique.bw
|
|
|
232
|
|
|
233 #
|
|
|
234 # --annotation
|
|
|
235 #
|
|
|
236
|
|
|
237 (time zstd ${{TMP}}.all.tsv -o {output.all_bw_count}) 2>&1 | tee -a {log}
|
|
|
238 size=$(wc -c < {output.all_bw_count})
|
|
|
239 echo "COUNT_BwQuantSize ${{size}}"
|
|
|
240 rm -f ${{TMP}}.all.tsv
|
|
|
241
|
|
|
242 (time zstd ${{TMP}}.unique.tsv -o {output.unique_bw_count}) 2>&1 | tee -a {log}
|
|
|
243 size=$(wc -c < {output.unique_bw_count})
|
|
|
244 echo "COUNT_BwQuantSize ${{size}}"
|
|
|
245 rm -f ${{TMP}}.unique.tsv
|
|
|
246
|
|
|
247 # Check that all temporaries were properly purged
|
|
|
248 set +o pipefail ; num_files=$(ls -d ${{TMP}}* 2>/dev/null | wc -l)
|
|
|
249 if (( $num_files > 0 )) ; then
|
|
|
250 echo "Failed to purge files (ignore . and ..): $(ls -ad ${{TMP}}*)"
|
|
|
251 exit 1
|
|
|
252 fi
|
|
|
253
|
|
|
254 echo "COUNT_BamcountComplete 1"
|
|
|
255 """
|
|
|
256
|
|
|
257 rule sort:
|
|
|
258 input:
|
|
|
259 config['temp'] + '/{quad}.bam'
|
|
|
260 wildcard_constraints:
|
|
|
261 quad="[^-]+"
|
|
|
262 output:
|
|
|
263 bam=temp(config['temp'] + '/{quad}-sorted.bam'),
|
|
|
264 bai=temp(config['temp'] + '/{quad}-sorted.bam.bai')
|
|
|
265 log:
|
|
|
266 config['output'] + '/{quad}.sort.log'
|
|
|
267 params:
|
|
|
268 srr=lambda wildcards: wildcards.quad.split('_')[0]
|
|
|
269 threads: 8
|
|
|
270 shell:
|
|
|
271 """
|
|
|
272 TMP="{config[temp]}/sort_temp.{params.srr}"
|
|
|
273 mkdir -p ${{TMP}}
|
|
|
274 time samtools sort \
|
|
|
275 -T ${{TMP}}/samtools_temp \
|
|
|
276 -@ {threads} \
|
|
|
277 -m 64M \
|
|
|
278 -o {output.bam} {input} 2>&1 | tee -a {log}
|
|
|
279 rm -rf ${{TMP}}
|
|
|
280 size=$(wc -c < {output.bam})
|
|
|
281 echo "COUNT_SortedBAMBytes ${{size}}"
|
|
|
282
|
|
|
283 time samtools index -@ {threads} {output.bam} 2>&1 | tee -a {log}
|
|
|
284 echo "COUNT_SortComplete 1"
|
|
|
285 """
|
|
|
286
|
|
|
287 rule align:
|
|
|
288 input:
|
|
|
289 reads0=config['temp'] + '/{quad}_0.fastq',
|
|
|
290 reads1=config['temp'] + '/{quad}_1.fastq',
|
|
|
291 reads2=config['temp'] + '/{quad}_2.fastq',
|
|
|
292 index1=lambda wildcards: '%s/%s/star_idx/SAindex' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
293 index2=lambda wildcards: '%s/%s/star_idx/SA' % (config['ref'], wildcards.quad.split('_')[2])
|
|
|
294 wildcard_constraints:
|
|
|
295 quad="[^-]+"
|
|
|
296 output:
|
|
|
297 bam=temp(config['temp'] + '/{quad}.bam'),
|
|
|
298 jxs=config['output'] + '/{quad}.sjout.zst',
|
|
|
299 chimeric=config['output'] + '/{quad}.Chimeric.out.junction.zst',
|
|
|
300 unmapped1=config['temp'] + '/{quad}_1.unmappedfastq',
|
|
|
301 unmapped2=config['temp'] + '/{quad}_2.unmappedfastq'
|
|
|
302 log:
|
|
|
303 config['output'] + '/{quad}.align.log'
|
|
|
304 params:
|
|
|
305 index_base=lambda wildcards: '%s/%s/star_idx' % (config['ref'], wildcards.quad.split('_')[2]),
|
|
|
306 srr=lambda wildcards: wildcards.quad.split('_')[0],
|
|
|
307 star_params="%s %s" % (config.get('star', ''), '--genomeLoad LoadAndRemove' if 'inputs' not in config else '')
|
|
|
308 threads: 16
|
|
|
309 shell:
|
|
|
310 """
|
|
|
311 READ_FILES="{input.reads0}"
|
|
|
312 if [[ -s {input.reads2} ]] ; then
|
|
|
313 READ_FILES="{input.reads1} {input.reads2}"
|
|
|
314 fi
|
|
|
315 TMP="{config[temp]}/align_temp.{params.srr}"
|
|
|
316 rm -rf ${{TMP}}
|
|
|
317 time STAR \
|
|
|
318 {params.star_params} \
|
|
|
319 --runMode alignReads \
|
|
|
320 --runThreadN {threads} \
|
|
|
321 --genomeDir {params.index_base} \
|
|
|
322 --readFilesIn ${{READ_FILES}} \
|
|
|
323 --twopassMode None \
|
|
|
324 --outTmpDir ${{TMP}} \
|
|
|
325 --outReadsUnmapped Fastx \
|
|
|
326 --outMultimapperOrder Random \
|
|
|
327 --outSAMreadID Number \
|
|
|
328 --outSAMtype BAM Unsorted \
|
|
|
329 --outSAMmode NoQS \
|
|
|
330 --outSAMattributes NH MD \
|
|
|
331 --chimOutType Junctions \
|
|
|
332 --chimOutJunctionFormat 1 \
|
|
|
333 --chimSegmentReadGapMax 3 \
|
|
|
334 --chimJunctionOverhangMin 12 \
|
|
|
335 --chimSegmentMin 12 2>&1 | tee -a {log}
|
|
|
336
|
|
|
337 # Full set of output files:
|
|
|
338 #
|
|
|
339 # Aligned.out.bam
|
|
|
340 # Chimeric.out.junction
|
|
|
341 # Log.final.out
|
|
|
342 # Log.out
|
|
|
343 # Log.progress.out
|
|
|
344 # SJ.out.tab
|
|
|
345 # Unmapped.out.mate1
|
|
|
346 # Unmapped.out.mate2 (if any reads were paired-end)
|
|
|
347
|
|
|
348 #
|
|
|
349 # Logs
|
|
|
350 #
|
|
|
351 rm -rf ${{TMP}}
|
|
|
352 cat Log.out >> {log}
|
|
|
353 cat Log.final.out >> {log}
|
|
|
354 rm -f Log*.out
|
|
|
355
|
|
|
356 #
|
|
|
357 # Junctions
|
|
|
358 #
|
|
|
359 test -f SJ.out.tab
|
|
|
360 time zstd SJ.out.tab -o {output.jxs} 2>&1 | tee -a {log}
|
|
|
361 rm -f SJ.out.tab
|
|
|
362 size=$(wc -c < {output.jxs})
|
|
|
363 echo "COUNT_CompressedJxBytes ${{size}}"
|
|
|
364
|
|
|
365 #
|
|
|
366 # Chimerics
|
|
|
367 #
|
|
|
368 test -f Chimeric.out.junction
|
|
|
369 test -s Chimeric.out.junction
|
|
|
370 sort -k1,1 -n -k2,2 Chimeric.out.junction > Chimeric.out.junction.sorted
|
|
|
371 time zstd Chimeric.out.junction.sorted -o {output.chimeric} 2>&1 | tee -a {log}
|
|
|
372 rm -f Chimeric.out.junction Chimeric.out.junction.sorted
|
|
|
373 size=$(wc -c < {output.chimeric})
|
|
|
374 echo "COUNT_ChimericBytes ${{size}}"
|
|
|
375
|
|
|
376 #
|
|
|
377 # Unmapped
|
|
|
378 #
|
|
|
379 touch {output.unmapped2}
|
|
|
380 test -f Unmapped.out.mate1
|
|
|
381 mv Unmapped.out.mate1 {output.unmapped1}
|
|
|
382 if [[ -f Unmapped.out.mate2 ]] ; then
|
|
|
383 mv Unmapped.out.mate2 {output.unmapped2}
|
|
|
384 fi
|
|
|
385
|
|
|
386 #
|
|
|
387 # Alignments
|
|
|
388 #
|
|
|
389 size=$(wc -c < Aligned.out.bam)
|
|
|
390 echo "COUNT_BAMBytes ${{size}}"
|
|
|
391 mv Aligned.out.bam {output.bam}
|
|
|
392 echo "COUNT_AlignComplete 1"
|
|
|
393 """
|
