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