New upstream version 3.0.2

# Change log

All notable changes to this project will be documented in this file.

## Version 3.0.2 2019-12-21

#### TRANSIT:

 - Mostly cosmetic fixes

 - Updated some command-line and GUI messages

 - Updated documentation (especially for GI and resampling)

 - Removed "warning: high stderr" from gene status in ZINB

 - Added LFCs in ZINB output

## Version 3.0.1 2019-08-01

#### TRANSIT:

 - Add check for python3 (TRANSIT 3+ requires python3.6+)

import sys,random,numpy

import pytransit.tnseq_tools

from statsmodels.stats.multitest import fdrcorrection

def read_wig(fname):

  coords,counts = [],[]

  for line in open(fname):

    if line[0] not in "0123456789": continue

    w = line.rstrip().split()

    coord,cnt = int(w[0]),int(w[1])

    coords.append(coord)

    counts.append(cnt)

  return coords,counts

# remove all runs of zeros of length >= W

def remove_essential_regions(wig,W):

  runs = []

  i,n = 0,len(wig)

  while i<n:

    if wig[i]>0: i += 1

    else:

      j = i

      while j<n and wig[j]==0: j += 1

      runs.append((i,j))

      i = j

  counts = []

  for k,(i,j) in enumerate(runs):

    if j-i<W: counts += wig[i:j]

    if k<len(runs)-1:

      next = runs[k+1][0]

      counts += wig[j:next]

  return counts

def get_counts(coords,counts,gene):

  cnts = []

  start,end = gene['start'],gene['end']

  for i,co in enumerate(coords):

    if co>=start and co<=end: cnts.append(counts[i])

  return cnts

def sample_counts(counts,size,times):

  samples = []

  for i in range(times):

    samples.append(random.sample(counts,size)) # without replacement

  return samples

###################################

if len(sys.argv)<3:

  print "usage: python fitness_defect.py <comma_separated_list_of_wig_files> <prot_table>"

  sys.exit(0)

print "# command: python",

for x in sys.argv: print x,

print

coords,counts = read_wig(sys.argv[1])

genes = pytransit.tnseq_tools.read_genes(sys.argv[2])

noness = remove_essential_regions(counts,5) # run length

print '# sites: %s, noness: %s' % (len(counts),len(noness))

noness_arr = numpy.array((noness))

noness_NZvals = noness_arr[noness_arr>0]

print '# noness: zeros=%s, NZmean=%0.1f' % (noness_NZvals.size,numpy.mean(noness_NZvals))

cache = {}

results = []

for gene in genes:

  sys.stderr.write("%s %s\n" % (gene['rv'],gene['gene']))

  cnts = get_counts(coords,counts,gene)

  n = len(cnts)

  if n==0: continue

  nonzeros = [cnts[x] for x in numpy.nonzero(cnts)[0]]

  zeros,NZmean = n-len(nonzeros),numpy.mean(nonzeros) if len(nonzeros)>0 else 0

  sat = len(nonzeros)/float(n)

  tot = sum(cnts)

  mn = tot/float(n)

  # determine p-value by comparing to sum of counts for random draws of sites from noness

  N,alpha = 10000,0.05

  if n in cache: sample = cache[n]

  else: 

    sample = sample_counts(noness,n,N)

    cache[n] = sample

  samplesums = [sum(lst) for lst in sample]

  meansum = numpy.mean(samplesums)

  PC = 1 # pseudo-counts

  rel = (sum(cnts)+PC)/float(meansum+PC)

  LFC = numpy.log2(rel)

  lesser = len(list(filter(lambda x: x<=tot,samplesums)))

  greater = len(list(filter(lambda x: x>=tot,samplesums)))

  pval = min(lesser,greater)/float(N)

  vals = [gene[x] for x in "rv gene start end strand".split()]

  vals += [len(cnts),"%s" % tot,"%0.1f" % mn,"%0.3f" % sat,len(nonzeros),"%0.1f" % NZmean]

  vals += [int(meansum),"%0.3f" % rel,"%0.3f" % LFC,pval]

  results.append(vals)

  #if gene['rv']=='Rv0020c': break

