s3_categorizePairs.py

#!/usr/bin/env python

## mHiC
## Author: Ye Zheng
## Contact: yezheng@stat.wisc.edu

'''
Script to identify Restriction Enzyme fragment where the each read end comes from.
Categorize each read pairs: SC, DE, Religation, Dump, valid read pairs.
Summarize categories.
Only valid read pairs with fragment assignment are saved in output.
Dec, 2016
'''

import sys
import os
import re
import pysam
import argparse
from bx.intervals.intersection import Intersecter, Interval

def get_args():
    '''Get arguments'''
    parser = argparse.ArgumentParser(description = '------------Usage Start------------',
                                     epilog = '------------Usage End------------')
    parser.add_argument('-f', '--fragment', help = 'Restriction Enzyme(RE) fragment file.', default = None)
    parser.add_argument('-r', '--read', help = 'Paired-ended read alignment file.', default = None)
    parser.add_argument('-o', '--outdir', help = 'Output directory to save category results.', default = None)
    parser.add_argument('-l', '--lower', help = 'Lower bound of the read pair distances summation from read end alignment position to its assigned RE fragment cutting site. Default value is None indicating no restriction.', default = None)
    parser.add_argument('-u', '--upper', help = 'Upper bound of the read pair distances summation from read end alignment position to its assigned RE fragment cutting site. Default value is None indicating no restriction.', default = None)
    parser.add_argument('-d', '--distance', help = 'Minimum distance between intrachromosomal read pair alignments positions. Default value is None indicating no restriction.', default = None)
    parser.add_argument('-m', '--method', help = "Binning methods: by fixed number of window size, i.e. window, or fixed number of RE fragment, i.e. fragment.", default = "fragment")
    parser.add_argument('-b', '--bin', help = "Binning resolution. If binning by fixed number of window size, it can be 10000 or any other suitable bin size. If binning by fixed number of RE fragment, it can be 10 representing 10RE fragments or other suitable bin size.", default = None)
    parser.add_argument('-s', '--summary', help = '(Optional) Summary of categories. Default is true.', default = True)
    parser.add_argument('-sf', '--summaryFile', help = '(Optional) Summary file path. Default is alignedFileName.readCategorySummary.', default = None)
    parser.add_argument('-v', '--verbose', help = '(Optional) Verbose. Default is true.', default = True)

    args = parser.parse_args()
    if args.fragment is None or args.read is None or args.outdir is None:
        parser.print_help()
        sys.exit()
    if not os.path.exists(args.fragment):
        print("Resctriction Enzyme fragment file does not exist!")
        sys.exit()
    if not os.path.exists(args.read):
        print("Paired-ended read alignmeng file does not exist!")
        sys.exit()
    if not os.path.exists(args.outdir):
        os.makedirs(args.outdir)
    methodList = ['fragment', 'window', 'frag', 'win', 'f', 'w']
    if not any(args.method in x for x in methodList):
        print("Binning method should either be fragment or window.")
        sys.exit()
    if args.bin is None:
        print("Please provide number representing binning resolution.")
        sys.exit()
    else:
        args.bin = int(args.bin)
    return args


def get_read_strand(read):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''


    """
    Conversion of read position to naive strand representation
    Parameters
    ----------
    read : list
        list of aligned reads
    """
    strand = "+"
    if read.is_reverse:
        strand = "-"
    return strand


def isIntraChrom(read1, read2):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    
    """
    Return true is the reads pair is intrachromosomal
    
    read1 : [AlignedRead]
    read2 : [AlignedRead]
    """
    if read1.tid == read2.tid:
        return True
    else:
        return False


def get_read_pos(read):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''

    """
    Return the read position (zero-based) used for the intersection with
    the restriction fragment
    The 5' end is not a good choice for the reverse reads (which contain part
    of the restriction site, and thus overlap the next restriction fragment)
    Using the left-most position (5' for forward, 3' for reverse) or the
    middle of the read should work but the middle of the reads might be more
    safe
    Parameters
    -----------
    read : list
        list of aligned reads
    """
    if read.alen is not None:
        pos = read.pos + read.alen/2
    elif read.cigarstring is not None:
        alignedLen = re.findall(r'(\d+)[M, D]', read.cigarstring)
        len = sum([int(i) for i in alignedLen])
        pos = read.pos + len/2
    else:
        print("Warning: neither read aligned length nor cigar string!")
        sys.exit()

    return pos


