<commit_before>"""

.. module:: findparams
    :synopsis: Using reference annotation and data to find parameters for assembler.

..  moduleauthor:: <NAME> <<EMAIL>>

"""
# system imports
import subprocess
import multiprocessing
import gzip
import os
import time
from functools import reduce
from operator import iadd, iand
from collections import Counter
from itertools import repeat
import logging
logging.basicConfig(level=logging.DEBUG)
LOG = logging.getLogger(__name__)
import json

# 3rd party imports
import pandas as PD
import numpy as N
import matplotlib.pyplot as P

# library imports
from jgem import utils as UT
from jgem import bedtools as BT
from jgem import bigwig as BW
from jgem import gtfgffbed as GGB

from jgem.bigwig import BWObj #, BWs 

from sklearn.linear_model import LogisticRegression

def find_maxgap(arr, emin, emax, th, win, gapmode):
    if (emin>th):
        return 0
    if (emax<=th):
        return win
    idx0 = N.nonzero(arr>th) # first find furthest point
    maxpos = idx0[0][-1] # after this all zero
    idx = N.nonzero(arr<=th)
    if len(idx[0])==0:
        return 0
    if (gapmode!='i')&(idx[0][0]>maxpos):
        return 0
    cmax = 1
    cst = idx[0][0]
    ced = cst
    for i in idx[0][1:]:
        if (i>maxpos)&(gapmode!='i'):
            break        
        if i==ced+1: # continuous
            ced = i
        else:
            cmax = max(cmax, ced-cst+1)
            cst = ced = i
    cmax = max(cmax, ced-cst+1)
    return cmax

def find_firstgap(arr, emin, emax, th, win):
    if (emin>th):
        return 0, win # no gap
    if (emax<=th):
        return win, 0 # max
    idx = N.nonzero(arr<=th) # gap pos
    cst = idx[0][0] # current gap start
    ced = cst # current gap end
    for i in idx[0][1:]:
        if i==ced+1: # continuous
            ced = i
        else: # end of first gap and pos
            break
    return ced-cst+1, cst
      
def average(arr, n):
    end =  n * int(len(arr)/n)
    return N.mean(arr[:end].reshape(-1, n), 1)

class BWs(object):

    def __init__(self, paths):
        self.bwobjs = [BWObj(p) for p in paths]

    def __enter__(self):
        for b in self.bwobjs:
            b.__enter__()

    def __exit__(self, exc_type, exc_value, traceback):
        for b in self.bwobjs:
            b.__exit__(exc_type, exc_value, traceback)

    def get(self, chrom, st, ed):
        a = self.bwobjs[0].get(chrom, st, ed)
        for b in self.bwobjs[1:]:
            a += b.get(chrom, st, ed)
        return a

COPYCOLS = ['chr','st','ed','_gidx','locus','gene_id','gene_type']
CALCFLUXCOLS = ['_id', 'sdelta','ecovavg','ecovmin','ecovmax',
                'sin','sout','ein','eout','sdin','sdout']

