1
2
3 import os, sys
4 import argparse
5 import subprocess
6
7 __doc__="""
8 Extract mapped or unmapped lines of SAM file
9 with or without header
10 and can sort result by query or subject
11
12 @requires: U{SamTools<http://sourceforge.net/projects/samtools/files/>}
13 """
14
16
17 parser = argparse.ArgumentParser(description='Extract mapped or unmapped lines of SAM file with or without header and can sort result by query or subject')
18
19 parser.add_argument('-sam', action="store", dest='samFile',
20 type=str, help='SAM file')
21
22 parser.add_argument('-head', action="store_true", dest='header',
23 help='keep header', default=False)
24
25 parser.add_argument('-map', action="store_true", dest='map',
26 help='keep only mapped lines', default=False)
27
28 parser.add_argument('-unmap', action="store_true", dest='unmap',
29 help='keep only unmapped lines', default=False)
30
31 parser.add_argument('-Q', action="store_true", dest='byQuery',
32 help='sort by query', default=False)
33
34 parser.add_argument('-S', action="store_true", dest='bySubject',
35 help='sort by subject', default=False)
36
37 return parser
38
39
41
42 parser=get_parser()
43
44 if len(sys.argv)==1:
45 parser.print_help()
46 sys.exit(1)
47
48 Arguments=parser.parse_args()
49
50 prefix = Arguments.samFile.split('.')
51 prefix = prefix[len(prefix)-2]
52
53 if(Arguments.map):
54 prefix = prefix + "_mapped.sam"
55 commandeLine = "samtools view -S -h -F 4 " + Arguments.samFile + " -o " + prefix
56 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
57 result=process.communicate()
58
59 elif(Arguments.unmap):
60 prefix = prefix + "_unmapped.sam"
61 commandeLine = "samtools view -S -h -f 4 " + Arguments.samFile + " -o " + prefix
62 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
63 result=process.communicate()
64
65 if(Arguments.byQuery):
66 if(Arguments.map or Arguments.unmap):
67 file = prefix
68 else:
69 file = Arguments.samFile
70 prefix = prefix.split('.')[0] + "_sortbyquery.sam"
71 if(Arguments.header):
72 commandeLine = "grep '^@' " + file + " > " + prefix
73 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
74 result=process.communicate()
75 commandeLine = "grep -v '^@' " + file + "| sort -k1 >> " + prefix
76 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
77 result=process.communicate()
78 if(Arguments.map):
79 commandeLine = "rm " + Arguments.samFile.split('.')[0] + "_mapped.sam"
80 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
81 result=process.communicate()
82 if(Arguments.unmap):
83 commandeLine = "rm " + Arguments.samFile.split('.')[0] + "_unmapped.sam"
84 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
85 result=process.communicate()
86
87 elif(Arguments.bySubject):
88 if(Arguments.map or Arguments.unmap):
89 file = prefix
90 else:
91 file = Arguments.samFile
92 prefix = prefix.split('.')[0] + "_sortbysubject.sam"
93 if(Arguments.header):
94 commandeLine = "grep '^@' " + file + " > " + prefix
95 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
96 result=process.communicate()
97 commandeLine = "grep -v '^@' " + file + "| sort -k3 >> " + prefix
98 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
99 result=process.communicate()
100 if(Arguments.map):
101 commandeLine = "rm " + Arguments.samFile.split('.')[0] + "_mapped.sam"
102 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
103 result=process.communicate()
104 if(Arguments.unmap):
105 commandeLine = "rm " + Arguments.samFile.split('.')[0] + "_unmapped.sam"
106 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
107 result=process.communicate()
108
109 elif(not Arguments.header and not Arguments.byQuery and not Arguments.bySubject):
110 if(Arguments.map or Arguments.unmap):
111 file = prefix
112 else:
113 file = Arguments.samFile
114 commandeLine = "grep -v '^@' " + file + " > " + file + ".tmp"
115 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
116 result=process.communicate()
117 commandeLine = "mv " + file + ".tmp " + file
118 process = subprocess.Popen(commandeLine,shell=True,stdout=subprocess.PIPE, stderr=subprocess.PIPE)
119 result=process.communicate()
120
121 if __name__ == "__main__":
122 main()
123