Commits

Davide Cittaro committed f4cdaa5

Full extract SM, to test

Comments (0)

Files changed (1)

 
 samples = args
 tokeep = {} #a dictionary files -> list of read groups 
+rg2sm = {}
+sample2fo = {}
 
 bamh = pysam.Samfile(bam_file, 'rb')
 for rg in bamh.header['RG']:
     # create one file for each
     if extract_by_rg:
       tokeep[rg['ID']] = [rg['ID']]
+      rg2sm[rg['ID']] = rg['ID']
     else:
       try:
         tokeep[rg['SM']].append(rg['ID'])
+        rg2sm[rg['ID']] = rg['SM']
       except KeyError:  
         tokeep[rg['SM']] = [rg['ID']]
+        rg2sm[rg['ID']] = rg['SM']
   else:
     if extract_by_rg:
       if rg['ID'] in samples:
         tokeep[rg['ID']] = [rg['ID']]
+        rg2sm[rg['ID']] = rg['ID']
     else:
       if rg['SM'] in samples:
         try:
           tokeep[rg['SM']].append(rg['ID'])
+          rg2sm[rg['ID']] = rg['SM']
         except KeyError:
           tokeep[rg['SM']] = [rg['ID']]
+          rg2sm[rg['ID']] = rg['SM']
 
+if not tokeep:
+  sys.stderr.write("No samples matched, exiting\n")
+  sys.exit(1)
 
-print tokeep          
+for sm in tokeep.keys():
+  fname = extract_prefix + sn + '.bam'
+  header = bamh.header
+  if extract_by_rg:
+    rh = [x for x in header['RG'] if x['RG'] == sm]
+  else:
+    rh = [x for x in header['RG'] if x['SM'] == sm]
+  header['RG'] = rh
+  sample2fo[sm] = pysam.Samfile(fname, 'wb', header=header)
+  
+  
+
+for r in bamh.fetch():
+  rg = [x[1] for x in r.tags if x[0] == 'RG'][0]
+  sm = rg2sm[rg]
+  sample2fo[sm].write(r)
+  
+for sm in sample2fo.keys():
+  sample2fo[sm].close()            
+
+
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.