class ParamFinder(object):
    """
    Args:
        refpre: pathprefix to ref (assume .ex.txt.gz, .sj.txt.gz)

    """
    def __init__(self, refpre, bwpre, refcode, genome, zoom=1, dstpre=None):
        self.refpre = refpre
        self.genome = genome
        self.bwpre = bwpre
        self.refcode = refcode
        if dstpre is None:
            dstpre = bwpre
        self.dstpre = dstpre
        self.zoom = zoom # multiply this factor to zoom in to small values
        self.ex = ex = UT.read_pandas(self.refpre+'.ex.txt.gz')
        self.sj = sj = UT.read_pandas(self.refpre+'.sj.txt.gz')
        self.set_bws(bwpre)

    def process(self, sdiffth=1, covfactor53=0, covfactorexon=0.05, np=10):
        self.extract_all()
        for x in ['ne_i','ne_5','ne_3','e5i','e3i','e53']:
            print('  #{0}:{1}'.format(x, len(getattr(self, x))))
        print('###### 53 params ########')
        self.calc_53_params(sdiffth=sdiffth, np=np)
        print('###### 53gap params ########')
        self.calc_53gap_params(covfactor=covfactor53, np=np)
        print('###### exon params ########')
        self.calc_exon_params(np=np, covfactor=covfactorexon)
        
    def set_bws(self, bwpre):
        self.bwpre = bwpre
        S2S = {'+':'.p','-':'.n','.':'.u'}
        self.bwpaths = bwp = {
            'ex': {s:bwpre+'.ex{0}.bw'.format(S2S[s]) for s in S2S},
            'sj': {s:bwpre+'.sj{0}.bw'.format(S2S[s]) for s in S2S},
        }
        self.bws = make_bws(bwp)
        
    def extract_all(self):
        self.extract_nonovl_exons()
        self.extract_exi53()
        self.extract_53_pair() # intergenic

    def extract_nonovl_exons(self):
        ex = self.ex
        sj = self.sj
        # nonovl exons    
        ex['gene_type'] = ex['extra'].str.split(';').str[2].str.split().str[1].str[1:-1]
        cols0 = ['chr','st','ed','_id']
        a = self.refpre+'.ex.bed.gz'
        a = UT.write_pandas(ex[cols0], a, '')
        b = self.refpre+'.sj.bed.gz'
        b = UT.write_pandas(sj[cols0], b, '')
        c1 = self.refpre+'.ex-ovl-sj.txt.gz'
        c2 = self.refpre+'.ex-ovl-ex.txt.gz'
        c1 = BT.bedtoolintersect(a,b,c1,wao=True)
        c2 = BT.bedtoolintersect(a,a,c2,wo=True)

        cols = cols0+['b_'+x for x in cols0]+['ovl']
        sov = UT.read_pandas(c1, names=cols)
        sov['len'] = sov['ed']-sov['st']
        sov['ovlratio'] = sov['ovl']/sov['len']
        sovg = sov.groupby('_id')['ovlratio'].max()
        snonov = sovg[sovg<1.] # not completely covered by junction

        eov = UT.read_pandas(c2, names=cols)
        eovsize = eov.groupby('_id').size()
        enonov = eovsize[eovsize==1] # only overlaps with self
        self.ne_i0 = ne_i0 = ex.set_index('_id').ix[enonov.index].sort_values(['chr','st','ed']).reset_index()
        self.ne_i0['len'] = ne_i0['ed']-ne_i0['st']

        LOG.info('#non-ex-ovl-ex={0}, #non-sj-ovl-ex={1}'.format(len(enonov), len(snonov)))
        ids = set(enonov.index).intersection(snonov.index)
        LOG.info('#non-ovl-ex={0}'.format(len(ids)))
        self.nov_ex = novex = ex.set_index('_id').ix[ids].sort_values(['chr','st','ed']).reset_index()
        novex['len'] = novex['ed']-novex['st']
        self.ne_i = novex[novex['cat']=='i']
        self.ne_5 = novex[novex['cat']=='5']
        self.ne_3 = novex[novex['cat']=='3']
        self.ne_s = novex[novex['cat']=='s']

    def extract_exi53(self):
        # internal exons overlapping with either 5 or 3 prime exons?
        cols0 = ['chr','st','ed','_id', 'sc1','strand']
        cols = cols0+['b_'+x for x in cols0]+['ovl']
        ex = self.ex

        exi = ex[ex['cat']=='i'] # internal exons
        ai = self.refpre + '.exi.bed.gz'
        ai = UT.write_pandas(exi[cols0], ai, '')

        e5 = ex[ex['cat']=='5']
        a5 = self.refpre + '.ex5.bed.gz'<commit_msg>ne_i0 (ok to overlap with junction) (ne_i did not get large Ttn exon)<commit_after>"""

.. module:: findparams
    :synopsis: Using reference annotation and data to find parameters for assembler.

..  moduleauthor:: <NAME> <<EMAIL>>

"""
# system imports
import subprocess
import multiprocessing
import gzip
import os
import time
from functools import reduce
from operator import iadd, iand
from collections import Counter
from itertools import repeat
import logging
logging.basicConfig(level=logging.DEBUG)
LOG = logging.getLogger(__name__)
import json