def get_cis_dist(read1, read2):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    """
    Calculte the contact distance between two intrachromosomal reads
    read1 : [AlignedRead]
    read2 : [AlignedRead]
    """
    # Get oriented reads
    ##r1, r2 = get_ordered_reads(read1, read2)
    dist = None
    if not r1.is_unmapped and not r2.is_unmapped:         
        ## Contact distances can be calculated for intrachromosomal reads only
        if isIntraChrom(read1, read2):
            r1pos = get_read_pos(read1)
            r2pos = get_read_pos(read2)
            dist = abs(r1pos - r2pos)
            return dist
        
        
def get_read_start(read):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    
    """
    Return the 5' end of the read
    """
    if read.is_reverse:
        if read.alen is not None:
            pos = int(read.pos + read.alen)
        elif read.cigarstring is not None:
            alignedLen = re.findall(r'(\d+)[M, D]', read.cigarstring)
            len = sum([int(i) for i in alignedLen])
            pos = int(read.pos + len)
        else:
            print("Warning: neither read aligned length nor cigar string!")
            sys.exit()
            
    else:
        pos = read.pos
    return pos

def get_ordered_reads(read1, read2):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    """
    Reorient reads
    The sequencing is usually not oriented. Reorient the reads so that r1 is
    always before r2.
   
    read1 = [AlignedRead]
    read2 = [AlignedRead]
    """
    if read1.tid == read2.tid:
        if get_read_pos(read1) < get_read_pos(read2):
            r1 = read1
            r2 = read2
        else:
            r1 = read2
            r2 = read1
    else:
        if read1.tid < read2.tid:
            r1 = read1
            r2 = read2
        else:
            r1 = read2
            r2 = read1
                
    return r1, r2


def get_overlapping_restriction_fragment(resFrag, chrom, read):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    """
    Intersect a given read with the set of restriction fragments
    ##
    resFrag = the restriction fragments [hash]
    chrom = the chromosome to look at [character]
    read = the read to intersect [AlignedRead]
    """
    # Get read position (middle or 5' end)
    pos = int(get_read_pos(read))
    
    if chrom in resFrag:
        # Overlap with the position of the read (zero-based)
        resfrag = resFrag[chrom].find(pos, pos+1)
        if len(resfrag) > 1:
            print("Warning : " +  len(resfrag) + " restriction fragments found for " + read.qname + "- skipped")
            return None
        elif len(resfrag) == 0:
            print("Warning - no restriction fragments for " + read.qname + " at " + chrom + ":" + str(pos))
            return None
        else:
            return resfrag[0]
    else:
        print("Warning - no restriction fragments for " + read.qname + " at " + chrom + ":" + str(pos))
        return None

def load_restriction_fragment(in_file, minfragsize=None, maxfragsize=None, verbose=False):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    '''
    Function load_restriction_fragment cite the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara 
    '''
    
    """
    Read a BED file and store the intervals in a tree
    Intervals are zero-based objects. The output object is a hash table with
    one search tree per chromosome
    in_file = input file [character]
    verbose = verbose mode [logical]
    """
    resFrag = {}
    if verbose:
        print("## Loading Restriction File Intervals '" + in_file + "'...")

    bed_handle = open(in_file)
    nline = 0
    for line in bed_handle:
        nline +=1
        bedtab = line.split("\t")
        try:
            chromosome, start, end, name = bedtab[:4]
        except ValueError:
            print("Warning : wrong input format in line" + nline + ". Not a BED file !?")
            continue

        # BED files are zero-based as Intervals objects
        start = int(start)  # + 1
        end = int(end)
        midPoint = (start + end)/2
        fragl = abs(end - start)
        name = name.strip()

        ## Discard fragments outside the size range
        if minfragsize != None and int(fragl) < int(minfragsize):
            print("Warning : fragment "+ name + " [" +  fragl + "] outside of range. Discarded")
            continue
        if maxfragsize != None and int(fragl) > int(maxfragsize):
            print("Warning : fragment " + name + " [" + fragl + "] outside of range. Discarded")
            continue
       
        if chromosome in resFrag:
            tree = resFrag[chromosome]
            tree.add_interval(Interval(start, end, value={'name': name, 'midPoint': midPoint}))
        else:
            tree = Intersecter()
            tree.add_interval(Interval(start, end, value={'name': name, 'midPoint': midPoint}))
            resFrag[chromosome] = tree
    
    bed_handle.close()
    return resFrag



