<commit_before>"""

.. module:: repeats
    :synopsis:  Repeats (transposon) related stuffs

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

"""

import csv
import subprocess
import os
import gzip
import glob
import logging
logging.basicConfig(level=logging.DEBUG)
LOG = logging.getLogger(__name__)
import uuid

import pandas as PD
import numpy as N
import matplotlib.pylab as P

from jgem import utils as UT
from jgem import fasta as FA
from jgem import filenames as FN
from jgem import bedtools as BT
from jgem import gtfgffbed as GGB

from jgem import assembler2 as A2


RMSKPARAMS = dict(
    np = 4,
    th_uexon=4,
    th_bp_ovl=50,
    th_ex_ovl=50,
    datacode='',
    gname='gname',
)

class RmskFilter(object):
    """Filter genes with by overlap to repeat masker.

    Args:
        sjexpre: path prefix to assembled ex.txt.gz, sj.txt.gz files (optionally unionex.txt.gz )
        code: identifier
        chromdir: direcotry which contains chromosomes sequences in FASTA format
        rmskviz: RepeatMasker viz track (UCSC) converted in BED7 (using jgem.repeats.rmskviz2bed7)
        outdir: output directory

    """

    def __init__(self, sjexpre, code, chromdir, rmskviz, outdir, **kw):
        self.sjexpre = sjexpre
        prefix = os.path.join(outdir, code)
        self.fnobj = FN.FileNamesBase(prefix)
        self.chromdir = chromdir
        self.rmskviz = rmskviz
        self.gfc = FA.GenomeFASTAChroms(chromdir)

        self.params = RMSKPARAMS.copy()
        self.params.update(kw)


        self.ex = UT.read_pandas(sjexpre+'.ex.txt.gz')
        self.sj = UT.read_pandas(sjexpre+'.sj.txt.gz')
        if 'glen' not in self.ex or 'tlen' not in self.ex:
            if not os.path.exists(sjexpre+'.ci.txt.gz'):
                ci = UT.chopintervals(ex, sjexpre+'.ci.txt.gz')
            else:
                ci = UT.read_pandas(sjexpre+'.ci.txt.gz')
            UT.set_glen_tlen(ex,ci,gidx='_gidx')
            UT.write_pandas(sjexpre+'.ex.txt.gz')
        uexpath = sjexpre+'.unionex.txt.gz'        
        if os.path.exists(uexpath):
            self.uex = UT.read_pandas(uexpath)
        else:
            LOG.info('making union exons...saving to {0}'.format(uexpath))
            self.uex = UT.make_unionex(self.ex, '_gidx')
            UT.write_pandas(self.uex, uexpath, 'h')


    def calculate(self):
        """ Calculate base pair overlap to repeat using UCSC genome mask of repeats to lower case, 
        and exon level overlap to repeat using UCSC RepeatMaskerViz track. 
        ALso make a dataframe containing summary. 
        """
        pr = self.params
        fn = self.fnobj

        uex = count_repeats_mp(self.uex, self.gfc, np=pr['np'], col='#repbp')
        uex = count_repeats_viz_mp(uex, self.rmskviz, np=pr['np'], idcol='_id', expand=0, col='repnames')
        self.ugb = ugb = self._make_gbed(self.ex, self.sj, uex, datacode=pr['datacode'], gname=pr['gname'])
        UT.write_pandas(ugb, fn.txtname('all.genes.stats', category='output'), 'h')


    def _make_gbed(self, ex, sj, ugb, datacode='', gname='gname'):
        # rep%
        if 'glen' not in ugb:

        gr = ugb.groupby('_gidx')
        gb2 = gr[['chr',gname,'tlen','glen']].first()
        gb2['#repbp'] = gr['#repbp'].sum()
        gb2['rep%'] = 100.*gb2['#repbp']/gb2['tlen']
        # rmskviz, exon%
        gb2['#uexons'] = gr.size()
        gbsub = ugb[ugb['repnames']!='.(-1)'] # .(-1) == overlap
        gb2['#uexons_rmsk'] = gbsub.groupby('_gidx').size() # num exons overlapping rmskviz
        gb2.loc[gb2['#uexons_rmsk'].isnull(),'#uexons_rmsk'] = 0 
        gb2['rviz%'] = 100.*gb2['#uexons_rmsk']/gb2['#uexons']
        gb2['repnames'] = gbsub.groupby('_gidx')['repnames'].apply(lambda x: ';'.join(list(x)))
        # locus
        gb2['st'] = gr['st'].min()
        gb2['ed'] = gr['ed'].max()
        gb2['glocus'] = UT.calc_locus(gb2,'chr','st','ed')
        # rmskviz, class=[Simple_repeats, LINE, SINE, LTR, DNA]
        rcols = ['Simple_repeat','LINE','SINE','LTR','DNA']
        for k in rcols:
            gb2[k] = gb2['repnames'].str.contains('#'+k)

        dc = '_'+datacode if datacode else ''
        egr = ex.groupby('_gidx')
        gb2['#exons'] = egr.size()
        gb2['avgecov'] = egr['ecov'+dc].mean()
        gb2['gcov'] = egr['gcov'+dc].first() 

        sgr = sj.groupby('_gidx')
        gb2['ucnt'] = sgr['ucnt'+dc].sum()<commit_msg>glen tlen<commit_after>"""

.. module:: repeats
    :synopsis:  Repeats (transposon) related stuffs

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

"""