# 3rd party imports
import pandas as PD
import numpy as N
import matplotlib.pyplot as P

# library imports
from jgem import utils as UT
from jgem import bedtools as BT
from jgem import bigwig as BW
from jgem import gtfgffbed as GGB

from jgem.bigwig import BWObj #, BWs 

from sklearn.linear_model import LogisticRegression

def find_maxgap(arr, emin, emax, th, win, gapmode):
    if (emin>th):
        return 0
    if (emax<=th):
        return win
    idx0 = N.nonzero(arr>th) # first find furthest point
    maxpos = idx0[0][-1] # after this all zero
    idx = N.nonzero(arr<=th)
    if len(idx[0])==0:
        return 0
    if (gapmode!='i')&(idx[0][0]>maxpos):
        return 0
    cmax = 1
    cst = idx[0][0]
    ced = cst
    for i in idx[0][1:]:
        if (i>maxpos)&(gapmode!='i'):
            break        
        if i==ced+1: # continuous
            ced = i
        else:
            cmax = max(cmax, ced-cst+1)
            cst = ced = i
    cmax = max(cmax, ced-cst+1)
    return cmax

def find_firstgap(arr, emin, emax, th, win):
    if (emin>th):
        return 0, win # no gap
    if (emax<=th):
        return win, 0 # max
    idx = N.nonzero(arr<=th) # gap pos
    cst = idx[0][0] # current gap start
    ced = cst # current gap end
    for i in idx[0][1:]:
        if i==ced+1: # continuous
            ced = i
        else: # end of first gap and pos
            break
    return ced-cst+1, cst
      
def average(arr, n):
    end =  n * int(len(arr)/n)
    return N.mean(arr[:end].reshape(-1, n), 1)

class BWs(object):

    def __init__(self, paths):
        self.bwobjs = [BWObj(p) for p in paths]

    def __enter__(self):
        for b in self.bwobjs:
            b.__enter__()

    def __exit__(self, exc_type, exc_value, traceback):
        for b in self.bwobjs:
            b.__exit__(exc_type, exc_value, traceback)

    def get(self, chrom, st, ed):
        a = self.bwobjs[0].get(chrom, st, ed)
        for b in self.bwobjs[1:]:
            a += b.get(chrom, st, ed)
        return a

COPYCOLS = ['chr','st','ed','_gidx','locus','gene_id','gene_type']
CALCFLUXCOLS = ['_id', 'sdelta','ecovavg','ecovmin','ecovmax',
                'sin','sout','ein','eout','sdin','sdout']