def get_overlapping_resFragBin(binTree, chrom, read):
    """
    Intersect a given read with the set of bin
    ##
    binTree = the bin tree [hash]
    chrom = the chromosome to look at [character]
    read = the read to intersect [AlignedRead]
    """
    # Get read position (middle or 5' end)
    pos = int(get_read_pos(read))
    
    if chrom in binTree.keys():
        # Overlap with the position of the read (zero-based)
        resFragBin = binTree[chrom].find(pos, pos+1)
        if len(resFragBin) > 1:
            print("Warning : " +  str(len(resFragBin)) + " bins found for " + read.qname + "- skipped")
            print(pos)
            print(resFragBin)
            return None
        elif len(resFragBin) == 0:
            print("Warning - no bin found for " + read.qname + " at " + chrom + ":" + str(pos))
            return None
        else:
            return resFragBin[0]
    else:
        print("Warning - no bin found for " + read.qname + " at " + chrom + ":" + str(pos))
        return None

def build_resFragBin_tree(in_file, resolution = None, verbose = False):
    """
    build restriction fragments bins tree based on the user defined resolution, e.g. 10fragment
    """
    binTree = {}
    if resolution is None:
        print("Please provide suitable resolution to bin contact! For example resolution = 10 for 10 RE fragment as one bin.")
        sys.exit()
    else:
        resolution = int(resolution)
        
    if verbose:
        print("## Building RE fragment bins tree from ordered'" + in_file + "'...")
    bed_handle = open(in_file)
    nline = 0
    flag = 0
    for line in bed_handle:
        nline += 1

        bedtab = line.split("\t")
        try:
            chromosome, start, end = bedtab[:3]
            start = int(start)
            end = int(end)
        except ValueError:
            print("Warning : wrong input format in line" + nline + ". Not a BED file !?")
            continue
            
        # First line
        if flag == 0:
            #Bin variables store current bin information
            startBin = start
            chromosomeBin = chromosome
            flag = 1
        elif chromosomePre != chromosome: #Start of another chromosome
            nline = 1
            if startBin != endPre:
                endBin = endPre
                midBin = round((startBin + endBin)/2)
                   
                if chromosomeBin in binTree.keys():
                    tree = binTree[chromosomeBin]
                    tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
                else:
                    tree = Intersecter()
                    tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
                    binTree[chromosomeBin] = tree
            startBin = start
            chromosomeBin = chromosome
        elif nline % resolution == 0:
            endBin = end
            midBin = round((startBin + endBin)/2)
                   
            if chromosomeBin in binTree.keys():
                tree = binTree[chromosomeBin]
                tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
            else:
                tree = Intersecter()
                tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
                binTree[chromosomeBin] = tree

            startBin = end
            chromosomeBin = chromosome

        # Update the Pre variable set
        startPre = start
        endPre = end
        chromosomePre = chromosome
        
    # for the last bin
    if nline % resolution != 0 and startBin != start and chromosomeBin == chromosome:
        endBin = end
        midBin = round((startBin + endBin)/2)
                   
        if chromosomeBin in binTree.keys():
            tree = binTree[chromosomeBin]
            tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
        else:
            tree = Intersecter()
            tree.add_interval(Interval(startBin, endBin, value={'midPoint': midBin}))
            binTree[chromosomeBin] = tree

    bed_handle.close()
    return binTree




            

def are_contiguous_fragments(frag1, frag2, chr1, chr2):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    '''
    Compare fragment positions to check if they are contiguous
    '''
    ret = False
    if chr1 == chr2:
        if int(frag1.start) < int(frag2.start):
            d = int(frag2.start) - int(frag1.end)
        else:
            d = int(frag1.start) - int(frag2.end)
        
        if d == 0:
            ret = True
            
            return ret

