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

"""
PBstat.py reads PBs frequency along the protein sequence
and computes Neq, distribution of PBs or logo representation of PBs.

2013 - <NAME>, <NAME> 
"""

#===============================================================================
# load modules
#===============================================================================
import PBlib as PB
import optparse 
# optparse in deprecated since Python 2.7 and has been replaced by argparse
# however many Python installations are steal using Python < 2.7
import os
import sys
import numpy 
import math
import subprocess

#===============================================================================
# functions
#===============================================================================
def array_to_string(ar):
    """convert numpy array to string for further import in R"""
    numpy.set_printoptions(threshold=numpy.inf)
    #  max_line_width : big enough to avoid unneeded newline
    #  precision : float with 4 digits
    return numpy.array_str(ar, max_line_width = 100000, precision = 4).translate(None, '[]')

#===============================================================================
# MAIN - program starts here
#===============================================================================

#-------------------------------------------------------------------------------
# manage parameters
#-------------------------------------------------------------------------------
parser = optparse.OptionParser(
    usage="%prog -f file.PB.count [options] -o output_root_name",
    version="1.0")
# mandatory arguments
mandatory_opts = optparse.OptionGroup(
    parser,
    'Mandatory arguments')
mandatory_opts.add_option("-f", action="store", type="string",
    help="name of file that contains PBs frequency (count)")
mandatory_opts.add_option("-o", action="store", type="string",
    help="root name for results")
parser.add_option_group(mandatory_opts)
# optional arguments
optional_opts = optparse.OptionGroup(
    parser,
    'Optional arguments')
optional_opts.add_option("--map", action="store_true", 
    dest = "mapdist", help = "generates map of the distribution of PBs along protein sequence")
optional_opts.add_option("--neq", action="store_true", 
    dest = "neq", help = "computes Neq and generates neq plot along protein sequence")
optional_opts.add_option("--logo", action="store_true", 
    dest = "logo", help = "generates logo representation of PBs frequency along protein sequence")
optional_opts.add_option("--residue-min", action="store", type="int",
    dest = "residue_min", help="defines lower bound of residue frame")
optional_opts.add_option("--residue-max", action="store", type="int",
    dest = "residue_max", help="defines upper bound of residue frame")
parser.add_option_group(optional_opts)
# get all parameters
(options, args) = parser.parse_args()

#  check options
if not options.f:
    parser.print_help()
    parser.error("option -f in mandatory")

if not options.o:
    parser.print_help()
    parser.error("option -o is mandatory")

#-------------------------------------------------------------------------------
# load and check data
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#  check count file is readable
if not os.path.isfile(options.f):
    sys.exit("ERROR: %s is not a valid file" % options.f)

# load count file 
# skip first row that contains PBs labels
try:
    freq = numpy.loadtxt(options.f, dtype=int, skiprows=1)
except:
    sys.exit("ERROR: wrong data format in %s" % options.f)

# check format
# 17 columns (residue number + 16 PBs) should be present
if len(freq[0,:]) != (PB.NUMBER + 1):
    sys.exit("ERROR: wrong data format in %s" % options.f)

# read residue numbers
residue_lst = list(freq[:, 0])

#-------------------------------------------------------------------------------
# check / define residue min / max
if options.residue_min:
    residue_min = options.residue_min
else:
    residue_min = min(freq[:,0])

if options.residue_max:
    residue_max = options.residue_max
else:
    residue_max = max(freq[:,0])

if residue_min < 0:
    sys.exit("ERROR: residue_min must be >= 0")

if residue_max < 0:
    sys.exit("ERROR: residue_max must be >= 0")

if residue_min >= residue_max:
    sys.exit("ERROR: residue_min must be > residue_max")

if residue_min not in residue_lst:
    sys.exit("ERROR: residue_min does not belong to the residue range in %s" % options.f)

if residue_max not in residue_lst:
    sys.exit("ERROR: residue_max does not belong to the residue range in %s" % options.f)


# get index of first residue
try:
    first_residue_index = int(freq[0, 0])
except:
    sys.exit("ERROR: cannot read index of first residue. Wrong data format in %s" % options.f)
print "Index of first residue is: ", first_residue_index

#  slice data to the required frame
freq = freq[residue_min - first_residue_index : residue_max - first_residue_index + 1, : ]

# determine number of sequences compiled
# use the sum of all residue at position 3
# since positions 1 and 2 have no PBs assignement
sequence_number = sum(freq[2, :])
if sequence_number == 0:
    sys.exit("ERROR: counting 0 sequences!")

# remove residue number    
# extract and normalize PBs frequencies 
freq = freq[:, 1:] / float(sequence_number)

#-------------------------------------------------------------------------------
# generates map of the distribution of PBs along protein sequence
#-------------------------------------------------------------------------------
if options.mapdist:
    #  convert array to string for further import in R
    freq_string = array_to_string(freq)

    # build R script
    #-------------------------------------------------------------------------------
    # data
    R_script="""
connector = textConnection("%s")
freq = read.table(connector, header = FALSE)
xnames = %d:%d<commit_msg>bug fix in computing of sequence_number, used to normalize PBs frequencies<commit_after>#! /usr/bin/env python

"""
PBstat.py reads PBs frequency along the protein sequence
and computes Neq, distribution of PBs or logo representation of PBs.

