#!/usr/bin/env python3

import re, os, sys
from argparse import ArgumentParser, RawDescriptionHelpFormatter
from collections import Counter

import pysam

#############
## Read Me ##
#############
# Complete redesign of the discard step by Tobias Neumann.
#
# 22/06/2020 Created
#
# This just scans through the BAM, and reports the number of unique IDs for mapped reads. Unpapped reads will be ignored.
#
# Last revised: 26/jun/2020   by: kimon.froussios@imp.ac.at

#################
## Parameters ##
################

usage = "Count unique read IDs in a BAM. Partially overlapping coordinates are included."
version = "0.3.0"
# Main Parsers
parser = ArgumentParser(description=usage, formatter_class=RawDescriptionHelpFormatter)
parser.add_argument("-b", "--bam", type=str, nargs='+', required=True, dest="allSamFile", help="BAM file with all the mapped reads")
parser.add_argument("-c", "--chr", type=str, required=False, dest="chromosome", help="chromosome")
parser.add_argument("-s", "--start", type=int, required=False, dest="coordstart", help="from coordinate, 0-based, inclusive")
parser.add_argument("-e", "--end", type=int, required=False, dest="coordend", help="to coordinate, 0-based ,not inclusive")
args = parser.parse_args()

##########
# Collect read names
##########

for f in args.allSamFile:
    readID = Counter()

    inSam = pysam.AlignmentFile(f, "rb")
    # Multiple alignments of a read are reported individually.
    for read in inSam:
        if not read.is_unmapped:
            flag = True   # keeping track of if all conditions are individually met, allowing for any combination of conditions to be specified
            if args.chromosome is not None and read.reference_name != args.chromosome:  # reference sequence outside specification
                flag = False
            if args.coordstart is not None and read.reference_end < args.coordstart:    # read mapped before specified location
                flag = False
            if args.coordend is not None and read.reference_start > args.coordend:      # read mapped after specified location
                flag = False        
            
            if flag:
                readID.update([read.query_name])
    inSam.close()

    i = 0
    for k in readID.items():
        if k[1] > 1:
            i += 1 

    sys.stdout.write(f + '\tReads:\t' + str(len(readID)) + '\tof which multimappers:\t' + str(i) + '\n')

@@ -395,16 +395,190 @@
 ignored.
+ %0A# Multimapping is always counted through the whole file, but the reads included in the report can be restricted to those with at least one alignment overlapping the region.
 %0A#%0A# Las
@@ -589,17 +589,17 @@
 vised: 2
-6
+9
 /jun/202
@@ -1451,16 +1451,113 @@
 ######%0A%0A
+# Count multimappers in the whole file, then fish out the ones that overlap the relevant region.%0A
 for f in
@@ -1596,16 +1596,37 @@
 ounter()
+%0A    relevant = set()
 %0A%0A    in
@@ -1779,16 +1779,62 @@
 mapped:%0A
+            readID.update(%5Bread.query_name%5D)%0A%0A
         
@@ -2454,21 +2454,8 @@
    %0A
-            %0A
     
@@ -2485,37 +2485,35 @@
           re
-adID.update(%5B
+levant.add(
 read.query_n
@@ -2515,17 +2515,16 @@
 ery_name
-%5D
 )%0A    in
@@ -2565,20 +2565,14 @@
 n re
-adID.items()
+levant
 :%0A  
@@ -2584,11 +2584,16 @@
  if 
-k%5B1
+readID%5Bk
 %5D %3E 
@@ -2649,16 +2649,26 @@
 '%5CtReads
+ in region
 :%5Ct' + s
@@ -2676,20 +2676,22 @@
 r(len(re
-adID
+levant
 )) + '%5Ct

Update CountMappedReadsBAM.py