def is_religation(read1, read2, frag1, frag2):
    '''
    This function cites the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara
    '''
    """
    Reads are expected to map adjacent fragments
    Check the orientation of reads -><-
    """
    ret=False
    if are_contiguous_fragments(frag1, frag2, read1.tid, read2.tid):
        ret=True
    return ret


def is_self_circle(read1, read2):
    '''
    Function cite the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara 
    '''
    """
    Both reads are expected to be on the same restriction fragments
    Check the orientation of reads <-->
    read1 : [AlignedRead]
    read2 : [AlignedRead]
    """
    ret = False
    # Get oriented reads
    r1, r2 = get_ordered_reads(read1, read2)
    # 1<- ->2 or 2<- ->1
    if get_read_strand(r1) == "-" and get_read_strand(r2) == "+":
        ret = True
    return ret


def is_dangling_end(read1, read2):
    '''
    Function cite the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara 
    '''
    """
    Both reads are expected to be on the same restriction fragments
    Check the orientation of reads -><-
    read1 : [AlignedRead]
    read2 : [AlignedRead]
    """
    ret = False
    # Get oriented reads
    r1, r2 = get_ordered_reads(read1, read2)
    # 1-> <-2 or 2-> <-1
    if get_read_strand(r1) == "+" and get_read_strand(r2) == "-":
        ret = True
    return ret


def get_interaction_type(read1, read1_chrom, resfrag1, read2, read2_chrom, resfrag2, verbose):
    '''
    Function cite the same function from Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara 
    '''
    """
    Returns the interaction type
    For a given reads pair and their related restriction fragment, classify
    the 3C products as :
    - Singleton Interaction (SI)
    - Valid Interaction (VI)
    - Self circle (SC)
    - Dangling end (DE)
    - Religation (RE)
    - Unknown (DUMP)

    ##
    read1 = the R1 read of the pair [AlignedRead]
    read1_chrom = the chromosome of R1 read [character]
    resfrag1 = restrictin fragment overlapping the R1 read [interval]
    read2 = the R2 read of the pair [AlignedRead]
    read2_chrom = the chromosome of R2 read [character]
    resfrag2 = restrictin fragment overlapping the R2 read [interval]
    verbose = verbose mode [logical]
    """

    # If returned InteractionType=None -> Same restriction fragment
    # and same strand = Dump
    interactionType = None
 
    if not r1.is_unmapped and not r2.is_unmapped and resfrag1 is not None and resfrag2 is not None:
        # same restriction fragment
        if resfrag1 == resfrag2:
            # Self_circle <- ->
            if is_self_circle(read1, read2):
                interactionType = "SC"
            # Dangling_end -> <-
            elif is_dangling_end(read1, read2):
                interactionType = "DE"
        elif is_religation(read1, read2, resfrag1, resfrag2):
            interactionType = "RE"
        else:
            interactionType = "VI"
    elif r1.is_unmapped or r2.is_unmapped:
        interactionType = "SI"

    return interactionType

def get_PE_fragment_size(read1, read2, resFrag1, resFrag2, interactionType):
    '''
    Function modified from function of Hi-C Pro (https://github.com/nservant/HiC-Pro/blob/master/scripts/mapped_2hic_fragments.py) by Nicolas Servant, Eric Viara 
    '''
    """
    Calculte the size of the DNA fragment library - (Ye: distances from read end alignment to its fragment cutting site.)
    read1 : [AlignedRead]
    read2 : [AlignedRead]
    resfrag1 = restrictin fragment overlapping the R1 read [interval]
    resfrag2 = restrictin fragment overlapping the R2 read [interval]
    interactionType : Type of interaction from get_interaction_type() [str]
    """

    fragmentsize = None

    # Get oriented reads
    r1, r2 = get_ordered_reads(read1, read2)
    if not r1.is_unmapped and not r2.is_unmapped:
        if r1 == read2:
            rfrag1 = resFrag2
            rfrag2 = resFrag1
        else:
            rfrag1 = resFrag1
            rfrag2 = resFrag2

        ## In this case use the read 5' end !
        r1pos = get_read_start(r1)
        r2pos = get_read_start(r2)

        if interactionType == "DE" or interactionType == "RE":
            fragmentsize = r2pos - r1pos
        elif interactionType == "SC":
            fragmentsize = (r1pos - rfrag1.start) + (rfrag2.end - r2pos)
        elif interactionType == "VI":
            dr1 = min(abs(rfrag1.end - r1pos), abs(r1pos - rfrag1.start))
            dr2 = min(abs(rfrag2.end - r2pos), abs(r2pos - rfrag2.start))
            # if get_read_strand(r1) == "+":
            #     dr1 = rfrag1.end - r1pos
            # else:
            #     dr1 = r1pos - rfrag1.start
            # if get_read_strand(r2) == "+":
            #     dr2 = rfrag2.end - r2pos
            # else:
            #     dr2 = r2pos - rfrag2.start
            fragmentsize = dr2 + dr1

    return fragmentsize

