<commit_before>#!/usr/bin/env python
## category General
## desc Calculate the inner mate-pair distance from two BAM files
'''
Calculate the inner mate-pair distance from two BAM files

With paired-end reads, knowing the inner mate-pair distance is key for some
downstream analyses or mapping. This value can be estimated using wet-lab
techniques, but can also be found empirically. If you map each fragment
separately to a common reference (genome or transcriptome), this command will
calculate the average distance between the fragments.

Reads where one (or both) fragment is unmapped are ignored, as are pairs that
map to different references.

(reference) ===============================================================
                 (frag 1) *********                      ******** (frag 2)
                                   |--------------------|
                                  inner mate-pair distance
'''
import sys
import os
from ngsutils.bam import bam_iter
from ngsutils.support.stats import counts_mean_stdev
import pysam


def bam_innerdist(bam1, bam2, out=sys.stdout):
    iter1 = bam_iter(bam1)
    iter2 = bam_iter(bam2, quiet=True)

    distances = {}
    total = 0
    proper = 0

    while True:
        try:
            read1 = iter1.next()
            read2 = iter2.next()
        except StopIteration:
            break

        if read1.qname != read2.qname:
            raise ValueError("Error: BAM files aren't properly paired! (%s, %s)\n" % (read1.qname, read2.qname))

        total += 1

        if read1.is_unmapped or read2.is_unmapped or read1.tid != read2.tid:
            continue

        proper += 1
        dist = read2.pos - read1.aend

        if not dist in distances:
            distances[dist] = 1
        else:<commit_msg>added code to allow for a read to map multiple places (still needs to be in correct sort order though)<commit_after>#!/usr/bin/env python
## category General
## desc Calculate the inner mate-pair distance from two BAM files
'''
Calculate the inner mate-pair distance from two BAM files

With paired-end reads, knowing the inner mate-pair distance is key for some
downstream analyses or mapping. This value can be estimated using wet-lab
techniques, but can also be found empirically. If you map each fragment
separately to a common reference (genome or transcriptome), this command will
calculate the average distance between the fragments.

Reads where one (or both) fragment is unmapped are ignored, as are pairs that
map to different references.

(reference) ===============================================================
                 (frag 1) *********                      ******** (frag 2)
                                   |--------------------|
                                  inner mate-pair distance
'''
import sys
import os
from ngsutils.bam import bam_iter
from ngsutils.support.stats import counts_mean_stdev
import pysam


def bam_innerdist(bam1, bam2, out=sys.stdout):
    iter1 = bam_iter(bam1)
    iter2 = bam_iter(bam2, quiet=True)

    distances = {}
    total = 0
    proper = 0

    read1_last = None
    read2_last = None
    read1 = None
    read2 = None

    while True:
        try:
            while not read1 or read1_last == read1.qname:
                read1 = iter1.next()
            while not read2 or read2_last == read2.qname:
                read2 = iter2.next()
        except StopIteration:
            break

        if read1.qname != read2.qname:
            raise ValueError("Error: BAM files aren't properly paired! (%s, %s)\n" % (read1.qname, read2.qname))

        read1_last = read1.qname
        read2_last = read2.qname

        total += 1

        if read1.is_unmapped or read2.is_unmapped or read1.tid != read2.tid:
            continue

        proper += 1
        dist = read2.pos - read1.aend

        if not dist in distances:
            distances[dist] = 1
        else:
{'fraction_non_alphanumeric': 0.11830985915492957, 'fraction_numerical': 0.01741357234314981, 'mean_word_length': 2.2604340567612686, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2979513', 'n_tokens_mistral': 1123, 'n_tokens_neox': 1058, 'n_words': 396}