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