def get_valid_pair(r1, r2, r1_chrom, r2_chrom, r1_frag, r2_frag, r1_bin, r2_bin, cisDist):
    # print valid read pairs out
    or1, or2 = get_ordered_reads(r1, r2)

    if or1 == r1 and or2 == r2:
        or1_chrom = r1_chrom
        or2_chrom = r2_chrom
        or1_fragname = r1_frag.value['name'] if r1_frag is not None else "NA"
        or2_fragname = r2_frag.value['name'] if r2_frag is not None else "NA"
        # or1_fragPos = r1_frag.value['midPoint'] if r1_frag is not None else 0
        # or2_fragPos = r2_frag.value['midPoint'] if r2_frag is not None else 0
        or1_bin = r1_bin
        or2_bin = r2_bin
    elif or1 == r2 and or2 == r1:
        or1_chrom = r2_chrom
        or2_chrom = r1_chrom
        or1_fragname = r2_frag.value['name'] if r2_frag is not None else "NA"
        or2_fragname = r1_frag.value['name'] if r1_frag is not None else "NA"
        # or1_fragPos = r2_frag.value['midPoint'] if r2_frag is not None else 0
        # or2_fragPos = r1_frag.value['midPoint'] if r1_frag is not None else 0
        or1_bin = r2_bin
        or2_bin = r1_bin

    validInfo = "\t".join([or1.qname, or1_chrom, str(get_read_pos(or1)), get_read_strand(or1), or1_fragname, str(or1_bin),\
                           or2_chrom, str(get_read_pos(or2)), get_read_strand(or2), or2_fragname, str(or2_bin),\
                           "InterChrom" if cisDist is None else str(cisDist)])
    return validInfo

def get_multi_read(readFile, readName, tag):
    read = pysam.AlignedSegment()

    read.qname = readName
    readInfo = tag.split(',')
    read.tid = readFile.get_tid(readInfo[0])
    read.is_reverse = True if readInfo[1].startswith('-') else False
    read.pos = int(readInfo[1].strip('-')) - 1
    read.cigarstring = readInfo[2].strip()
    return read
    