class ParamFinder(object):
    """
    Args:
        refpre: pathprefix to ref (assume .ex.txt.gz, .sj.txt.gz)

    """
    def __init__(self, refpre, bwpre, refcode, genome, zoom=1, dstpre=None):
        self.refpre = refpre
        self.genome = genome
        self.bwpre = bwpre
        self.refcode = refcode
        if dstpre is None:
            dstpre = bwpre
        self.dstpre = dstpre
        self.zoom = zoom # multiply this factor to zoom in to small values
        self.ex = ex = UT.read_pandas(self.refpre+'.ex.txt.gz')
        self.sj = sj = UT.read_pandas(self.refpre+'.sj.txt.gz')
        self.set_bws(bwpre)

    def process(self, sdiffth=1, covfactor53=0, covfactorexon=0.05, np=10):
        self.extract_all()
        for x in ['ne_i','ne_5','ne_3','e5i','e3i','e53']:
            print('  #{0}:{1}'.format(x, len(getattr(self, x))))
        print('###### 53 params ########')
        self.calc_53_params(sdiffth=sdiffth, np=np)
        print('###### 53gap params ########')
        self.calc_53gap_params(covfactor=covfactor53, np=np)
        print('###### exon params ########')
        self.calc_exon_params(np=np, covfactor=covfactorexon)
        
    def set_bws(self, bwpre):
        self.bwpre = bwpre
        S2S = {'+':'.p','-':'.n','.':'.u'}
        self.bwpaths = bwp = {
            'ex': {s:bwpre+'.ex{0}.bw'.format(S2S[s]) for s in S2S},
            'sj': {s:bwpre+'.sj{0}.bw'.format(S2S[s]) for s in S2S},
        }
        self.bws = make_bws(bwp)
        
    def extract_all(self):
        self.extract_nonovl_exons()
        self.extract_exi53()
        self.extract_53_pair() # intergenic

    def extract_nonovl_exons(self):
        ex = self.ex
        sj = self.sj
        # nonovl exons    
        ex['gene_type'] = ex['extra'].str.split(';').str[2].str.split().str[1].str[1:-1]
        cols0 = ['chr','st','ed','_id']
        a = self.refpre+'.ex.bed.gz'
        a = UT.write_pandas(ex[cols0], a, '')
        b = self.refpre+'.sj.bed.gz'
        b = UT.write_pandas(sj[cols0], b, '')
        c1 = self.refpre+'.ex-ovl-sj.txt.gz'
        c2 = self.refpre+'.ex-ovl-ex.txt.gz'
        c1 = BT.bedtoolintersect(a,b,c1,wao=True)
        c2 = BT.bedtoolintersect(a,a,c2,wo=True)

        cols = cols0+['b_'+x for x in cols0]+['ovl']
        sov = UT.read_pandas(c1, names=cols)
        sov['len'] = sov['ed']-sov['st']
        sov['ovlratio'] = sov['ovl']/sov['len']
        sovg = sov.groupby('_id')['ovlratio'].max()
        snonov = sovg[sovg<1.] # not completely covered by junction

        eov = UT.read_pandas(c2, names=cols)
        eovsize = eov.groupby('_id').size()
        enonov = eovsize[eovsize==1] # only overlaps with self
        self.ne_i0 = ne_i0 = ex.set_index('_id').ix[enonov.index].sort_values(['chr','st','ed']).reset_index()
        self.ne_i0['len'] = ne_i0['ed']-ne_i0['st']
        LOG.info('#non-ovl-ex0={0}'.format(len(enonov)))

        LOG.info('#non-ex-ovl-ex={0}, #non-sj-ovl-ex={1}'.format(len(enonov), len(snonov)))
        ids = set(enonov.index).intersection(snonov.index)
        LOG.info('#non-ovl-ex={0}'.format(len(ids)))
        self.nov_ex = novex = ex.set_index('_id').ix[ids].sort_values(['chr','st','ed']).reset_index()
        novex['len'] = novex['ed']-novex['st']
        self.ne_i = novex[novex['cat']=='i']
        self.ne_5 = novex[novex['cat']=='5']
        self.ne_3 = novex[novex['cat']=='3']
        self.ne_s = novex[novex['cat']=='s']

    def extract_exi53(self):
        # internal exons overlapping with either 5 or 3 prime exons?
        cols0 = ['chr','st','ed','_id', 'sc1','strand']
        cols = cols0+['b_'+x for x in cols0]+['ovl']
        ex = self.ex

        exi = ex[ex['cat']=='i'] # internal exons
        ai = self.refpre + '.exi.bed.gz'
        ai = UT.write_pandas(exi[cols0], ai, '')

        e5 = ex[ex['cat']=='5']
        a5 = self.refpre + '.ex5.bed.gz'
{'fraction_non_alphanumeric': 0.17115954134710304, 'fraction_numerical': 0.01769306705140861, 'mean_word_length': 2.4915164369034994, 'pattern_counts': {'":': 0, '<': 19, '<?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': '2836245', 'n_tokens_mistral': 5156, 'n_tokens_neox': 4960, 'n_words': 1060}