2013 - <NAME>, <NAME> 
"""

#===============================================================================
# load modules
#===============================================================================
import PBlib as PB
import optparse 
# optparse in deprecated since Python 2.7 and has been replaced by argparse
# however many Python installations are steal using Python < 2.7
import os
import sys
import numpy 
import math
import subprocess

#===============================================================================
# functions
#===============================================================================
def array_to_string(ar):
    """convert numpy array to string for further import in R"""
    numpy.set_printoptions(threshold=numpy.inf)
    #  max_line_width : big enough to avoid unneeded newline
    #  precision : float with 4 digits
    return numpy.array_str(ar, max_line_width = 100000, precision = 4).translate(None, '[]')

#===============================================================================
# MAIN - program starts here
#===============================================================================

#-------------------------------------------------------------------------------
# manage parameters
#-------------------------------------------------------------------------------
parser = optparse.OptionParser(
    usage="%prog -f file.PB.count [options] -o output_root_name",
    version="1.0")
# mandatory arguments
mandatory_opts = optparse.OptionGroup(
    parser,
    'Mandatory arguments')
mandatory_opts.add_option("-f", action="store", type="string",
    help="name of file that contains PBs frequency (count)")
mandatory_opts.add_option("-o", action="store", type="string",
    help="root name for results")
parser.add_option_group(mandatory_opts)
# optional arguments
optional_opts = optparse.OptionGroup(
    parser,
    'Optional arguments')
optional_opts.add_option("--map", action="store_true", 
    dest = "mapdist", help = "generates map of the distribution of PBs along protein sequence")
optional_opts.add_option("--neq", action="store_true", 
    dest = "neq", help = "computes Neq and generates neq plot along protein sequence")
optional_opts.add_option("--logo", action="store_true", 
    dest = "logo", help = "generates logo representation of PBs frequency along protein sequence")
optional_opts.add_option("--residue-min", action="store", type="int",
    dest = "residue_min", help="defines lower bound of residue frame")
optional_opts.add_option("--residue-max", action="store", type="int",
    dest = "residue_max", help="defines upper bound of residue frame")
parser.add_option_group(optional_opts)
# get all parameters
(options, args) = parser.parse_args()

#  check options
if not options.f:
    parser.print_help()
    parser.error("option -f in mandatory")

if not options.o:
    parser.print_help()
    parser.error("option -o is mandatory")

#-------------------------------------------------------------------------------
# load and check data
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
#  check count file is readable
if not os.path.isfile(options.f):
    sys.exit("ERROR: %s is not a valid file" % options.f)

# load count file 
# skip first row that contains PBs labels
try:
    freq = numpy.loadtxt(options.f, dtype=int, skiprows=1)
except:
    sys.exit("ERROR: wrong data format in %s" % options.f)

# check format
# 17 columns (residue number + 16 PBs) should be present
if len(freq[0,:]) != (PB.NUMBER + 1):
    sys.exit("ERROR: wrong data format in %s" % options.f)

# read residue numbers
residue_lst = list(freq[:, 0])

#-------------------------------------------------------------------------------
# check / define residue min / max
if options.residue_min:
    residue_min = options.residue_min
else:
    residue_min = min(freq[:,0])

if options.residue_max:
    residue_max = options.residue_max
else:
    residue_max = max(freq[:,0])

if residue_min < 0:
    sys.exit("ERROR: residue_min must be >= 0")

if residue_max < 0:
    sys.exit("ERROR: residue_max must be >= 0")

if residue_min >= residue_max:
    sys.exit("ERROR: residue_min must be > residue_max")

if residue_min not in residue_lst:
    sys.exit("ERROR: residue_min does not belong to the residue range in %s" % options.f)

if residue_max not in residue_lst:
    sys.exit("ERROR: residue_max does not belong to the residue range in %s" % options.f)


# get index of first residue
try:
    first_residue_index = int(freq[0, 0])
except:
    sys.exit("ERROR: cannot read index of first residue. Wrong data format in %s" % options.f)
print "Index of first residue is: ", first_residue_index

#  slice data to the required frame
freq = freq[residue_min - first_residue_index : residue_max - first_residue_index + 1, : ]

# determine number of sequences compiled
# use the sum of all residue at position 3
# since positions 1 and 2 have no PBs assignement
# and begin at 1 to not sum the index of the line (here is 3)
sequence_number = sum(freq[2, 1:])
if sequence_number == 0:
    sys.exit("ERROR: counting 0 sequences!")

# remove residue number    
# extract and normalize PBs frequencies
freq = freq[:, 1:] / float(sequence_number)

#-------------------------------------------------------------------------------
# generates map of the distribution of PBs along protein sequence
#-------------------------------------------------------------------------------
if options.mapdist:
    #  convert array to string for further import in R
    freq_string = array_to_string(freq)

    # build R script
    #-------------------------------------------------------------------------------
    # data
    R_script="""
connector = textConnection("%s")
freq = read.table(connector, header = FALSE)
xnames = %d:%d
{'fraction_non_alphanumeric': 0.2848189989280119, 'fraction_numerical': 0.007174074379483797, 'mean_word_length': 2.686322188449848, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 15, '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': '1100175', 'n_tokens_mistral': 3321, 'n_tokens_neox': 3077, 'n_words': 1035}