import csv
import subprocess
import os
import gzip
import glob
import logging
logging.basicConfig(level=logging.DEBUG)
LOG = logging.getLogger(__name__)
import uuid

import pandas as PD
import numpy as N
import matplotlib.pylab as P

from jgem import utils as UT
from jgem import fasta as FA
from jgem import filenames as FN
from jgem import bedtools as BT
from jgem import gtfgffbed as GGB

from jgem import assembler2 as A2


RMSKPARAMS = dict(
    np = 4,
    th_uexon=4,
    th_bp_ovl=50,
    th_ex_ovl=50,
    datacode='',
    gname='gname',
)

class RmskFilter(object):
    """Filter genes with by overlap to repeat masker.

    Args:
        sjexpre: path prefix to assembled ex.txt.gz, sj.txt.gz files (optionally unionex.txt.gz )
        code: identifier
        chromdir: direcotry which contains chromosomes sequences in FASTA format
        rmskviz: RepeatMasker viz track (UCSC) converted in BED7 (using jgem.repeats.rmskviz2bed7)
        outdir: output directory

    """

    def __init__(self, sjexpre, code, chromdir, rmskviz, outdir, **kw):
        self.sjexpre = sjexpre
        prefix = os.path.join(outdir, code)
        self.fnobj = FN.FileNamesBase(prefix)
        self.chromdir = chromdir
        self.rmskviz = rmskviz
        self.gfc = FA.GenomeFASTAChroms(chromdir)

        self.params = RMSKPARAMS.copy()
        self.params.update(kw)


        self.ex = UT.read_pandas(sjexpre+'.ex.txt.gz')
        self.sj = UT.read_pandas(sjexpre+'.sj.txt.gz')
        if 'glen' not in self.ex or 'tlen' not in self.ex:
            if not os.path.exists(sjexpre+'.ci.txt.gz'):
                ci = UT.chopintervals(ex, sjexpre+'.ci.txt.gz')
            else:
                ci = UT.read_pandas(sjexpre+'.ci.txt.gz')
            UT.set_glen_tlen(ex,ci,gidx='_gidx')
            UT.write_pandas(sjexpre+'.ex.txt.gz')
        uexpath = sjexpre+'.unionex.txt.gz'        
        if os.path.exists(uexpath):
            self.uex = UT.read_pandas(uexpath)
        else:
            LOG.info('making union exons...saving to {0}'.format(uexpath))
            self.uex = UT.make_unionex(self.ex, '_gidx')
            UT.write_pandas(self.uex, uexpath, 'h')


    def calculate(self):
        """ Calculate base pair overlap to repeat using UCSC genome mask of repeats to lower case, 
        and exon level overlap to repeat using UCSC RepeatMaskerViz track. 
        ALso make a dataframe containing summary. 
        """
        pr = self.params
        fn = self.fnobj

        uex = count_repeats_mp(self.uex, self.gfc, np=pr['np'], col='#repbp')
        uex = count_repeats_viz_mp(uex, self.rmskviz, np=pr['np'], idcol='_id', expand=0, col='repnames')
        self.ugb = ugb = self._make_gbed(self.ex, self.sj, uex, datacode=pr['datacode'], gname=pr['gname'])
        UT.write_pandas(ugb, fn.txtname('all.genes.stats', category='output'), 'h')


    def _make_gbed(self, ex, sj, ugb, datacode='', gname='gname'):
        # rep%
        gr = ugb.groupby('_gidx')
        gb2 = gr[['chr',gname,'tlen','glen']].first()
        gb2['#repbp'] = gr['#repbp'].sum()
        gb2['rep%'] = 100.*gb2['#repbp']/gb2['tlen']
        # rmskviz, exon%
        gb2['#uexons'] = gr.size()
        gbsub = ugb[ugb['repnames']!='.(-1)'] # .(-1) == overlap
        gb2['#uexons_rmsk'] = gbsub.groupby('_gidx').size() # num exons overlapping rmskviz
        gb2.loc[gb2['#uexons_rmsk'].isnull(),'#uexons_rmsk'] = 0 
        gb2['rviz%'] = 100.*gb2['#uexons_rmsk']/gb2['#uexons']
        gb2['repnames'] = gbsub.groupby('_gidx')['repnames'].apply(lambda x: ';'.join(list(x)))
        # locus
        gb2['st'] = gr['st'].min()
        gb2['ed'] = gr['ed'].max()
        gb2['glocus'] = UT.calc_locus(gb2,'chr','st','ed')
        # rmskviz, class=[Simple_repeats, LINE, SINE, LTR, DNA]
        rcols = ['Simple_repeat','LINE','SINE','LTR','DNA']
        for k in rcols:
            gb2[k] = gb2['repnames'].str.contains('#'+k)

        dc = '_'+datacode if datacode else ''
        egr = ex.groupby('_gidx')
        gb2['#exons'] = egr.size()
        gb2['avgecov'] = egr['ecov'+dc].mean()
        gb2['gcov'] = egr['gcov'+dc].first() 

        sgr = sj.groupby('_gidx')
        gb2['ucnt'] = sgr['ucnt'+dc].sum()
{'fraction_non_alphanumeric': 0.15637530072173217, 'fraction_numerical': 0.010310459388246077, 'mean_word_length': 2.6588432523051133, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, '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': '2172688', 'n_tokens_mistral': 3431, 'n_tokens_neox': 3255, 'n_words': 671}