python/etc/separateSTARbam.py

   1 #!/usr/bin/env python3
   2
   3 import sys
   4 import os
   5 import pysam
   6
   7 def main():
   8     if len(sys.argv) < 3 :
   9         print('Usage:',sys.argv[0],'<bam file> <output prefix>',file=sys.stderr,flush=True);
  10         exit(0);
  11     bamFile = sys.argv[1]
  12     outPrefix = sys.argv[2]
  13     if os.path.dirname(outPrefix): os.makedirs(os.path.dirname(outPrefix), exist_ok=True)
  14     uniqBamName = ''.join((outPrefix,'.uniq.bam'))
  15     multiBamName = ''.join((outPrefix,'.multi.bam'))
  16     nfoName = ''.join((outPrefix,'.info'))
  17     readCnt = {'Uniq':0, 'Multi':0, 'Secondary':0, 'Total':0}
  18
  19     with pysam.AlignmentFile(bamFile, "rb", threads=4) as samIn:
  20         uniqBam = pysam.AlignmentFile(uniqBamName, "wb", template=samIn,threads=4)
  21         multiBam = pysam.AlignmentFile(multiBamName, "wb", template=samIn,threads=4)
  22         for aln in samIn:
  23             readCnt['Total'] += 1
  24             if aln.mapq == 255:
  25                 readCnt['Uniq'] += 1
  26                 uniqBam.write(aln)
  27             elif not aln.is_secondary:
  28                 readCnt['Multi'] += 1
  29                 multiBam.write(aln)
  30             else:
  31                 readCnt['Secondary'] += 1
  32         uniqBam.close()
  33         multiBam.close()
  34
  35     with open(nfoName,'w') as nfo:
  36         nfo.write(str(readCnt))
  37
  38 if __name__ == "__main__":
  39     main()  # ./separateSTARbam.py Aligned.sortedByCoord.out.bam Separated &