pvals = [x[-1] for x in results]

qvals = list(fdrcorrection(pvals)[1])

results = [x+["%0.6f" % y] for x,y in zip(results,qvals)]

print '\t'.join("ORF gene start end strand TAs sum mean sat NZsites NZmean expec_sum FC LFC pval qval".split())

for vals in results:

  print '\t'.join([str(x) for x in vals])

__all__ = ["transit_tools", "tnseq_tools", "norm_tools", "stat_tools"]

__version__ = "v3.0.1"

__version__ = "v3.0.2"

prefix = "[TRANSIT]"

        (args, kwargs) = transit_tools.cleanargs(rawargs)

        # ctrl-vs-exp = condition 1-vs-2

        # originally, MAD defined order of CL args this way: strA/cond1, strB/cond1, strA/cond2, strB/cond

        #ctrldataA = args[0].split(",")

        #ctrldataB = args[1].split(",")

        #expdataA = args[2].split(",")

        #expdataB = args[3].split(",")

        # TRI changed order of args this way: strA/cond1, strA/cond2, strB/cond1, strB/cond

        ctrldataA = args[0].split(",")

        ctrldataB = args[1].split(",")

        expdataA = args[2].split(",")

        expdataA = args[1].split(",")

        ctrldataB = args[2].split(",")

        expdataB = args[3].split(",")

        annotationPath = args[4]

        output_path = args[5]

        output_file = open(output_path, "w")

        Na2 = len(self.expdataA)

        Nb2 = len(self.expdataB)

        # Get data

        self.transit_message("Getting Data")

        (data, position) = transit_tools.get_validated_data(wiglist, wxobj=self.wxobj)

    @classmethod

    def usage_string(self):

        return """python %s GI <comma-separated .wig control files condition A> <comma-separated .wig control files condition B> <comma-separated .wig experimental files condition A> <comma-separated .wig experimental files condition B> <annotation .prot_table or GFF3> <output file> [Optional Arguments]

        #return """python %s GI <comma-separated .wig control files condition A> <comma-separated .wig control files condition B> <comma-separated .wig experimental files condition A> <comma-separated .wig experimental files condition B> <annotation .prot_table or GFF3> <output file> [Optional Arguments]

        return """python %s GI <wigs_for_strA_cond1> <wigs_for_strA_cond2> <wigs_for_strB_cond1> <wigs_for_strB_cond2> <annotation .prot_table or GFF3> <output file> [Optional Arguments]

        GI performs a comparison among 4 groups of datasets, strain A and B assessed in conditions 1 and 2 (e.g. control vs treatment).

        It looks for interactions where the response to the treatment (i.e. effect on insertion counts) depends on the strain (output variable: delta_LFC).

        Provide replicates in each group as a comma-separated list of wig files.

        Significant interactions are those with "HDI outside ROPE?"=TRUE, and all genes are sorted by significance using BFDR."

        Optional Arguments:

        -s <integer>    :=  Number of samples. Default: -s 10000

        -pc             :=  Pseudocounts to be added at each site.

        -l              :=  Perform LOESS Correction; Helps remove possible genomic position bias.

                            Default: Turned Off.

        -iN <float>     :=  Ignore TAs occuring within given percentage (as integer) of the N terminus. Default: -iN 0

        -iC <float>     :=  Ignore TAs occuring within given percentage (as integer) of the C terminus. Default: -iC 0

        -iN <int>       :=  Ignore TAs occuring within given percentage (as integer) of the N terminus. Default: -iN 0

        -iC <int>       :=  Ignore TAs occuring within given percentage (as integer) of the C terminus. Default: -iC 0

        --ctrl_lib      :=  String of letters representing library of control files in order

                            e.g. 'AABB'. Default empty. Letters used must also be used in --exp_lib

                            If non-empty, resampling will limit permutations to within-libraries.