<commit_before>#!/usr/bin/env python
#
# -----------------------------------------------------------------------------
# Copyright (c) 2016   <NAME> <<EMAIL>>
# Copyright (c) 2016   Indiana University
#
# This file is part of genhub (http://github.com/standage/genhub) and is
# licensed under the BSD 3-clause license: see LICENSE.txt.
# -----------------------------------------------------------------------------

from __future__ import print_function
from __future__ import division
import argparse
import re
import sys
import genhub


def cli():
    """Define the command-line interface of the program."""
    desc = 'Calculate measures of "compactness" for the specified genome(s).'
    parser = argparse.ArgumentParser(description=desc)
    parser.add_argument('-v', '--version', action='version',
                        version='GenHub v%s' % genhub.__version__)
    parser.add_argument('-c', '--cfgdir', default=None, metavar='DIR',
                        help='directory (or comma-separated list of '
                        'directories) from which to load user-supplied genome '
                        'configuration files')
    parser.add_argument('-w', '--workdir', metavar='WD', default='./species',
                        help='working directory for data files; default is '
                        '"./species"')
    parser.add_argument('-l', '--length', metavar='LEN', type=int,
                        default=1000000, help='minimum length threshold; '
                        'default is 1000000 (1Mb)')
    parser.add_argument('species', nargs='+', help='species label(s)')
    return parser


def ilocus_type(record):
    typematch = re.search('iLocus_type=([^;\n]+)', record)
    assert typematch
    return typematch.group(1)


def seq_info(db):
    with open(db.ilocusfile) as infile:
        seqlengths = dict()
        gilocus_counts = dict()
        for line in infile:
            if line.startswith('##sequence-region'):
                seqreg = re.search('##sequence-region\s+(\S+)\s+(\d+)\s+(\d+)',
                                   line)
                assert seqreg
                seqid = seqreg.group(1)
                start = int(seqreg.group(2))
                end = int(seqreg.group(3))
                assert start == 1
                seqlengths[seqid] = end
            elif '\tlocus\t' in line:
                values = line.split('\t')
                assert len(values) == 9
                assert values[2] == 'locus'
                seqid = values[0]
                attrs = values[8]
                locustype = ilocus_type(attrs)
                if locustype in ['siLocus', 'ciLocus', 'niLocus']:
                    if seqid not in gilocus_counts:
                        gilocus_counts[seqid] = 0
                    gilocus_counts[seqid] += 1

    return seqlengths, gilocus_counts


def milocus_info(db):
    with open(db.milocusfile) as infile:
        milocus_occ = dict()
        singleton_count = dict()
        for line in infile:
            values = line.split('\t')
            assert len(values) == 9
            assert values[2] == 'locus'
            seqid = values[0]
            start = int(values[3]) - 1
            end = int(values[4])
            attrs = values[8]
            locustype = ilocus_type(attrs)
            if locustype == 'miLocus':
                length = end - start
                if seqid not in milocus_occ:
                    milocus_occ[seqid] = 0
                milocus_occ[seqid] += length
            elif locustype in ['siLocus', 'ciLocus', 'niLocus']:
                if seqid not in singleton_count:
                    singleton_count[seqid] = 0
                singleton_count[seqid] += 1
    print(milocus_occ)
    return milocus_occ, singleton_count


def main(args):
    print('Species', 'SeqID', 'Sigma', 'Phi', sep='\t')

    registry = genhub.registry.Registry()
    if args.cfgdir:
        for cfgdirpath in args.cfgdir.split(','):
            registry.update(cfgdirpath)
    conf = registry.genomes(args.species)

    for species in args.species:
        config = conf[species]
        db = genhub.genomedb.GenomeDB(species, config, workdir=args.workdir)
        seqlengths, gilocus_counts = seq_info(db)
        milocus_occ, singleton_count = milocus_info(db)
        for seqid in seqlengths:
            length = seqlengths[seqid]
            if length < args.length:
                continue
            milocus_space = 0
            if seqid in milocus_occ:
                milocus_space = milocus_occ[seqid]
            giloci = gilocus_counts[seqid]
            singletons = 0
            if seqid in singleton_count:
                singletons = singleton_count[seqid]
            nonsingletons = giloci - singletons
            sigma = milocus_space / length
            phi = nonsingletons / giloci

            print(species, seqid, sigma, phi, sep='\t')

if __name__ == '__main__':
    main(args=cli().parse_args())<commit_msg>Improved compactness script, filtering implemented<commit_after>#!/usr/bin/env python
