annotate call_amr_mutations.py @ 4:2774bdf30122 draft

Uploaded
author greg
date Tue, 21 Mar 2023 19:32:25 +0000
parents 30b17a17c1aa
children bafbed02fdd2
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3
30b17a17c1aa Uploaded
greg
parents:
diff changeset
1 #!/usr/bin/env python
30b17a17c1aa Uploaded
greg
parents:
diff changeset
2
30b17a17c1aa Uploaded
greg
parents:
diff changeset
3 # NOTE: This tool provides the functionality of the PIMA filter_varsacn()
30b17a17c1aa Uploaded
greg
parents:
diff changeset
4 # function here https://github.com/appliedbinf/pima_md/blob/main/pima.py#L3012
30b17a17c1aa Uploaded
greg
parents:
diff changeset
5
30b17a17c1aa Uploaded
greg
parents:
diff changeset
6 import argparse
30b17a17c1aa Uploaded
greg
parents:
diff changeset
7 import re
30b17a17c1aa Uploaded
greg
parents:
diff changeset
8 import subprocess
30b17a17c1aa Uploaded
greg
parents:
diff changeset
9 import sys
30b17a17c1aa Uploaded
greg
parents:
diff changeset
10
30b17a17c1aa Uploaded
greg
parents:
diff changeset
11
30b17a17c1aa Uploaded
greg
parents:
diff changeset
12 def run_command(self, command):
30b17a17c1aa Uploaded
greg
parents:
diff changeset
13 try:
30b17a17c1aa Uploaded
greg
parents:
diff changeset
14 return re.split('\\n', subprocess.check_output(command, shell=True).decode('utf-8'))
30b17a17c1aa Uploaded
greg
parents:
diff changeset
15 except Exception:
30b17a17c1aa Uploaded
greg
parents:
diff changeset
16 message = 'Command %s failed: exiting...' % command
30b17a17c1aa Uploaded
greg
parents:
diff changeset
17 sys.exit(message)
30b17a17c1aa Uploaded
greg
parents:
diff changeset
18
30b17a17c1aa Uploaded
greg
parents:
diff changeset
19
30b17a17c1aa Uploaded
greg
parents:
diff changeset
20 def filter_varscan(varscan_raw, output):
30b17a17c1aa Uploaded
greg
parents:
diff changeset
21 cmd = ' '.join(['cat', varscan_raw,
30b17a17c1aa Uploaded
greg
parents:
diff changeset
22 '| awk \'(NR > 1 && $9 == 2 && $5 + $6 >= 15)',
30b17a17c1aa Uploaded
greg
parents:
diff changeset
23 '{OFS = "\\t";f = $6 / ($5 + $6); gsub(/.*\\//, "", $4);s = $4;gsub(/[+\\-]/, "", s);$7 = sprintf("%.2f%%", f * 100);'
30b17a17c1aa Uploaded
greg
parents:
diff changeset
24 'min = 1 / log(length(s) + 2) / log(10) + 2/10;if(f > min){print}}\'',
30b17a17c1aa Uploaded
greg
parents:
diff changeset
25 '1>' + output])
30b17a17c1aa Uploaded
greg
parents:
diff changeset
26 output = run_command(cmd)
30b17a17c1aa Uploaded
greg
parents:
diff changeset
27
30b17a17c1aa Uploaded
greg
parents:
diff changeset
28
30b17a17c1aa Uploaded
greg
parents:
diff changeset
29 if __name__ == '__main__':
30b17a17c1aa Uploaded
greg
parents:
diff changeset
30 parser = argparse.ArgumentParser()
30b17a17c1aa Uploaded
greg
parents:
diff changeset
31
30b17a17c1aa Uploaded
greg
parents:
diff changeset
32 parser.add_argument('--varscan_raw', action='store', dest='varscan_raw', help='Raw varscan mpileup VCF file')
30b17a17c1aa Uploaded
greg
parents:
diff changeset
33 parser.add_argument('--output', action='store', dest='output', help='Output filtered VCF file')
30b17a17c1aa Uploaded
greg
parents:
diff changeset
34
30b17a17c1aa Uploaded
greg
parents:
diff changeset
35 args = parser.parse_args()
30b17a17c1aa Uploaded
greg
parents:
diff changeset
36
30b17a17c1aa Uploaded
greg
parents:
diff changeset
37 filter_varscan(args.varscan_raw, args.output)