if __name__ == "__main__":
    args = get_args()

    if args.verbose:
        print("RE fragment = ", args.fragment)
        print("Paired-ended read alignment file = ", args.read)
        print("Output directory = ", args.outdir)
        print("RE cutting site distance lower bound = ", args.lower)
        print("RE cutting site distance upper bound = ", args.upper)
        print("Minimum intrachromosomal contact distance = ", args.distance)
        print("Binning method = ", args.method)
        print("Binning resolution = ", args.bin)
        print("Summary = ", args.summary)
        print("SummaryFile = ", args.summaryFile)
        print("Verbose = ", args.verbose)

    # count for uni-mapping read pairs
    sc_count = 0
    de_count = 0
    re_count = 0
    ud_count = 0 #uni dump
    si_count = 0
    uv_count = 0
    uni_count = 0
    # count for multi-mapping read pairs/multi-mapping unique valid read pairs/multi-mapping multiple valid read pairs
    multi_count = 0
    md_count = 0 #multi dump
    multi_uv_count = 0
    multi_mv_count = 0
    multi_mv_uf_count = 0
    multi_mv_ub_count = 0
    
    fileName = re.sub("\.bam$|\.sam$", "", os.path.basename(args.read))
    outValid = open(args.outdir + "/" + fileName + ".validPairs", "w")
    outContactNum = open(args.outdir + "/" + fileName + ".uniMulti", "w")
    
    fragTree = load_restriction_fragment(in_file = args.fragment, verbose = args.verbose)
    print(fragTree.keys())

    fragMethodList = ['fragment', 'frag', 'f']
    if any(args.method in x for x in fragMethodList):
        binTree = build_resFragBin_tree(in_file = args.fragment, resolution = args.bin, verbose = args.verbose)
    else:
        binTree = None

    fragMethodList = ['window', 'win', 'w', 'windows']
    if any(args.method in x for x in fragMethodList):
        binWindow = int(args.bin)
    else:
        binWindow = None
        
    if args.verbose:
        print("------------Begin reading paired-ended file------------")

    if args.read.endswith(".bam"):
        readFile = pysam.Samfile(args.read, "rb")
    elif args.read.endswith(".sam"):
        readFile = pysam.Samfile(args.read, "r")
    else:
        print("Input paired-ended alignment file should either be in bam format or sam format.")
        sys.exit()

    if args.verbose:
        print("------------Categorize aligned read pairs------------")

    # new one with all kinds of potential contacts
    for read in readFile.fetch(until_eof = True):
        if read.is_read1:
            if not read.is_unmapped:
                r1 = read
            else:
                r1 = None
        elif read.is_read2:
            if not read.is_unmapped:
                r2 = read
            else:
                r2 = None

            if r1 is not None and r2 is not None:
                m1 = [r1]
                if r1.has_tag('XA'):
                    for xaTag in r1.get_tag('XA').strip(';').split(';'):
                        m1.append(get_multi_read(readFile, r1.qname, xaTag))
                m2 = [r2]
                if r2.has_tag('XA'):
                    for xaTag in r2.get_tag('XA').strip(';').split(';'):
                        m2.append(get_multi_read(readFile, r2.qname, xaTag))
                valid_count = 0
                fragPair = []
                binPair = []
                validInfoList = []
                for end1 in m1:
                    for end2 in m2:
                        
                        end1_chrom = readFile.getrname(end1.tid)
                        end2_chrom = readFile.getrname(end2.tid)
                        end1_frag = get_overlapping_restriction_fragment(fragTree, end1_chrom, end1)
                        if end1_frag is None:
                            interactionType = "DUMP"
                            break

                        end2_frag = get_overlapping_restriction_fragment(fragTree, end2_chrom, end2)
                        if end2_frag is None:
                            interactionType = "DUMP"
                            break
                        # get interaction type, alignmenet position to fragment cutting site distance sum, and intrachromosal contact distance
                        if end1_frag is not None and end2_frag is not None:
                            interactionType = get_interaction_type(end1, end1_chrom, end1_frag, end2, end2_chrom, end2_frag, args.verbose)
                            fragDist = get_PE_fragment_size(end1, end2, end1_frag, end2_frag, interactionType)
                            cisDist = get_cis_dist(end1, end2)
                            # If the alignment position to fragment cutting sites distance sum is out of range allowed
                            if (args.lower is not None and fragDist is not None and fragDist < int(args.lower)) or \
                               (args.upper is not None and fragDist is not None and fragDist > int(args.upper)):
                                interactionType = "DUMP"


                            # If intrachromosal contact distance is too near
                            if interactionType == "VI" and args.distance is not None and cisDist is not None and cisDist < int(args.distance):
                                interactionType = "DUMP"

                            if interactionType == "VI": #if this potential contact is valid - write out
                                valid_count += 1
                                # get the overlapped bin info
                                if binTree != None:
                                    end1_bin = get_overlapping_resFragBin(binTree, end1_chrom, end1)
                                    end2_bin = get_overlapping_resFragBin(binTree, end2_chrom, end2)
                                    mid1 = end1_bin.value["midPoint"]
                                    mid2 = end2_bin.value["midPoint"]
                                elif binWindow != None:
                                    mid1 = int(int(get_read_pos(end1))/binWindow) * binWindow + int(binWindow/2)
                                    mid2 = int(int(get_read_pos(end2))/binWindow) * binWindow + int(binWindow/2)
                                else:
                                    print("No binning method specified!")
                                    sys.exit()

                                if  end1_chrom != end2_chrom or mid1 != mid2:
                                    # Summarize valid pair information and prepare to write out
                                    validInfo = get_valid_pair(end1, end2, end1_chrom, end2_chrom, end1_frag, end2_frag, mid1, mid2, cisDist)
                                    validInfoList.append(validInfo)
                                    
                                    # valid interactions cover how many distinct fragment pairs
                                    fragPair.append([end1.qname, end1_frag.value["name"], end2_frag.value["name"]])
                                    # valid interactions cover how many distinct bins
                                    binPair.append([end1.qname, end1_chrom, mid1, end2_chrom, mid2])
                                else:
                                    interactionType = "DUMP"
                            
                # Identify read pair type: Uni or Multi
                if len(validInfoList) == 1:
                    outValid.write(validInfoList[0] + "\t" + "UNI" + "\n")
                elif len(validInfoList) > 1:
                    for eachInfo in validInfoList:
                        outValid.write(eachInfo + "\t" + "MULTI" + "\n")

                # Identify fragment pair type: Uni or Multi
                if len(fragPair) > 0:
                    fragPairUni = [list(x) for x in set(tuple(x) for x in fragPair)]
                else:
                    fragPairUni = []

                # Identify bin pair type: Uni or Multi
                if len(binPair) > 0:
                    binPairUni = [list(x) for x in set(tuple(x) for x in binPair)]
                else:
                    binPairUni = []

                # Count
                if len(m1) == 1 and len(m2) == 1: #uni-read pair
                    uni_count += 1
                    outContactNum.write('\t'.join([r1.qname, "Uni", "1"]) + "\n")
                    if interactionType == "VI":
                        uv_count += 1
                    elif interactionType == "DE":
                        de_count += 1
                    elif interactionType == "SC":
                        sc_count += 1
                    elif interactionType == "RE":
                        re_count += 1
                    else:
                        interactionType = "DUMP"
                        ud_count += 1
                else: #multi-read pair
                    multi_count +=1
                    if valid_count == 1: #indicating multi-mapping read pair has unique valid contact
                        multi_uv_count += 1
                        outContactNum.write('\t'.join([r1.qname, "MultiToUni", "1"]) + "\n")
                    elif valid_count == 0:
                        md_count += 1
                    else:
                        multi_mv_count += 1
                        outContactNum.write('\t'.join([r1.qname, "Multi", str(valid_count)]) + "\n")
                        if len(fragPairUni) == 1:
                            multi_mv_uf_count += 1
                        if len(binPairUni) == 1:
                            multi_mv_ub_count += 1
            else:
                interactionType = "DUMP"
                ud_count += 1
        else:
            print("Warning! Read is either read1 nor read2")

    
    if args.summary:
        if args.summaryFile is None:
            outSummary = open(args.outdir + "/" + fileName + ".readCategorySummary", "w")
        else:
            outSummary = open(args.summaryFile, "a")
        outSummary.write("\n" + "\n")
        outSummary.write("Step: Categorize read pair type.\n" + "\n")
        outSummary.write("Total number of read pairs:\t" + str(uni_count + multi_count) + "\n")
        outSummary.write("  Uniquely mapping read pairs:\t" + str(uni_count) + "\n")
        outSummary.write("    Uniquely mapping valid read pairs:\t" + str(uv_count) + "\n")
        outSummary.write("    Uniquely mapping dangling end read pairs:\t" + str(de_count) + "\n")
        outSummary.write("    Uniquely mapping self cycle read pairs:\t" +str(sc_count) + "\n")
        outSummary.write("    Uniquely mapping religation read pairs:\t" + str(re_count) + "\n")
        outSummary.write("    Uniquely mapping dumped read pairs:\t" + str(ud_count) + "\n")
        outSummary.write("  Multi-mapping read pairs:\t" + str(multi_count) + "\n")
        outSummary.write("    Multi-mapping unique valid read pairs:\t" + str(multi_uv_count) + "\n")
        outSummary.write("    Multi-mapping multiple valid read pairs:\t" + str(multi_mv_count) + "\n")
        outSummary.write("    Multi-mapping multiple valid read pair but unique fragment pairs:\t" + str(multi_mv_uf_count) + "\n")
        outSummary.write("    Multi-mapping multiple valid read pair but unique bin pairs:\t" + str(multi_mv_ub_count) + "\n")
        outSummary.write("    Multi-mapping dumped read pairs:\t" + str(md_count) + "\n")

        outSummary.close()
    outContactNum.close()
    outValid.close()