#
# -----------------------------------------------------------------------------
# Copyright (c) 2016   <NAME> <<EMAIL>>
# Copyright (c) 2016   Indiana University
#
# This file is part of genhub (http://github.com/standage/genhub) and is
# licensed under the BSD 3-clause license: see LICENSE.txt.
# -----------------------------------------------------------------------------

from __future__ import print_function
from __future__ import division
import argparse
import pandas
import re
import sys
import genhub


def cli():
    """Define the command-line interface of the program."""
    desc = 'Calculate measures of compactness for the specified genome(s).'
    parser = argparse.ArgumentParser(description=desc)
    parser.add_argument('-v', '--version', action='version',
                        version='GenHub v%s' % genhub.__version__)
    parser.add_argument('-c', '--cfgdir', default=None, metavar='DIR',
                        help='directory (or comma-separated list of '
                        'directories) from which to load user-supplied genome '
                        'configuration files')
    parser.add_argument('-w', '--workdir', metavar='WD', default='./species',
                        help='working directory for data files; default is '
                        '"./species"')
    parser.add_argument('-l', '--length', metavar='LEN', type=int,
                        default=1000000, help='minimum length threshold; '
                        'default is 1000000 (1Mb)')
    parser.add_argument('species', nargs='+', help='species label(s)')
    return parser


def thresholds(seqid, iloci, iqnt=0.95, gqnt=0.05):
    seqloci = iloci.loc[iloci.SeqID == seqid]
    ithresh = None
    if iqnt:
        iiloci = seqloci.loc[seqloci.LocusClass == 'iiLocus']
        ithresh = int(iiloci['Length'].quantile(iqnt))
    gthresh = None
    if gqnt:
        gilocus_types = ['siLocus', 'ciLocus', 'niLocus']
        giloci = seqloci.loc[seqloci.LocusClass.isin(gilocus_types)]
        gthresh = int(giloci['Length'].quantile(iqnt))
    return ithresh, gthresh


def seqlen(seqid, iloci, ithresh=None, gthresh=None):
    seqloci = iloci.loc[(iloci.SeqID == seqid) &
                        (seqloci.LocusClass != 'fiLocus')]
    effsize = seqloci['EffectiveLength'].sum()
    if ithresh:
        longiiloci = seqloci.loc[(seqloci.LocusClass == 'iiLocus') &
                                 (seqloci.Length > ithresh)]
        effsize -= longiiloci['Length'].sum()
    if gthresh:
        gilocus_types = ['siLocus', 'ciLocus', 'niLocus']
        shortgiloci = seqloci.loc[(seqloci.LocusClass.isin(gilocus_types)) &
                                  (seqloci.Length < gthresh)]
        effsize -= shortgiloci['Length'].sum()
    return effsize


def calc_phi(seqid, iloci, miloci):
    gilocus_types = ['siLocus', 'ciLocus', 'niLocus']
    giloci = iloci.loc[(iloci.SeqID == seqid) &
                       (iloci.LocusClass.isin(gilocus_types))]
    singletons = miloci.loc[(miloci.SeqID == seqid) &
                            (miloci.LocusClass.isin(gilocus_types))]
    if gthresh:
        giloci = giloci.loc[giloci.Length >= gthresh]
        singletons = singletons.loc[singletons.Length >= gthresh]
    merged = len(giloci) - len(singletons)
    return merged / len(giloci)


def main(args):
    print('Species', 'SeqID', 'Sigma', 'Phi', sep='\t')

    registry = genhub.registry.Registry()
    if args.cfgdir:
        for cfgdirpath in args.cfgdir.split(','):
            registry.update(cfgdirpath)
    conf = registry.genomes(args.species)

    for species in args.species:
        config = conf[species]
        db = genhub.genomedb.GenomeDB(species, config, workdir=args.workdir)
        iloci = pandas.read_table(db.ilocustable)
        miloci = pandas.read_table(db.milocustable)
        for seqid in iloci['SeqID']:
            ithresh, gthresh = thresholds(seqid, iloci)
            length = seqlen(seqid, iloci, ithresh, gthresh)
            phi = calc_phi(seqid, iloci, miloci)
            milocus_occ = miloci.loc[miloci.SeqID == seqid]['Length'].sum()
            sigma = milocus_occ / length
            print(species, seqid, sigma, phi, sep='\t')

if __name__ == '__main__':
    main(args=cli().parse_args())
{'fraction_non_alphanumeric': 0.15041962556488056, 'fraction_numerical': 0.008284914998924036, 'mean_word_length': 1.9498571881942242, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7288159', 'n_tokens_mistral': 2917, 'n_tokens_neox': 2713, 'n_words': 628}