<commit_before>#!/usr/bin/env python

from common import *
from collections import OrderedDict as od
import subprocess


def countBaseAtPos(bamfile,chrom,pos,mutid='null'):
    """ return list of bases at position chrom,pos
    """
    locstr = chrom + ":" + str(pos) + "-" + str(pos)
    args = ['samtools','mpileup',bamfile,'-r',locstr]

    p = subprocess.Popen(args,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
    p.wait()
    pout = p.stdout.readlines()

    pileup = None 

    for line in pout:
        try:
            c = line.strip().split()
            assert len(c) > 5
            pileup = c[4].upper()
        except AssertionError:
            sys.stderr.write("INFO\t" + now() + "\t" + mutid + "\tmpileup failed, no coverage for base: " + chrom + ":" + str(pos) + "\n")
            return []
    bases = []
    if pileup:
        for b in pileup:
            if b in ['A','T','C','G']:
                bases.append(b)

    return bases


def makeins(read, start, ins, debug=False):
    assert len(read.seq) > len(ins) + 2
    
    if debug:
        print "DEBUG: INS: read.pos:", read.pos
        print "DEBUG: INS: start:   ", start
        print "DEBUG: INS: ins:     ", ins
        print "DEBUG: DEL: cigar:     ", read.cigarstring

    orig_len = len(read.seq)
    pos_in_read = start - read.pos + read.qstart


    if pos_in_read > 0: # insertion start in read span
        if debug:
            print "DEBUG: INS: pos_in_read:", pos_in_read
        left  = read.seq[:pos_in_read]
        right = read.seq[pos_in_read:]

        newseq = left + ins + right
        newseq = newseq[:orig_len]

    else: # insertion continues to the left of read
        right = read.seq[pos_in_read:]
        newseq = ins + right
        newseq = newseq[-orig_len:]

    if debug:
        print "DEBUG: INS: orig seq:", read.seq
        print "DEBUG: INS: newseq:  ", newseq
    return newseq


def makedel(read, chrom, start, end, ref, debug=False):
    assert len(read.seq) > end-start-2
    
    if debug:
        print "DEBUG: DEL: read.pos:", read.pos
        print "DEBUG: DEL: start:   ", start
        print "DEBUG: DEL: ins:     ", end
        print "DEBUG: DEL: cigar:     ", read.cigarstring
        print "DEBUG: DEL: orig seq:     ", read.seq

    orig_len = len(read.seq)
    #orig_end = read.pos + orig_len
    start_in_read = start - read.pos + read.qstart
    end_in_read = end - read.pos + read.qstart

    if debug:
        print "DEBUG: DEL: start_in_read:", start_in_read
        print "DEBUG: DEL: end_in_read:  ", end_in_read

    if start_in_read < 0: # deletion begins to the left of the read
        if debug:
            print "DEBUG: DEL: del begins to left of read." 

        assert end_in_read < orig_len
        right = read.seq[end_in_read:]
        left  = ref.fetch(chrom, start-(len(read.seq) - len(right)), start)

    elif end_in_read > orig_len: # deletion ends to the right of the read
        if debug:
            print "DEBUG: DEL: del ends to right of read."

        assert start_in_read > 0
        left  = read.seq[:start_in_read]
        right = ref.fetch(chrom, end, end+(len(read.seq) - len(left)))

    else:
        if debug:
            print "DEBUG: DEL: del starts and ends within read." 

        assert end_in_read <= orig_len and start_in_read >= 0 # deletion contained within the read
        left  = read.seq[:start_in_read]
        right = read.seq[end_in_read:]
        right += ref.fetch(chrom, read.pos+len(read.seq), read.pos+len(read.seq)+(len(read.seq)-len(left)-len(right)))

    if debug:
        print "DEBUG: DEL:  newseq:     ", left + right
    return left + right


def find_mate(read, bam):
    ''' AlignmentFile.mate() can return a non-primary alignment, so use this function instead '''
    chrom = bam.getrname(read.tid)
    for rec in bam.fetch(chrom, read.pnext, read.pnext+1):
        if rec.qname == read.qname and rec.pos == read.pnext:
            if not rec.is_secondary and bin(rec.flag & 2048) != bin(2048):
                if rec.is_read1 != read.is_read1:
                    return rec
    return None


def mutate(args, log, bamfile, bammate, chrom, mutstart, mutend, mutpos_list, avoid=None, mutid_list=None, is_snv=False, mutbase_list=None, is_insertion=False, is_deletion=False, ins_seq=None, reffile=None, indel_start=None, indel_end=None):
    assert mutend > mutstart, "mutation start must occur before mutation end: " + mutid

    hasSNP = False

    outreads = od()
    mutreads = od()
    mutmates = od()

    region = 'haplo_' + chrom + '_' + str(mutstart) + '_' + str(mutend)<commit_msg>update how read mate is found to allow for mates on different chromosomes<commit_after>#!/usr/bin/env python

from common import *
from collections import OrderedDict as od
import subprocess


def countBaseAtPos(bamfile,chrom,pos,mutid='null'):
    """ return list of bases at position chrom,pos
    """
    locstr = chrom + ":" + str(pos) + "-" + str(pos)
    args = ['samtools','mpileup',bamfile,'-r',locstr]

    p = subprocess.Popen(args,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
    p.wait()
    pout = p.stdout.readlines()

    pileup = None 

    for line in pout:
        try:
            c = line.strip().split()
            assert len(c) > 5
            pileup = c[4].upper()
        except AssertionError:
            sys.stderr.write("INFO\t" + now() + "\t" + mutid + "\tmpileup failed, no coverage for base: " + chrom + ":" + str(pos) + "\n")
            return []
    bases = []
    if pileup:
        for b in pileup:
            if b in ['A','T','C','G']:
                bases.append(b)

    return bases


def makeins(read, start, ins, debug=False):
    assert len(read.seq) > len(ins) + 2
    
    if debug:
        print "DEBUG: INS: read.pos:", read.pos
        print "DEBUG: INS: start:   ", start
        print "DEBUG: INS: ins:     ", ins
        print "DEBUG: DEL: cigar:     ", read.cigarstring

    orig_len = len(read.seq)
    pos_in_read = start - read.pos + read.qstart


    if pos_in_read > 0: # insertion start in read span
        if debug:
            print "DEBUG: INS: pos_in_read:", pos_in_read
        left  = read.seq[:pos_in_read]
        right = read.seq[pos_in_read:]

        newseq = left + ins + right
        newseq = newseq[:orig_len]

    else: # insertion continues to the left of read
        right = read.seq[pos_in_read:]
        newseq = ins + right
        newseq = newseq[-orig_len:]

    if debug:
        print "DEBUG: INS: orig seq:", read.seq
        print "DEBUG: INS: newseq:  ", newseq
    return newseq


def makedel(read, chrom, start, end, ref, debug=False):
    assert len(read.seq) > end-start-2
    
    if debug:
        print "DEBUG: DEL: read.pos:", read.pos
        print "DEBUG: DEL: start:   ", start
        print "DEBUG: DEL: ins:     ", end
        print "DEBUG: DEL: cigar:     ", read.cigarstring
        print "DEBUG: DEL: orig seq:     ", read.seq

    orig_len = len(read.seq)
    #orig_end = read.pos + orig_len
    start_in_read = start - read.pos + read.qstart
    end_in_read = end - read.pos + read.qstart

    if debug:
        print "DEBUG: DEL: start_in_read:", start_in_read
        print "DEBUG: DEL: end_in_read:  ", end_in_read

    if start_in_read < 0: # deletion begins to the left of the read
        if debug:
            print "DEBUG: DEL: del begins to left of read." 

        assert end_in_read < orig_len
        right = read.seq[end_in_read:]
        left  = ref.fetch(chrom, start-(len(read.seq) - len(right)), start)

    elif end_in_read > orig_len: # deletion ends to the right of the read
        if debug:
            print "DEBUG: DEL: del ends to right of read."

        assert start_in_read > 0
        left  = read.seq[:start_in_read]
        right = ref.fetch(chrom, end, end+(len(read.seq) - len(left)))

    else:
        if debug:
            print "DEBUG: DEL: del starts and ends within read." 

        assert end_in_read <= orig_len and start_in_read >= 0 # deletion contained within the read
        left  = read.seq[:start_in_read]
        right = read.seq[end_in_read:]
        right += ref.fetch(chrom, read.pos+len(read.seq), read.pos+len(read.seq)+(len(read.seq)-len(left)-len(right)))

    if debug:
        print "DEBUG: DEL:  newseq:     ", left + right
    return left + right


def find_mate(read, bam):
    ''' AlignmentFile.mate() can return a non-primary alignment, so use this function instead '''
    chrom = read.next_reference_name
    for rec in bam.fetch(chrom, read.next_reference_start, read.next_reference_start+1):
        if rec.query_name == read.query_name and rec.reference_start == read.next_reference_start:
            if not rec.is_secondary and bin(rec.flag & 2048) != bin(2048):
                if rec.is_read1 != read.is_read1:
                    return rec
    return None


def mutate(args, log, bamfile, bammate, chrom, mutstart, mutend, mutpos_list, avoid=None, mutid_list=None, is_snv=False, mutbase_list=None, is_insertion=False, is_deletion=False, ins_seq=None, reffile=None, indel_start=None, indel_end=None):
    assert mutend > mutstart, "mutation start must occur before mutation end: " + mutid

    hasSNP = False

    outreads = od()
    mutreads = od()
    mutmates = od()

    region = 'haplo_' + chrom + '_' + str(mutstart) + '_' + str(mutend)
{'fraction_non_alphanumeric': 0.12219718914279584, 'fraction_numerical': 0.004076815792296964, 'mean_word_length': 2.399708242159008, 'pattern_counts': {'":': 4, '<': 9, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7353586', 'n_tokens_mistral': 3201, 'n_tokens_neox': 3024, 'n_words': 857}