-
Notifications
You must be signed in to change notification settings - Fork 15
/
bam_coverage_join.py
executable file
·49 lines (37 loc) · 1.33 KB
/
bam_coverage_join.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#!/usr/bin/python
import sys
from subprocess import call
from os import listdir
from os.path import isfile, join
print "bam_coverage_join.py FastaReference ListOfBams [Maximum Depth (Default: 8000)]"
try:
ref = sys.argv[1]
except:
ref = raw_input("Introduce FASTA reference: ")
try:
li = sys.argv[2]
except:
li = raw_input("Introduce list of BAM files: ")
try:
depth = sys.argv[3]
except:
depth = 8000
bams = open(li).readlines()
to_join = []
onlyfiles = [f for f in listdir(".") if isfile(join(".",f))]
for bam in bams:
bam = bam[:-1]
if bam+".var2" not in onlyfiles:
if bam+".var" not in onlyfiles:
print "Not found file %s.var, generating..." % bam
call("pysamstats --max-depth=%s -f %s --type variation %s > %s.var" % (depth, ref,bam,bam), shell=True)
call("""awk {'print $1"popeye"$2"\t"$4'} %s.var > %s.var2""" % (bam,bam), shell=True)
print "Not found file %s.var2, generating..." % bam
else:
print "Found file %s.var2, skipping." % bam
to_join.append(bam+".var2")
print "Joining lists..."
call("join_multiple_lists.py %s " % " ".join(to_join), shell=True)
call("""sed 's/popeye/\t/g' toico.txt > toico2.txt""", shell=True)
call("(head -n 1 toico2.txt && tail -n +2 toico2.txt | sort -k 1,1 -k 2n) > toico3.txt ", shell=True)
print "We're done!!!\n"