<s> '''<STR_LIT>''' <EOL> import math <EOL> import numpy <EOL> import sys <EOL> import collections <EOL> import itertools <EOL> import re <EOL> import pandas <EOL> import ggplot <EOL> import copy <EOL> import numpy as np <EOL> from scipy . stats import ttest_ind <EOL> import matplotlib <EOL> import matplotlib . pyplot as plt <EOL> import rpy2 <EOL> from rpy2 . robjects import r as R <EOL> import rpy2 . robjects as ro <EOL> from rpy2 . robjects import pandas2ri <EOL> from rpy2 . robjects . packages import importr <EOL> from rpy2 . robjects . vectors import FloatVector <EOL> import os <EOL> try : <EOL> import CGAT . Experiment as E <EOL> import CGAT . IOTools as IOTools <EOL> import CGAT . Stats as Stats <EOL> except ImportError : <EOL> import Experiment as E <EOL> import IOTools <EOL> import Stats <EOL> pandas2ri . activate ( ) <EOL> def runDETest ( raw_DataFrame , <EOL> design_file , <EOL> outfile , <EOL> de_caller , <EOL> ** kwargs ) : <EOL> '''<STR_LIT>''' <EOL> if de_caller . lower ( ) == "<STR_LIT>" : <EOL> pass <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> def splitModel ( model ) : <EOL> '''<STR_LIT>''' <EOL> return [ x for x in <EOL> re . split ( "<STR_LIT>" , re . sub ( "<STR_LIT>" , "<STR_LIT>" , model ) ) if <EOL> len ( x ) > <NUM_LIT:0> ] <EOL> def adjustPvalues ( p_values ) : <EOL> '''<STR_LIT>''' <EOL> stats = importr ( '<STR_LIT>' ) <EOL> adj_pvalues = list ( stats . p_adjust ( FloatVector ( p_values ) , method = '<STR_LIT>' ) ) <EOL> return adj_pvalues <EOL> def pvaluesToSignficant ( p_values , fdr ) : <EOL> '''<STR_LIT>''' <EOL> return [ int ( x < fdr ) for x in p_values ] <EOL> def makeMAPlot ( resultsTable , title , outfile ) : <EOL> '''<STR_LIT>''' <EOL> plotter = R ( '''<STR_LIT>''' ) <EOL> r_resultsTable = pandas2ri . py2ri ( resultsTable ) <EOL> plotter ( r_resultsTable , title , outfile ) <EOL> class ExperimentalDesign ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , filename_or_table ) : <EOL> if isinstance ( filename_or_table , str ) : <EOL> self . table = pandas . read_csv ( filename_or_table , sep = "<STR_LIT:\t>" , <EOL> index_col = <NUM_LIT:0> ) <EOL> else : <EOL> self . table = filename_or_table <EOL> assert self . table . shape , "<STR_LIT>" <EOL> self . _update ( ) <EOL> def _update ( self ) : <EOL> """<STR_LIT>""" <EOL> self . table = self . table [ self . table [ "<STR_LIT>" ] != <NUM_LIT:0> ] <EOL> self . conditions = self . table [ '<STR_LIT>' ] . tolist ( ) <EOL> self . pairs = self . table [ '<STR_LIT>' ] . tolist ( ) <EOL> self . groups = ( list ( collections . OrderedDict . fromkeys ( <EOL> self . conditions ) ) ) <EOL> self . samples = self . table . index . tolist ( ) <EOL> if len ( self . table . columns ) > <NUM_LIT:3> : <EOL> self . factors = self . table . iloc [ : , <NUM_LIT:3> : ] <EOL> else : <EOL> self . factors = None <EOL> max_per_group = max ( [ self . conditions . count ( x ) for x in self . groups ] ) <EOL> self . has_replicates = max_per_group >= <NUM_LIT:2> <EOL> npairs = len ( set ( self . pairs ) ) <EOL> has_pairs = npairs == <NUM_LIT:2> <EOL> if has_pairs : <EOL> min_per_pair = min ( [ self . pairs . count ( x ) for x in set ( self . pairs ) ] ) <EOL> self . has_pairs = min_per_pair >= <NUM_LIT:2> <EOL> else : <EOL> self . has_pairs = False <EOL> def validate ( self , counts , model = None ) : <EOL> missing = set ( self . samples ) . difference ( set ( counts . table . columns ) ) <EOL> if len ( missing ) > <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % "<STR_LIT:U+002CU+0020>" . join ( missing ) ) <EOL> if model is not None : <EOL> model_terms = splitModel ( model ) <EOL> missing = set ( model_terms ) . difference ( <EOL> set ( self . table . columns . tolist ( ) ) ) <EOL> if len ( missing ) > <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" % <EOL> "<STR_LIT:U+002CU+0020>" . join ( missing ) ) <EOL> for term in model_terms : <EOL> levels = set ( self . factors . ix [ : , term ] ) <EOL> if len ( levels ) < <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" % <EOL> ( term , "<STR_LIT:U+002CU+0020>" . join ( levels ) ) ) <EOL> def restrict ( self , counts ) : <EOL> '''<STR_LIT>''' <EOL> self . table = self . table . ix [ counts . table . columns , : ] <EOL> def revalidate ( self , counts , model = None ) : <EOL> '''<STR_LIT>''' <EOL> if len ( set ( self . samples ) . symmetric_difference ( <EOL> set ( counts . table . columns ) ) ) > <NUM_LIT:0> : <EOL> self . restrict ( counts ) <EOL> self . _update ( ) <EOL> self . validate ( counts , model ) <EOL> else : <EOL> pass <EOL> def firstPairOnly ( self ) : <EOL> '''<STR_LIT>''' <EOL> if not self . pairs : <EOL> self . pairs = self . table [ '<STR_LIT>' ] . tolist ( ) <EOL> self . table = self . table . ix [ self . table [ '<STR_LIT>' ] == min ( self . pairs ) , ] <EOL> def getSamplesInGroup ( self , group ) : <EOL> """<STR_LIT>""" <EOL> if group not in self . groups : <EOL> raise KeyError ( "<STR_LIT>" ) <EOL> return self . table [ self . table [ "<STR_LIT>" ] == group ] . index . tolist ( ) <EOL> def getGroupForSample ( self , sample ) : <EOL> """<STR_LIT>""" <EOL> return self . table . loc [ sample ] [ "<STR_LIT>" ] <EOL> def getGroups2Samples ( self ) : <EOL> """<STR_LIT>""" <EOL> groups_to_tracks = { } <EOL> for group in self . groups : <EOL> match_group = ( self . table [ '<STR_LIT>' ] == group ) . tolist ( ) <EOL> subset = self . table . iloc [ match_group , ] <EOL> groups_to_tracks [ group ] = subset . index . tolist ( ) <EOL> return groups_to_tracks <EOL> def mapGroupsSuffix ( self , shuffle_suffix , keep_suffix ) : <EOL> '''<STR_LIT>''' <EOL> groups_to_keep_tracks = { } <EOL> groups_to_spike_tracks = { } <EOL> keep_suffix = keep_suffix . split ( "<STR_LIT:U+002C>" ) <EOL> for group in self . groups : <EOL> match_group = ( self . table [ '<STR_LIT>' ] == group ) . tolist ( ) <EOL> tmp_design = self . table . iloc [ match_group , ] <EOL> groups_to_spike_tracks [ group ] = [ <EOL> x + shuffle_suffix for x in tmp_design . index . tolist ( ) ] <EOL> groups_to_keep_tracks [ group ] = copy . copy ( <EOL> groups_to_spike_tracks [ group ] ) <EOL> groups_to_keep_tracks [ group ] . extend ( <EOL> [ x + y for x in tmp_design . index . tolist ( ) for y in keep_suffix ] ) <EOL> return groups_to_keep_tracks , groups_to_spike_tracks <EOL> class DEResult ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self , testTable = None ) : <EOL> self . table = testTable <EOL> def getResults ( self ) : <EOL> '''<STR_LIT>''' <EOL> pass <EOL> def summariseDEResults ( self ) : <EOL> '''<STR_LIT>''' <EOL> self . Summary = { } <EOL> control_names = set ( self . table [ '<STR_LIT>' ] ) <EOL> treatment_names = set ( self . table [ '<STR_LIT>' ] ) <EOL> for control , treatment in itertools . product ( control_names , <EOL> treatment_names ) : <EOL> tmp_table = self . table [ self . table [ '<STR_LIT>' ] == control ] <EOL> tmp_table = tmp_table [ tmp_table [ '<STR_LIT>' ] == treatment ] <EOL> tmp_table . reset_index ( inplace = True ) <EOL> n_rows = tmp_table . shape [ <NUM_LIT:0> ] <EOL> if n_rows > <NUM_LIT:0> : <EOL> if control != treatment : <EOL> label = control + "<STR_LIT:_>" + treatment <EOL> else : <EOL> label = control <EOL> label = re . sub ( "<STR_LIT::>" , "<STR_LIT>" , label ) <EOL> counts = E . Counter ( ) <EOL> counts . signficant = sum ( tmp_table [ '<STR_LIT>' ] ) <EOL> counts . insignficant = ( len ( tmp_table [ '<STR_LIT>' ] ) - <EOL> counts . signficant ) <EOL> counts . all_over = sum ( [ x > <NUM_LIT:0> for x in tmp_table [ '<STR_LIT>' ] ] ) <EOL> counts . all_under = sum ( [ x < <NUM_LIT:0> for x in tmp_table [ '<STR_LIT>' ] ] ) <EOL> counts . signficant_over = sum ( <EOL> [ tmp_table [ '<STR_LIT>' ] [ x ] == <NUM_LIT:1> and <EOL> tmp_table [ '<STR_LIT>' ] [ x ] > <NUM_LIT:1> for x in range ( <NUM_LIT:0> , n_rows ) ] ) <EOL> counts . signficant_under = sum ( <EOL> [ tmp_table [ '<STR_LIT>' ] [ x ] == <NUM_LIT:1> and <EOL> tmp_table [ '<STR_LIT>' ] [ x ] < <NUM_LIT:1> for x in range ( <NUM_LIT:0> , n_rows ) ] ) <EOL> self . Summary [ label ] = counts <EOL> def plotMA ( self , contrast = None , outfile_prefix = None , <EOL> point_alpha = <NUM_LIT:1> , point_size = <NUM_LIT:1> ) : <EOL> '''<STR_LIT>''' <EOL> makeMAPlot = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> makeMAPlot ( pandas2ri . py2ri ( self . table ) ) <EOL> def plotVolcano ( self , contrast = None , outfile_prefix = None ) : <EOL> '''<STR_LIT>''' <EOL> makeVolcanoPlot = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> makeVolcanoPlot ( pandas2ri . py2ri ( self . table ) ) <EOL> class DEExperiment ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __init__ ( self ) : <EOL> pass <EOL> def __call__ ( self ) : <EOL> '''<STR_LIT>''' <EOL> self . callDifferentialExpression ( ) <EOL> def run ( self ) : <EOL> '''<STR_LIT>''' <EOL> class DEExperiment_TTest ( DEExperiment ) : <EOL> '''<STR_LIT>''' <EOL> def run ( self , counts , design , normalise = True , <EOL> normalise_method = "<STR_LIT>" ) : <EOL> if normalise is True : <EOL> counts . normalise ( method = normalise_method ) <EOL> df_dict = collections . defaultdict ( list ) <EOL> for combination in itertools . combinations ( design . groups , <NUM_LIT:2> ) : <EOL> control , treatment = combination <EOL> n_rows = counts . table . shape [ <NUM_LIT:0> ] <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( control , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( treatment , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( counts . table . index . tolist ( ) ) <EOL> df_dict [ "<STR_LIT:status>" ] . extend ( ( "<STR_LIT:OK>" , ) * n_rows ) <EOL> c_keep = [ x == control for <EOL> x in design . conditions ] <EOL> control_counts = counts . table . iloc [ : , c_keep ] <EOL> t_keep = [ x == treatment for <EOL> x in design . conditions ] <EOL> treatment_counts = counts . table . iloc [ : , t_keep ] <EOL> c_mean = control_counts . mean ( axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( c_mean ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( control_counts . std ( axis = <NUM_LIT:1> ) ) <EOL> t_mean = treatment_counts . mean ( axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( t_mean ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( treatment_counts . std ( axis = <NUM_LIT:1> ) ) <EOL> t , prob = ttest_ind ( control_counts , treatment_counts , axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( prob ) <EOL> result = DEResult_TTest ( testTable = pandas . DataFrame ( df_dict ) ) <EOL> result . table . set_index ( "<STR_LIT>" , inplace = True ) <EOL> return result <EOL> class DEResult_TTest ( DEResult ) : <EOL> def getResults ( self , fdr ) : <EOL> '''<STR_LIT>''' <EOL> self . table [ "<STR_LIT>" ] = ( <EOL> self . table [ "<STR_LIT>" ] / self . table [ "<STR_LIT>" ] ) <EOL> self . table [ "<STR_LIT>" ] = adjustPvalues ( self . table [ "<STR_LIT>" ] ) <EOL> self . table [ "<STR_LIT>" ] = pvaluesToSignficant ( <EOL> self . table [ "<STR_LIT>" ] , fdr ) <EOL> self . table [ "<STR_LIT>" ] = list ( numpy . log2 ( self . table [ "<STR_LIT>" ] ) ) <EOL> self . table [ "<STR_LIT>" ] = self . table [ "<STR_LIT>" ] <EOL> self . table [ "<STR_LIT>" ] = "<STR_LIT>" . join ( ( self . table [ '<STR_LIT>' ] , <EOL> self . table [ '<STR_LIT>' ] ) ) <EOL> class DEExperiment_edgeR ( DEExperiment ) : <EOL> '''<STR_LIT>''' <EOL> def run ( self , <EOL> counts , <EOL> design , <EOL> model = None , <EOL> dispersion = None , <EOL> ref_group = None , <EOL> contrasts = None , <EOL> outfile_prefix = None ) : <EOL> if not design . has_replicates and dispersion is None : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> r_counts = pandas2ri . py2ri ( counts . table ) <EOL> r_groups = ro . StrVector ( design . conditions ) <EOL> r_pairs = ro . StrVector ( design . pairs ) <EOL> r_has_pairs = ro . default_py2ri ( design . has_pairs ) <EOL> r_has_replicates = ro . default_py2ri ( design . has_replicates ) <EOL> if dispersion is not None : <EOL> r_dispersion = ro . default_py2ri ( dispersion ) <EOL> else : <EOL> r_dispersion = ro . default_py2ri ( False ) <EOL> if model is not None : <EOL> r_factors_df = pandas2ri . py2ri ( design . factors ) <EOL> else : <EOL> r_factors_df = ro . default_py2ri ( False ) <EOL> if ref_group is not None : <EOL> r_ref_group = ro . default_py2ri ( ref_group ) <EOL> else : <EOL> r_ref_group = ro . default_py2ri ( design . groups [ <NUM_LIT:0> ] ) <EOL> if contrasts is not None : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> r_contrasts = ro . default_py2ri ( contrasts ) <EOL> E . info ( '<STR_LIT>' % <EOL> ( design . groups , design . pairs , design . has_replicates , <EOL> design . has_pairs ) ) <EOL> buildDGEList = R ( '''<STR_LIT>''' ) <EOL> r_countsTable = buildDGEList ( r_counts , r_groups , <EOL> r_ref_group , r_factors_df ) <EOL> if model is None : <EOL> buildDesign = R ( '''<STR_LIT>''' ) <EOL> r_design = buildDesign ( r_countsTable , r_has_pairs ) <EOL> else : <EOL> buildDesign = R ( '''<STR_LIT>''' % model ) <EOL> r_design = buildDesign ( r_factors_df ) <EOL> E . info ( "<STR_LIT>" ) <EOL> E . info ( r_design ) <EOL> fitModel = R ( '''<STR_LIT>''' ) <EOL> r_fit = fitModel ( r_countsTable , r_design , <EOL> r_has_replicates , r_dispersion ) <EOL> E . info ( "<STR_LIT>" ) <EOL> if contrasts is None : <EOL> lrtTest = R ( '''<STR_LIT>''' ) <EOL> r_lrt_table = lrtTest ( r_fit , outfile_prefix , <EOL> r_countsTable , r_design ) <EOL> else : <EOL> pass <EOL> E . info ( "<STR_LIT>" ) <EOL> outputCPMTable = R ( '''<STR_LIT>''' ) <EOL> outputCPMTable ( r_countsTable , outfile_prefix ) <EOL> result = DEResult_edgeR ( testTable = pandas2ri . ri2py ( r_lrt_table ) ) <EOL> return result <EOL> class DEResult_edgeR ( DEResult ) : <EOL> def getResults ( self , fdr , DEtype = "<STR_LIT>" ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" ) <EOL> df_dict = collections . defaultdict ( ) <EOL> n_rows = self . table . shape [ <NUM_LIT:0> ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = adjustPvalues ( self . table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] = pvaluesToSignficant ( <EOL> df_dict [ "<STR_LIT>" ] , fdr ) <EOL> df_dict [ "<STR_LIT>" ] = list ( numpy . log2 ( self . table [ '<STR_LIT>' ] ) ) <EOL> df_dict [ "<STR_LIT>" ] = df_dict [ "<STR_LIT>" ] <EOL> df_dict [ "<STR_LIT>" ] = [ math . pow ( <NUM_LIT:2> , float ( x ) ) for <EOL> x in self . table [ '<STR_LIT>' ] ] <EOL> df_dict [ "<STR_LIT:status>" ] = ( "<STR_LIT:OK>" , ) * n_rows <EOL> self . table = pandas . DataFrame ( df_dict ) <EOL> self . table . set_index ( "<STR_LIT>" , inplace = True ) <EOL> def plotMAplot ( self , design , outfile_prefix ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> class DEExperiment_DESeq ( DEExperiment ) : <EOL> '''<STR_LIT>''' <EOL> def run ( self , <EOL> counts , <EOL> design , <EOL> model = None , <EOL> dispersion_method = "<STR_LIT>" , <EOL> ref_group = None , <EOL> contrasts = None , <EOL> outfile_prefix = None , <EOL> sharing_mode = "<STR_LIT>" , <EOL> fit_type = "<STR_LIT>" ) : <EOL> r_counts = pandas2ri . py2ri ( counts . table ) <EOL> r_groups = ro . StrVector ( design . conditions ) <EOL> r_pairs = ro . StrVector ( design . pairs ) <EOL> r_has_pairs = ro . default_py2ri ( design . has_pairs ) <EOL> r_has_replicates = ro . default_py2ri ( design . has_replicates ) <EOL> r_dispersion_method = ro . default_py2ri ( dispersion_method ) <EOL> r_sharing_mode = ro . default_py2ri ( sharing_mode ) <EOL> r_fit_type = ro . default_py2ri ( fit_type ) <EOL> if model is not None : <EOL> r_factors_df = pandas2ri . py2ri ( design . factors ) <EOL> else : <EOL> r_factors_df = ro . default_py2ri ( False ) <EOL> if ref_group is not None : <EOL> r_ref_group = ro . default_py2ri ( ref_group ) <EOL> else : <EOL> r_ref_group = ro . default_py2ri ( design . groups [ <NUM_LIT:0> ] ) <EOL> if contrasts is not None : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> r_contrasts = ro . default_py2ri ( contrasts ) <EOL> if dispersion_method == "<STR_LIT>" and not design . has_replicates : <EOL> E . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> r_dispersion_method = ro . default_py2ri ( "<STR_LIT>" ) <EOL> E . info ( '<STR_LIT>' % <EOL> ( design . groups , design . pairs , design . has_replicates , <EOL> design . has_pairs ) ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( dispersion_method , fit_type , sharing_mode ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> buildCountDataSet = R ( '''<STR_LIT>''' ) <EOL> buildCountDataSet ( r_counts , r_groups , r_dispersion_method , <EOL> r_fit_type , r_sharing_mode ) <EOL> result = None <EOL> return result <EOL> class DEResult_DEResult ( DEResult ) : <EOL> def getResults ( self , fdr , DEtype = "<STR_LIT>" ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" ) <EOL> df_dict = collections . defaultdict ( ) <EOL> n_rows = self . table . shape [ <NUM_LIT:0> ] <EOL> if DEtype == "<STR_LIT>" : <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> else : <EOL> pass <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = adjustPvalues ( self . table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] = pvaluesToSignficant ( <EOL> df_dict [ "<STR_LIT>" ] , fdr ) <EOL> df_dict [ "<STR_LIT>" ] = list ( numpy . log2 ( self . table [ '<STR_LIT>' ] ) ) <EOL> df_dict [ "<STR_LIT>" ] = df_dict [ "<STR_LIT>" ] <EOL> df_dict [ "<STR_LIT>" ] = [ math . pow ( <NUM_LIT:2> , float ( x ) ) for <EOL> x in self . table [ '<STR_LIT>' ] ] <EOL> df_dict [ "<STR_LIT:status>" ] = ( "<STR_LIT:OK>" , ) * n_rows <EOL> self . table = pandas . DataFrame ( df_dict ) <EOL> self . table . set_index ( "<STR_LIT>" , inplace = True ) <EOL> def plotMAplot ( self , design , outfile_prefix ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> class DEExperiment_DESeq2 ( DEExperiment ) : <EOL> '''<STR_LIT>''' <EOL> def run ( self , <EOL> counts , <EOL> design , <EOL> model = None , <EOL> contrasts = None , <EOL> outfile_prefix = None , <EOL> fdr = <NUM_LIT:0.1> ) : <EOL> counts . table = counts . table [ design . table . index ] <EOL> r_counts = pandas2ri . py2ri ( counts . table ) <EOL> r_groups = ro . StrVector ( design . conditions ) <EOL> r_pairs = ro . StrVector ( design . pairs ) <EOL> r_has_pairs = ro . default_py2ri ( design . has_pairs ) <EOL> r_has_replicates = ro . default_py2ri ( design . has_replicates ) <EOL> if design . factors is not None : <EOL> r_factors_df = pandas2ri . py2ri ( design . factors ) <EOL> else : <EOL> r_factors_df = ro . default_py2ri ( False ) <EOL> if contrasts is not None : <EOL> DEtype = "<STR_LIT>" <EOL> if not model : <EOL> if design . factors is not None : <EOL> model = "<STR_LIT>" + "<STR_LIT:+>" . join ( design . factors . columns ) <EOL> model_terms = design . factors . columns . values . tolist ( ) <EOL> else : <EOL> E . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> else : <EOL> model_terms = [ x for x in re . split ( "<STR_LIT>" , model ) [ <NUM_LIT:1> : ] <EOL> if x != "<STR_LIT:0>" ] <EOL> r_model = ro . default_py2ri ( model ) <EOL> else : <EOL> DEtype = "<STR_LIT>" <EOL> r_DEtype = ro . default_py2ri ( DEtype ) <EOL> E . info ( '<STR_LIT>' <EOL> '<STR_LIT>' % <EOL> ( design . groups , design . pairs , design . has_replicates , <EOL> design . has_pairs ) ) <EOL> E . info ( design . factors ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if DEtype == "<STR_LIT>" : <EOL> buildDesign = R ( '''<STR_LIT>''' ) <EOL> r_design = buildDesign ( r_counts , r_groups ) <EOL> buildCountDataSet = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> r_dds = buildCountDataSet ( r_counts , r_design ) <EOL> performDifferentialTesting = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> results = pandas2ri . ri2py ( performDifferentialTesting ( r_dds ) ) <EOL> else : <EOL> r_design = pandas2ri . py2ri ( design . table ) <EOL> buildCountDataSet = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> r_dds = buildCountDataSet ( r_counts , r_design , r_model ) <EOL> results = pandas . DataFrame ( ) <EOL> n = <NUM_LIT:0> <EOL> for contrast in contrasts : <EOL> assert contrast in design . table . columns , ( <EOL> "<STR_LIT>" % contrast ) <EOL> model = [ x for x in model_terms if x != contrast ] <EOL> model = "<STR_LIT>" + "<STR_LIT:+>" . join ( model ) <EOL> performDifferentialTesting = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> tmp_results = pandas2ri . ri2py ( performDifferentialTesting ( r_dds ) ) <EOL> n2 = n + tmp_results . shape [ <NUM_LIT:0> ] <EOL> tmp_results . index = range ( n , n2 ) <EOL> n = n2 <EOL> results = results . append ( tmp_results ) <EOL> final_result = DEResult_DESeq2 ( testTable = results ) <EOL> return final_result <EOL> class DEResult_DESeq2 ( DEResult ) : <EOL> def getResults ( self , fdr ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" ) <EOL> df_dict = collections . defaultdict ( ) <EOL> n_rows = self . table . shape [ <NUM_LIT:0> ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = adjustPvalues ( self . table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] = pvaluesToSignficant ( <EOL> df_dict [ "<STR_LIT>" ] , fdr ) <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . table [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = [ math . pow ( <NUM_LIT:2> , float ( x ) ) for <EOL> x in self . table [ '<STR_LIT>' ] ] <EOL> df_dict [ "<STR_LIT:status>" ] = ( "<STR_LIT:OK>" , ) * n_rows <EOL> self . table = pandas . DataFrame ( df_dict ) <EOL> class DEExperiment_Sleuth ( DEExperiment ) : <EOL> '''<STR_LIT>''' <EOL> def run ( self , <EOL> design , <EOL> base_dir , <EOL> model = None , <EOL> contrast = None , <EOL> outfile_prefix = None , <EOL> counts = None , <EOL> tpm = None , <EOL> fdr = <NUM_LIT:0.1> , <EOL> lrt = False , <EOL> reduced_model = None ) : <EOL> if lrt : <EOL> E . info ( "<STR_LIT>" ) <EOL> assert reduced_model is not None , ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> design . table [ '<STR_LIT>' ] = design . table . index <EOL> r_design_df = pandas2ri . py2ri ( design . table ) <EOL> E . info ( '<STR_LIT>' <EOL> '<STR_LIT>' % <EOL> ( design . groups , design . pairs , design . has_replicates , <EOL> design . has_pairs ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> model_terms = [ x for x in re . split ( "<STR_LIT>" , model ) [ <NUM_LIT:1> : ] <EOL> if x != "<STR_LIT:0>" ] <EOL> variates = "<STR_LIT>" % "<STR_LIT:U+002C>" . join ( model_terms ) <EOL> createSleuthObject = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> so = createSleuthObject ( r_design_df ) <EOL> if counts : <EOL> makeCountsTable = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> makeCountsTable ( so ) <EOL> if tpm : <EOL> makeTPMTable = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> makeTPMTable ( so ) <EOL> if lrt : <EOL> differentialTesting = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> else : <EOL> differentialTesting = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> results = differentialTesting ( so ) <EOL> final_result = DEResult_Sleuth ( so = results , contrast = contrast , lrt = lrt ) <EOL> return final_result <EOL> class DEResult_Sleuth ( DEResult ) : <EOL> def __init__ ( self , so = None , contrast = None , lrt = False ) : <EOL> '''<STR_LIT>''' <EOL> self . contrast = contrast <EOL> self . lrt = lrt <EOL> self . so = so <EOL> R ( '''<STR_LIT>''' ) <EOL> if self . lrt : <EOL> extractSleuthResults = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> else : <EOL> extractSleuthResults = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> self . sleuthResults = pandas2ri . ri2py ( extractSleuthResults ( self . so ) ) <EOL> def getResults ( self , fdr ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" ) <EOL> df_dict = collections . defaultdict ( ) <EOL> n_rows = self . sleuthResults . shape [ <NUM_LIT:0> ] <EOL> df_dict [ "<STR_LIT>" ] = ( "<STR_LIT>" , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = ( "<STR_LIT>" , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = self . sleuthResults [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = self . contrast <EOL> df_dict [ "<STR_LIT>" ] = [ math . exp ( float ( x ) ) for <EOL> x in self . sleuthResults [ '<STR_LIT>' ] ] <EOL> df_dict [ "<STR_LIT>" ] = df_dict [ "<STR_LIT>" ] <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = ( <NUM_LIT:0> , ) * n_rows <EOL> df_dict [ "<STR_LIT>" ] = self . sleuthResults [ '<STR_LIT>' ] <EOL> df_dict [ "<STR_LIT>" ] = adjustPvalues ( self . sleuthResults [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] = pvaluesToSignficant ( df_dict [ "<STR_LIT>" ] , <EOL> fdr ) <EOL> df_dict [ "<STR_LIT>" ] = [ math . exp ( float ( x ) ) for <EOL> x in self . sleuthResults [ '<STR_LIT:b>' ] ] <EOL> df_dict [ "<STR_LIT>" ] = [ math . log ( float ( x ) , <NUM_LIT:2> ) for x in df_dict [ '<STR_LIT>' ] ] <EOL> df_dict [ "<STR_LIT>" ] = df_dict [ "<STR_LIT>" ] <EOL> df_dict [ "<STR_LIT:status>" ] = ( "<STR_LIT:OK>" , ) * n_rows <EOL> self . table = pandas . DataFrame ( df_dict ) <EOL> def buildProbeset2Gene ( infile , <EOL> outfile , <EOL> database = "<STR_LIT>" , <EOL> mapping = "<STR_LIT>" ) : <EOL> '''<STR_LIT>''' <EOL> R . library ( database ) <EOL> m = R ( mapping ) <EOL> result = R . toTable ( m ) <EOL> outf = open ( outfile , "<STR_LIT:w>" ) <EOL> outf . write ( "<STR_LIT>" ) <EOL> for probeset_id , gene_id in zip ( result [ "<STR_LIT>" ] , <EOL> result [ "<STR_LIT>" ] ) : <EOL> outf . write ( "<STR_LIT>" % ( probeset_id , gene_id ) ) <EOL> outf . close ( ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( len ( result ) , <EOL> outfile , <EOL> len ( set ( result [ "<STR_LIT>" ] ) ) , <EOL> len ( set ( result [ "<STR_LIT>" ] ) ) ) ) <EOL> GeneExpressionResult = collections . namedtuple ( <EOL> "<STR_LIT>" , <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> def writeExpressionResults ( outfile , result ) : <EOL> '''<STR_LIT>''' <EOL> if outfile == sys . stdout : <EOL> outf = outfile <EOL> else : <EOL> outf = IOTools . openFile ( outfile , "<STR_LIT:w>" ) <EOL> outf . write ( "<STR_LIT>" % "<STR_LIT:\t>" . join ( GeneExpressionResult . _fields ) ) <EOL> for x in sorted ( result ) : <EOL> outf . write ( "<STR_LIT>" % "<STR_LIT:\t>" . join ( map ( str , x ) ) ) <EOL> if outf != sys . stdout : <EOL> outf . close ( ) <EOL> class WelchsTTest ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __call__ ( self , <EOL> probesets , <EOL> treatments , <EOL> controls ) : <EOL> assert len ( probesets ) == len ( treatments [ <NUM_LIT:0> ] ) <EOL> assert len ( probesets ) == len ( controls [ <NUM_LIT:0> ] ) <EOL> nskipped = <NUM_LIT:0> <EOL> results = [ ] <EOL> for probeset , treatment , control in zip ( <EOL> probesets , zip ( * treatments ) , zip ( * controls ) ) : <EOL> nval1 , nval2 = len ( treatment ) , len ( control ) <EOL> mean1 , mean2 = numpy . mean ( treatment ) , numpy . mean ( control ) <EOL> stddev1 , stddev2 = numpy . std ( treatment ) , numpy . std ( control ) <EOL> try : <EOL> s = Stats . doWelchsTTest ( nval1 , mean1 , stddev1 , <EOL> nval2 , mean2 , stddev2 , <EOL> alpha = <NUM_LIT> ) <EOL> except ValueError : <EOL> E . warn ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" % probeset ) <EOL> nskipped += <NUM_LIT:1> <EOL> continue <EOL> s . mProbeset = probeset <EOL> results . append ( s ) <EOL> qvalues = Stats . doFDR ( [ x . mPValue for x in results ] ) . mQValues <EOL> for s , qvalue in zip ( results , qvalues ) : <EOL> s . mQValue = qvalue <EOL> return results , nskipped <EOL> class SAMR ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __call__ ( self , probesets , <EOL> treatments , <EOL> controls , <EOL> pattern = None , <EOL> fdr = <NUM_LIT> , <EOL> ngenes = None , <EOL> npermutations = <NUM_LIT:1000> , <EOL> ndelta = <NUM_LIT:10> , <EOL> method = "<STR_LIT>" ) : <EOL> if ngenes and fdr : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> R . library ( "<STR_LIT>" ) <EOL> m = numpy . matrix ( treatments + controls ) <EOL> m = numpy . transpose ( m ) <EOL> labels = numpy . array ( [ <NUM_LIT:1> ] * len ( treatments ) + [ <NUM_LIT:2> ] * len ( controls ) ) <EOL> R . assign ( "<STR_LIT:x>" , numpy . array ( m ) ) <EOL> R . assign ( "<STR_LIT:y>" , labels ) <EOL> R . assign ( "<STR_LIT>" , probesets ) <EOL> data = R ( <EOL> '''<STR_LIT>''' ) <EOL> result = R ( <EOL> '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> class SAM ( object ) : <EOL> '''<STR_LIT>''' <EOL> def __call__ ( self , probesets , <EOL> treatments , <EOL> controls , <EOL> pattern = None , <EOL> fdr = <NUM_LIT> , <EOL> ngenes = None , <EOL> npermutations = <NUM_LIT:1000> , <EOL> ndelta = <NUM_LIT:10> , <EOL> method = "<STR_LIT>" , <EOL> use_excel_sam = False , <EOL> treatment_label = "<STR_LIT>" , <EOL> control_label = "<STR_LIT>" ) : <EOL> if ngenes and fdr : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> R . library ( "<STR_LIT>" ) <EOL> m = numpy . matrix ( treatments + controls ) <EOL> m = numpy . transpose ( m ) <EOL> E . debug ( "<STR_LIT>" % m . shape ) <EOL> labels = numpy . array ( [ <NUM_LIT:1> ] * len ( treatments ) + [ <NUM_LIT:0> ] * len ( controls ) ) <EOL> kwargs = { } <EOL> if use_excel_sam : <EOL> kwargs . update ( <EOL> { "<STR_LIT>" : <EOL> R ( '''<STR_LIT>''' % <EOL> locals ( ) ) , <EOL> "<STR_LIT>" : True , <EOL> "<STR_LIT>" : True } ) <EOL> else : <EOL> kwargs . update ( { "<STR_LIT>" : <EOL> R ( '''<STR_LIT>''' % <EOL> locals ( ) ) } , ) <EOL> if method == "<STR_LIT>" : <EOL> kwargs [ "<STR_LIT>" ] = R ( '''<STR_LIT>''' ) <EOL> kwargs [ "<STR_LIT:B>" ] = npermutations <EOL> elif method == "<STR_LIT>" : <EOL> kwargs [ "<STR_LIT>" ] = R ( '''<STR_LIT>''' ) <EOL> elif method == "<STR_LIT>" : <EOL> kwargs [ "<STR_LIT>" ] = R ( '''<STR_LIT>''' ) <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" % method ) <EOL> E . info ( "<STR_LIT>" % str ( kwargs ) ) <EOL> a = R . sam ( numpy . array ( m ) , <EOL> labels , <EOL> gene_names = numpy . array ( probesets ) , <EOL> ** kwargs ) <EOL> R . assign ( "<STR_LIT:a>" , a ) <EOL> fdr_data = collections . namedtuple ( "<STR_LIT>" , ( <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT:false>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> cutoff_data = collections . namedtuple ( <EOL> "<STR_LIT>" , ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> gene_data = collections . namedtuple ( <EOL> "<STR_LIT>" , ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> def _totable ( robj ) : <EOL> '''<STR_LIT>''' <EOL> s = numpy . matrix ( robj ) <EOL> t = [ numpy . array ( x ) . reshape ( - <NUM_LIT:1> , ) for x in s ] <EOL> return t <EOL> t = _totable ( a . do_slot ( '<STR_LIT>' ) ) <EOL> assert len ( t [ <NUM_LIT:0> ] ) == len ( fdr_data . _fields ) <EOL> fdr_values = [ fdr_data ( * x ) for x in t ] <EOL> if fdr is not None and fdr > <NUM_LIT:0> : <EOL> s = numpy . matrix ( R . findDelta ( a , fdr ) ) <EOL> try : <EOL> cutoffs = [ cutoff_data ( * numpy . array ( x ) . reshape ( - <NUM_LIT:1> , ) ) <EOL> for x in s ] <EOL> E . debug ( "<STR_LIT>" % ( fdr , str ( cutoffs ) ) ) <EOL> cutoff = cutoffs [ - <NUM_LIT:1> ] <EOL> except TypeError : <EOL> E . debug ( "<STR_LIT>" ) <EOL> cutoff = None <EOL> elif ngenes : <EOL> s = numpy . matrix ( R . findDelta ( a , ngenes ) ) <EOL> try : <EOL> cutoffs = [ cutoff_data ( * numpy . array ( x ) . reshape ( - <NUM_LIT:1> , ) ) <EOL> for x in s ] <EOL> E . debug ( "<STR_LIT>" % ( fdr , str ( cutoffs ) ) ) <EOL> cutoff = cutoffs [ - <NUM_LIT:1> ] <EOL> except TypeError : <EOL> E . debug ( "<STR_LIT>" ) <EOL> cutoff = None <EOL> else : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> pvalues = dict ( zip ( probesets , R ( '''<STR_LIT>''' ) ) ) <EOL> qvalues = dict ( zip ( probesets , R ( '''<STR_LIT>''' ) ) ) <EOL> if pattern : <EOL> outfile = pattern % "<STR_LIT>" <EOL> R . pdf ( outfile ) <EOL> if cutoff : <EOL> R . plot ( a , cutoff . delta ) <EOL> else : <EOL> R . plot ( a ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> siggenes = { } <EOL> significant_genes = set ( ) <EOL> if cutoff is not None : <EOL> E . debug ( "<STR_LIT>" % str ( cutoff ) ) <EOL> summary = R ( '''<STR_LIT>''' % cutoff . delta ) <EOL> R . assign ( "<STR_LIT>" , summary ) <EOL> significant_genes = set ( <EOL> [ probesets [ int ( x ) - <NUM_LIT:1> ] for x <EOL> in R ( '''<STR_LIT>''' ) ] ) <EOL> r_result = zip ( * _totable ( summary . do_slot ( '<STR_LIT>' ) ) ) <EOL> if len ( r_result ) > <NUM_LIT:0> : <EOL> assert len ( r_result [ <NUM_LIT:0> ] ) == <NUM_LIT:6> , "<STR_LIT>" % len ( r_result [ <NUM_LIT:0> ] ) <EOL> for x in r_result : <EOL> if x [ <NUM_LIT:4> ] > fdr : <EOL> E . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( str ( x ) , x [ <NUM_LIT:4> ] ) ) <EOL> siggenes [ probesets [ int ( x [ <NUM_LIT:0> ] ) - <NUM_LIT:1> ] ] = gene_data ( * x ) <EOL> else : <EOL> E . debug ( "<STR_LIT>" ) <EOL> genes = [ ] <EOL> for probeset , treatment , control in zip ( <EOL> probesets , zip ( * treatments ) , zip ( * controls ) ) : <EOL> mean1 , mean2 = numpy . mean ( treatment ) , numpy . mean ( control ) <EOL> if probeset in siggenes : <EOL> s = siggenes [ probeset ] <EOL> pvalue = s . rawp <EOL> qvalue = s . qvalue <EOL> else : <EOL> pvalue = pvalues [ probeset ] <EOL> qvalue = qvalues [ probeset ] <EOL> significant = ( <NUM_LIT:0> , <NUM_LIT:1> ) [ probeset in significant_genes ] <EOL> genes . append ( GeneExpressionResult . _make ( ( probeset , <EOL> treatment_label , <EOL> mean1 , <EOL> numpy . std ( treatment ) , <EOL> control_label , <EOL> mean2 , <EOL> numpy . std ( control ) , <EOL> pvalue , <EOL> qvalue , <EOL> mean1 - mean2 , <EOL> math . pow ( <EOL> <NUM_LIT:2> , mean1 - mean2 ) , <EOL> math . pow ( <EOL> <NUM_LIT:2> , mean1 - mean2 ) , <EOL> significant , <EOL> "<STR_LIT:OK>" ) ) ) <EOL> return genes , cutoff , fdr_values <EOL> def loadTagData ( tags_filename , design_filename ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" % tags_filename ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> E . info ( "<STR_LIT>" % <EOL> tuple ( R ( '''<STR_LIT>''' ) ) ) <EOL> E . debug ( "<STR_LIT>" % R ( '''<STR_LIT>''' ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . debug ( "<STR_LIT>" % R ( '''<STR_LIT>''' ) ) <EOL> pheno = R ( <EOL> '''<STR_LIT>''' ) <EOL> missing = R ( '''<STR_LIT>''' ) <EOL> if missing : <EOL> E . warn ( "<STR_LIT>" % <EOL> missing ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . debug ( "<STR_LIT>" % <EOL> R ( '''<STR_LIT>''' ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" % <EOL> tuple ( R ( '''<STR_LIT>''' ) ) ) <EOL> def filterTagData ( filter_min_counts_per_row = <NUM_LIT:1> , <EOL> filter_min_counts_per_sample = <NUM_LIT:10> , <EOL> filter_percentile_rowsums = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_row ) <EOL> E . info ( "<STR_LIT>" % <EOL> tuple ( R ( '''<STR_LIT>''' ) ) ) <EOL> observations , samples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( observations , samples ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> empty_samples = tuple ( <EOL> R ( '''<STR_LIT>''' % filter_min_counts_per_sample ) ) <EOL> sample_names = R ( '''<STR_LIT>''' ) <EOL> nempty_samples = sum ( empty_samples ) <EOL> if nempty_samples : <EOL> E . warn ( "<STR_LIT>" % <EOL> ( nempty_samples , <EOL> "<STR_LIT:U+002C>" . join ( [ sample_names [ x ] <EOL> for x , y in enumerate ( empty_samples ) if y ] ) ) ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_sample ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_sample ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_sample ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_sample ) <EOL> observations , samples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> if filter_percentile_rowsums > <NUM_LIT:0> : <EOL> percentile = float ( filter_percentile_rowsums ) / <NUM_LIT> <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' % <EOL> locals ( ) ) <EOL> discard , keep = R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( filter_percentile_rowsums , <EOL> keep , discard ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> observations , samples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> return observations , samples <EOL> def groupTagData ( ref_group = None , ref_regex = None ) : <EOL> '''<STR_LIT>''' <EOL> if ref_regex is not None and ref_group is None : <EOL> groups = R ( '''<STR_LIT>''' ) <EOL> for g in groups : <EOL> rx = re . compile ( ref_regex ) <EOL> if rx . search ( g ) : <EOL> ref_group = g <EOL> if ref_group is not None : <EOL> E . info ( "<STR_LIT>" % ref_group ) <EOL> R ( '''<STR_LIT>''' % ref_group ) <EOL> groups = R ( '''<STR_LIT>''' ) <EOL> pairs = R ( '''<STR_LIT>''' ) <EOL> factors = R ( '''<STR_LIT>''' ) <EOL> if isinstance ( factors , rpy2 . robjects . vectors . DataFrame ) : <EOL> E . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" % factors . r_repr ( ) ) <EOL> max_per_group = R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] <EOL> has_replicates = max_per_group >= <NUM_LIT:2> <EOL> npairs = R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] <EOL> has_pairs = npairs == <NUM_LIT:2> <EOL> if has_pairs : <EOL> min_per_pair = R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] <EOL> has_pairs = min_per_pair >= <NUM_LIT:2> <EOL> return groups , pairs , has_replicates , has_pairs <EOL> def plotCorrelationHeatmap ( method = "<STR_LIT>" ) : <EOL> '''<STR_LIT>''' <EOL> if method == "<STR_LIT>" : <EOL> R ( '''<STR_LIT>''' ) <EOL> else : <EOL> R ( '''<STR_LIT>''' % <EOL> method ) <EOL> R ( '''<STR_LIT>''' ) <EOL> def plotPairs ( ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> try : <EOL> R ( '''<STR_LIT>''' ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> def plotPCA ( groups = True ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if groups is False : <EOL> R ( '''<STR_LIT>''' ) <EOL> nrows , nlevels = R ( '''<STR_LIT>''' ) <EOL> if nlevels > <NUM_LIT:1> : <EOL> R ( '''<STR_LIT>''' ) <EOL> if nlevels > <NUM_LIT:2> : <EOL> R ( '''<STR_LIT>''' ) <EOL> try : <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> def runEdgeR ( outfile , <EOL> outfile_prefix = "<STR_LIT>" , <EOL> fdr = <NUM_LIT:0.1> , <EOL> prefix = "<STR_LIT>" , <EOL> dispersion = None , <EOL> ref_group = None , <EOL> ref_regex = None , <EOL> ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> groups , pairs , has_replicates , has_pairs = groupTagData ( ref_group , <EOL> ref_regex ) <EOL> R . png ( '<STR_LIT>' % locals ( ) ) <EOL> plotCorrelationHeatmap ( ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> E . info ( '<STR_LIT>' % <EOL> ( groups , pairs , has_replicates , has_pairs ) ) <EOL> if has_pairs : <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> first = True <EOL> for g1 , g2 in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> R ( '''<STR_LIT>''' % g1 ) <EOL> R ( '''<STR_LIT>''' % g2 ) <EOL> if first : <EOL> R ( '''<STR_LIT>''' ) <EOL> first = False <EOL> else : <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> first = True <EOL> legend = [ ] <EOL> for pair in pairs : <EOL> for g1 , g2 in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> key = re . sub ( "<STR_LIT:->" , "<STR_LIT:_>" , "<STR_LIT>" % ( pair , g1 , g2 ) ) <EOL> legend . append ( key ) <EOL> R ( '''<STR_LIT>''' % ( <EOL> pair , g1 ) ) <EOL> R ( '''<STR_LIT>''' % ( <EOL> pair , g2 ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if first : <EOL> data = R ( '''<STR_LIT>''' % key ) <EOL> first = False <EOL> else : <EOL> R ( '''<STR_LIT>''' % key ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> try : <EOL> R . ggsave ( '''<STR_LIT>''' % locals ( ) ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if ref_group is not None : <EOL> R ( '''<STR_LIT>''' % ref_group ) <EOL> else : <EOL> R ( '''<STR_LIT>''' % groups [ <NUM_LIT:0> ] ) <EOL> E . info ( "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" ) <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> try : <EOL> R ( '''<STR_LIT>''' ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> if has_pairs : <EOL> R ( '''<STR_LIT>''' ) <EOL> else : <EOL> R ( '''<STR_LIT>''' ) <EOL> if has_replicates : <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> else : <EOL> if dispersion is None : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> E . warn ( "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' % <EOL> dispersion ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> rtype = collections . namedtuple ( "<STR_LIT>" , "<STR_LIT>" ) <EOL> if len ( groups ) == <NUM_LIT:2> : <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' % "<STR_LIT>" . join ( groups ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> results = [ ] <EOL> counts = E . Counter ( ) <EOL> df = R ( '''<STR_LIT>''' ) <EOL> df [ "<STR_LIT>" ] = R ( '''<STR_LIT>''' ) <EOL> counts . significant = sum ( df . padj <= fdr ) <EOL> counts . insignificant = sum ( df . padj > fdr ) <EOL> counts . significant_over = sum ( ( df . padj <= fdr ) & ( df . logFC > <NUM_LIT:0> ) ) <EOL> counts . significant_under = sum ( ( df . padj <= fdr ) & ( df . logFC < <NUM_LIT:0> ) ) <EOL> counts . input = len ( df ) <EOL> counts . all_over = sum ( df . logFC > <NUM_LIT:0> ) <EOL> counts . all_under = sum ( df . logFC < <NUM_LIT:0> ) <EOL> counts . fail = sum ( df . PValue . isnull ( ) ) <EOL> counts . ok = counts . input - counts . fail <EOL> df [ "<STR_LIT>" ] = df . logFC . pow ( <NUM_LIT> ) <EOL> df [ "<STR_LIT>" ] = df . padj <= fdr <EOL> for interval , d in df . iterrows ( ) : <EOL> results . append ( GeneExpressionResult . _make ( ( <EOL> interval , <EOL> groups [ <NUM_LIT:1> ] , <EOL> d . logCPM , <EOL> <NUM_LIT:0> , <EOL> groups [ <NUM_LIT:0> ] , <EOL> d . logCPM , <EOL> <NUM_LIT:0> , <EOL> d . PValue , <EOL> d . padj , <EOL> d . logFC , <EOL> d . fold , <EOL> d . logFC , <EOL> str ( int ( d . significant ) ) , <EOL> "<STR_LIT:OK>" ) ) ) <EOL> writeExpressionResults ( outfile , results ) <EOL> outf = IOTools . openFile ( "<STR_LIT>" % locals ( ) , "<STR_LIT:w>" ) <EOL> outf . write ( "<STR_LIT>" % counts . asTable ( ) ) <EOL> outf . close ( ) <EOL> def deseqPlotSizeFactors ( outfile ) : <EOL> '''<STR_LIT>''' <EOL> R . png ( outfile ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqOutputSizeFactors ( outfile ) : <EOL> '''<STR_LIT>''' <EOL> size_factors = R ( '''<STR_LIT>''' ) <EOL> samples = R ( '''<STR_LIT>''' ) <EOL> with IOTools . openFile ( outfile , "<STR_LIT:w>" ) as outf : <EOL> outf . write ( "<STR_LIT>" ) <EOL> for name , x in zip ( samples , size_factors ) : <EOL> outf . write ( "<STR_LIT>" % ( name , str ( x ) ) ) <EOL> def deseqPlotCorrelationHeatmap ( outfile , vsd ) : <EOL> '''<STR_LIT>''' <EOL> dists = R ( '''<STR_LIT>''' ) <EOL> R . png ( outfile ) <EOL> R [ '<STR_LIT>' ] ( <EOL> dists , <EOL> trace = '<STR_LIT:none>' , <EOL> margin = ro . IntVector ( ( <NUM_LIT:10> , <NUM_LIT:10> ) ) ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqPlotGeneHeatmap ( outfile , <EOL> data , <EOL> Rowv = False , <EOL> Colv = False ) : <EOL> '''<STR_LIT>''' <EOL> if len ( data ) == <NUM_LIT:0> : <EOL> return <EOL> if min ( R . dim ( data ) ) < <NUM_LIT:2> : <EOL> return <EOL> R . png ( outfile , width = <NUM_LIT> , height = <NUM_LIT> ) <EOL> hmcol = R . colorRampPalette ( R [ '<STR_LIT>' ] ( <NUM_LIT:9> , "<STR_LIT>" ) ) ( <NUM_LIT:100> ) <EOL> R [ '<STR_LIT>' ] ( <EOL> data , <EOL> col = hmcol , <EOL> trace = "<STR_LIT:none>" , <EOL> dendrogram = "<STR_LIT:none>" , <EOL> Rowv = Rowv , <EOL> Colv = Colv , <EOL> labRow = False , <EOL> margin = ro . IntVector ( ( <NUM_LIT:5> , <NUM_LIT:5> ) ) , <EOL> lhei = ro . IntVector ( ( <NUM_LIT:1> , <NUM_LIT:10> ) ) , <EOL> key = False ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqPlotPCA ( outfile , vsd , max_genes = <NUM_LIT> ) : <EOL> '''<STR_LIT>''' <EOL> R . png ( outfile ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> try : <EOL> R ( '''<STR_LIT>''' ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqPlotPairs ( outfile ) : <EOL> '''<STR_LIT>''' <EOL> R . png ( outfile , width = <NUM_LIT> , height = <NUM_LIT> ) <EOL> plotPairs ( ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqPlotPvaluesAgainstRowsums ( outfile ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> R . png ( outfile ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def deseqParseResults ( control_name , treatment_name , fdr , vsd = False ) : <EOL> '''<STR_LIT>''' <EOL> results = [ ] <EOL> isna = R [ "<STR_LIT>" ] <EOL> names = list ( R [ '<STR_LIT>' ] . names ) <EOL> m = dict ( [ ( x , x ) for x in names ] ) <EOL> m . update ( dict ( <EOL> pval = "<STR_LIT>" , <EOL> baseMeanA = "<STR_LIT>" , <EOL> baseMeanB = "<STR_LIT>" , <EOL> id = "<STR_LIT>" , <EOL> log2FoldChange = "<STR_LIT>" ) ) <EOL> rtype = collections . namedtuple ( "<STR_LIT>" , names ) <EOL> counts = E . Counter ( ) <EOL> for data in zip ( * R [ '<STR_LIT>' ] ) : <EOL> counts . input += <NUM_LIT:1> <EOL> d = rtype . _make ( data ) <EOL> if d . padj <= fdr : <EOL> signif = <NUM_LIT:1> <EOL> counts . significant += <NUM_LIT:1> <EOL> if d . log2FoldChange > <NUM_LIT:0> : <EOL> counts . significant_over += <NUM_LIT:1> <EOL> else : <EOL> counts . significant_under += <NUM_LIT:1> <EOL> else : <EOL> signif = <NUM_LIT:0> <EOL> counts . insignificant += <NUM_LIT:1> <EOL> if d . log2FoldChange > <NUM_LIT:0> : <EOL> counts . all_over += <NUM_LIT:1> <EOL> else : <EOL> counts . all_under += <NUM_LIT:1> <EOL> if d . baseMeanA == <NUM_LIT:0.0> and d . baseMeanB == <NUM_LIT:0.0> : <EOL> d = d . _replace ( foldChange = <NUM_LIT:0> , log2FoldChange = <NUM_LIT:0> ) <EOL> if d . pval != R ( '''<STR_LIT>''' ) : <EOL> status = "<STR_LIT:OK>" <EOL> else : <EOL> status = "<STR_LIT>" <EOL> counts [ status ] += <NUM_LIT:1> <EOL> counts . output += <NUM_LIT:1> <EOL> assert ( d . foldChange > <NUM_LIT:1> ) == ( d . baseMeanB > d . baseMeanA ) <EOL> results . append ( GeneExpressionResult . _make ( ( <EOL> d . id , <EOL> treatment_name , <EOL> d . baseMeanB , <EOL> <NUM_LIT:0> , <EOL> control_name , <EOL> d . baseMeanA , <EOL> <NUM_LIT:0> , <EOL> d . pval , <EOL> d . padj , <EOL> d . log2FoldChange , <EOL> d . foldChange , <EOL> d . transformed_log2FoldChange , <EOL> str ( signif ) , <EOL> status ) ) ) <EOL> return results , counts <EOL> def deseq2ParseResults ( fdr , vsd = False ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' % fdr ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> D = R ( '<STR_LIT>' ) <EOL> D . index = D [ '<STR_LIT>' ] <EOL> D = D . drop ( '<STR_LIT>' , <NUM_LIT:1> ) <EOL> return D <EOL> def runDESeq ( outfile , <EOL> outfile_prefix = "<STR_LIT>" , <EOL> fdr = <NUM_LIT:0.1> , <EOL> prefix = "<STR_LIT>" , <EOL> fit_type = "<STR_LIT>" , <EOL> dispersion_method = "<STR_LIT>" , <EOL> sharing_mode = "<STR_LIT>" , <EOL> ref_group = None , <EOL> ref_regex = None , <EOL> ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> groups , pairs , has_replicates , has_pairs = groupTagData ( ref_group , <EOL> ref_regex ) <EOL> E . info ( "<STR_LIT>" % ( has_replicates ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> no_size_factors = R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] <EOL> if no_size_factors : <EOL> E . warn ( "<STR_LIT>" ) <EOL> return <EOL> if has_replicates : <EOL> E . info ( "<STR_LIT>" ) <EOL> else : <EOL> E . info ( "<STR_LIT>" ) <EOL> dispersion_method = "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % <EOL> ( dispersion_method , fit_type , sharing_mode ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> cds = R ( '''<STR_LIT>''' ) <EOL> if dispersion_method == "<STR_LIT>" : <EOL> R . png ( '<STR_LIT>' % <EOL> outfile_prefix ) <EOL> R . plotDispEsts ( cds ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> elif not has_replicates : <EOL> pass <EOL> else : <EOL> dispersions = R ( '''<STR_LIT>''' ) <EOL> for dispersion in dispersions : <EOL> R . png ( '<STR_LIT>' % <EOL> ( outfile_prefix , dispersion ) ) <EOL> R . plotDispEsts ( cds , name = dispersion ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> deseqPlotSizeFactors ( '<STR_LIT>' % locals ( ) ) <EOL> deseqOutputSizeFactors ( "<STR_LIT>" % locals ( ) ) <EOL> deseqPlotPairs ( '<STR_LIT>' % locals ( ) ) <EOL> if dispersion_method not in ( "<STR_LIT>" , ) : <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> else : <EOL> R ( '''<STR_LIT>''' ) <EOL> vsd = R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> deseqPlotCorrelationHeatmap ( <EOL> '<STR_LIT>' % outfile_prefix , <EOL> vsd ) <EOL> deseqPlotPCA ( '<STR_LIT>' % outfile_prefix , <EOL> vsd ) <EOL> select = R . order ( R . rowMeans ( R . counts ( cds ) ) , decreasing = True ) <EOL> deseqPlotGeneHeatmap ( <EOL> '<STR_LIT>' % outfile_prefix , <EOL> R [ '<STR_LIT>' ] ( R . exprs ( vsd ) . rx ( select ) ) ) <EOL> deseqPlotGeneHeatmap ( <EOL> '<STR_LIT>' % outfile_prefix , <EOL> R [ '<STR_LIT>' ] ( R . exprs ( vsd ) . rx ( select [ : <NUM_LIT:200> ] ) ) ) <EOL> all_results = [ ] <EOL> for combination in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> control , treatment = combination <EOL> gfix = "<STR_LIT>" % ( control , treatment ) <EOL> outfile_groups_prefix = outfile_prefix + gfix <EOL> E . info ( ( "<STR_LIT>" <EOL> "<STR_LIT>" ) % <EOL> ( control , treatment ) ) <EOL> res = R ( '''<STR_LIT>''' % ( control , treatment ) ) <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> deseqPlotPvaluesAgainstRowsums ( <EOL> '<STR_LIT>' % locals ( ) ) <EOL> E . info ( "<STR_LIT>" % ( control , treatment ) ) <EOL> R ( '''<STR_LIT>''' % <EOL> ( treatment , control ) ) <EOL> R . png ( '''<STR_LIT>''' % <EOL> locals ( ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> select = R ( '''<STR_LIT>''' % fdr ) <EOL> if R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] > <NUM_LIT:0> : <EOL> E . info ( '<STR_LIT>' % <EOL> ( treatment , control , len ( select ) ) ) <EOL> data = R . exprs ( vsd ) . rx ( select ) <EOL> if not isinstance ( data , rpy2 . robjects . vectors . FloatVector ) : <EOL> order = R . order ( R . rowMeans ( data ) , decreasing = True ) <EOL> deseqPlotGeneHeatmap ( <EOL> '<STR_LIT>' % outfile_groups_prefix , <EOL> R [ '<STR_LIT>' ] ( data [ order ] ) , <EOL> Colv = False , <EOL> Rowv = True ) <EOL> else : <EOL> E . warn ( '<STR_LIT>' ) <EOL> else : <EOL> E . warn ( '<STR_LIT>' % fdr ) <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> if has_replicates : <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if R ( '''<STR_LIT>''' ) [ <NUM_LIT:0> ] > <NUM_LIT:0> : <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> else : <EOL> E . warn ( '<STR_LIT>' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> results , counts = deseqParseResults ( control , <EOL> treatment , <EOL> fdr = fdr ) <EOL> all_results += results <EOL> E . info ( counts ) <EOL> outf = IOTools . openFile ( <EOL> "<STR_LIT>" % locals ( ) , "<STR_LIT:w>" ) <EOL> outf . write ( "<STR_LIT>" % counts . asTable ( ) ) <EOL> outf . close ( ) <EOL> writeExpressionResults ( outfile , all_results ) <EOL> def runDESeq2 ( outfile , <EOL> outfile_prefix = "<STR_LIT>" , <EOL> fdr = <NUM_LIT:0.1> , <EOL> ref_group = None , <EOL> model = None , <EOL> contrasts = None , <EOL> plot = <NUM_LIT:1> , <EOL> ) : <EOL> """<STR_LIT>""" <EOL> R ( '''<STR_LIT>''' ) <EOL> if isinstance ( R ( '''<STR_LIT>''' ) , rpy2 . robjects . vectors . DataFrame ) : <EOL> E . info ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> mdata = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> else : <EOL> R ( '''<STR_LIT>''' ) <EOL> mdata = "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % mdata ) <EOL> if model : <EOL> assert contrasts , "<STR_LIT>" <EOL> terms = set ( [ x for x in re . split ( "<STR_LIT>" , model ) if x != '<STR_LIT>' ] ) <EOL> assert terms . issubset ( mdata ) , "<STR_LIT>" <EOL> else : <EOL> if mdata != "<STR_LIT>" : <EOL> E . warn ( "<STR_LIT>" ) <EOL> terms = [ "<STR_LIT>" , ] <EOL> model = "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % model ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> if plot == <NUM_LIT:1> : <EOL> for factor in terms : <EOL> outf = outfile_prefix + factor + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outf ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R . data ( '<STR_LIT>' ) <EOL> rlog_out = R ( '<STR_LIT>' ) <EOL> rlogoutf = outfile_prefix + "<STR_LIT>" <EOL> rlog_out . to_csv ( rlogoutf , sep = "<STR_LIT:\t>" , index = False ) <EOL> os . system ( "<STR_LIT>" % rlogoutf ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" ) <EOL> if contrasts : <EOL> contrasts = ( x . split ( "<STR_LIT::>" ) for x in contrasts . split ( "<STR_LIT:U+002C>" ) ) <EOL> else : <EOL> groups = R ( '''<STR_LIT>''' ) <EOL> contrasts = ( ( "<STR_LIT>" , ) + x for x in itertools . combinations ( groups , <NUM_LIT:2> ) ) <EOL> df_final = pandas . DataFrame ( ) <EOL> raw_final = pandas . DataFrame ( ) <EOL> all_results = [ ] <EOL> for combination in contrasts : <EOL> variable , control , treatment = combination <EOL> gfix = "<STR_LIT>" % ( variable , control , treatment ) <EOL> outfile_groups_prefix = outfile_prefix + gfix + "<STR_LIT>" <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> E . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( treatment , control , variable ) ) <EOL> if plot == <NUM_LIT:1> : <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> R ( '''<STR_LIT>''' % fdr ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( treatment , control , variable ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R . data ( '<STR_LIT>' ) <EOL> raw_out = R ( '<STR_LIT>' ) <EOL> df_out = deseq2ParseResults ( fdr , vsd = False ) <EOL> raw_out [ "<STR_LIT>" ] = [ treatment , ] * len ( df_out . index ) <EOL> raw_out [ "<STR_LIT>" ] = [ control , ] * len ( df_out . index ) <EOL> raw_out [ "<STR_LIT>" ] = [ variable , ] * len ( df_out . index ) <EOL> raw_final = raw_final . append ( raw_out , ignore_index = True ) <EOL> df_out . to_csv ( IOTools . openFile ( outfile_prefix + gfix + "<STR_LIT>" , "<STR_LIT:w>" ) , <EOL> sep = "<STR_LIT:\t>" , <EOL> index_label = "<STR_LIT>" ) <EOL> E . info ( "<STR_LIT>" <EOL> "<STR_LIT>" % ( treatment , control , variable ) ) <EOL> df_out . reset_index ( inplace = True ) <EOL> df_out . rename ( columns = { "<STR_LIT:index>" : "<STR_LIT>" } , inplace = True ) <EOL> df_final = df_final . append ( df_out , ignore_index = True ) <EOL> results = df_final . values . tolist ( ) <EOL> writeExpressionResults ( outfile , results ) <EOL> rawoutf = outfile_prefix + "<STR_LIT>" <EOL> raw_final . to_csv ( rawoutf , sep = "<STR_LIT:\t>" , index = False ) <EOL> os . system ( "<STR_LIT>" % rawoutf ) <EOL> Design = collections . namedtuple ( "<STR_LIT>" , ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) ) <EOL> def readDesignFile ( design_file ) : <EOL> '''<STR_LIT>''' <EOL> design = collections . OrderedDict ( ) <EOL> with IOTools . openFile ( design_file ) as inf : <EOL> for line in inf : <EOL> if line . startswith ( "<STR_LIT>" ) : <EOL> continue <EOL> track , include , group , pair = line . split ( "<STR_LIT:\t>" ) [ : <NUM_LIT:4> ] <EOL> if track in design : <EOL> raise ValueError ( "<STR_LIT>" % track ) <EOL> design [ track ] = Design . _make ( ( int ( include ) , group , pair ) ) <EOL> return design <EOL> def plotTagStats ( infile , design_file , outfile_prefix ) : <EOL> '''<STR_LIT>''' <EOL> loadTagData ( infile , design_file ) <EOL> nobservations , nsamples = filterTagData ( ) <EOL> if nobservations == <NUM_LIT:0> : <EOL> E . warn ( "<STR_LIT>" ) <EOL> return <EOL> if nsamples == <NUM_LIT:0> : <EOL> E . warn ( "<STR_LIT>" ) <EOL> return <EOL> groups , pairs , has_replicates , has_pairs = groupTagData ( ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R . png ( outfile_prefix + "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> R . png ( outfile_prefix + "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R ( '''<STR_LIT>''' ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def plotDETagStats ( infile , outfile_prefix , <EOL> additional_file = None , <EOL> join_columns = None , <EOL> additional_columns = None ) : <EOL> '''<STR_LIT>''' <EOL> table = pandas . read_csv ( IOTools . openFile ( infile ) , <EOL> sep = "<STR_LIT:\t>" ) <EOL> if additional_file is not None : <EOL> additional_table = pandas . read_csv ( <EOL> IOTools . openFile ( additional_file ) , <EOL> sep = "<STR_LIT:\t>" ) <EOL> table = pandas . merge ( table , <EOL> additional_table , <EOL> on = join_columns , <EOL> how = "<STR_LIT:left>" , <EOL> sort = False ) <EOL> table . reset_index ( inplace = True ) <EOL> table [ '<STR_LIT>' ] = numpy . log10 ( table [ '<STR_LIT>' ] + <NUM_LIT:1> ) <EOL> table [ '<STR_LIT>' ] = numpy . log10 ( table [ '<STR_LIT>' ] + <NUM_LIT:1> ) <EOL> table [ '<STR_LIT>' ] = numpy . array ( [ "<STR_LIT>" ] * len ( table ) ) <EOL> table . loc [ <EOL> ( table [ "<STR_LIT>" ] > <NUM_LIT:0> ) & ( table [ "<STR_LIT>" ] == <NUM_LIT:1> ) , "<STR_LIT>" ] = "<STR_LIT>" <EOL> table . loc [ <EOL> ( table [ "<STR_LIT>" ] < <NUM_LIT:0> ) & ( table [ "<STR_LIT>" ] == <NUM_LIT:1> ) , "<STR_LIT>" ] = "<STR_LIT>" <EOL> def _dplot ( table , outfile , column ) : <EOL> plot = ggplot . ggplot ( <EOL> ggplot . aes ( column , <EOL> colour = '<STR_LIT>' , <EOL> fill = '<STR_LIT>' ) , <EOL> data = table ) + ggplot . geom_density ( alpha = <NUM_LIT:0.5> ) <EOL> try : <EOL> ggplot . ggsave ( filename = outfile , plot = plot ) <EOL> except Exception , msg : <EOL> E . warn ( "<STR_LIT>" % ( column , msg ) ) <EOL> def _bplot ( table , outfile , column ) : <EOL> plot = ggplot . ggplot ( <EOL> ggplot . aes ( x = column , y = '<STR_LIT>' ) , <EOL> data = table ) + ggplot . geom_boxplot ( ) <EOL> try : <EOL> ggplot . ggsave ( filename = outfile , plot = plot ) <EOL> except ValueError , msg : <EOL> E . warn ( msg ) <EOL> _dplot ( table , <EOL> outfile_prefix + "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> _dplot ( table , <EOL> outfile_prefix + "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> _bplot ( table , <EOL> outfile_prefix + "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> _bplot ( table , <EOL> outfile_prefix + "<STR_LIT>" , <EOL> "<STR_LIT>" ) <EOL> if additional_columns : <EOL> for column in additional_columns : <EOL> _dplot ( table , <EOL> outfile_prefix + "<STR_LIT>" % column , <EOL> column ) <EOL> _bplot ( table , <EOL> outfile_prefix + "<STR_LIT>" % column , <EOL> column ) <EOL> return <EOL> def runMockAnalysis ( outfile , <EOL> outfile_prefix , <EOL> ref_group = None , <EOL> ref_regex = None , <EOL> pseudo_counts = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> groups , pairs , has_replicates , has_pairs = groupTagData ( ref_group , <EOL> ref_regex ) <EOL> all_results = [ ] <EOL> for combination in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> control , treatment = combination <EOL> R ( '''<STR_LIT>''' % <EOL> control ) <EOL> R ( '''<STR_LIT>''' % <EOL> treatment ) <EOL> if pseudo_counts : <EOL> R ( '''<STR_LIT>''' % pseudo_counts ) <EOL> R ( '''<STR_LIT>''' % pseudo_counts ) <EOL> R ( '''<STR_LIT>''' ) <EOL> results = [ ] <EOL> for identifier , treatment_count , control_count , foldchange in zip ( R ( '''<STR_LIT>''' ) , <EOL> R ( '''<STR_LIT>''' ) , <EOL> R ( '''<STR_LIT>''' ) , <EOL> R ( '''<STR_LIT>''' ) ) : <EOL> try : <EOL> log2fold = math . log ( foldchange ) <EOL> except ValueError : <EOL> log2fold = "<STR_LIT>" <EOL> results . append ( GeneExpressionResult . _make ( ( <EOL> identifier , <EOL> treatment , <EOL> treatment_count , <EOL> <NUM_LIT:0> , <EOL> control , <EOL> control_count , <EOL> <NUM_LIT:0> , <EOL> <NUM_LIT:1> , <EOL> <NUM_LIT:1> , <EOL> log2fold , <EOL> foldchange , <EOL> log2fold , <EOL> "<STR_LIT:0>" , <EOL> "<STR_LIT:OK>" ) ) ) <EOL> all_results . extend ( results ) <EOL> writeExpressionResults ( outfile , all_results ) <EOL> def outputTagSummary ( filename_tags , <EOL> outfile , <EOL> output_filename_pattern , <EOL> filename_design = None ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" % filename_tags ) <EOL> if filename_design is not None : <EOL> loadTagData ( filename_tags , filename_design ) <EOL> nobservations , nsamples = filterTagData ( ) <EOL> else : <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> nobservations , nsamples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( nobservations , nsamples ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> filter_min_counts_per_sample = <NUM_LIT:1> <EOL> empty_samples = tuple ( <EOL> R ( '''<STR_LIT>''' % filter_min_counts_per_sample ) ) <EOL> sample_names = R ( '''<STR_LIT>''' ) <EOL> nempty_samples = sum ( empty_samples ) <EOL> if nempty_samples : <EOL> E . warn ( "<STR_LIT>" % <EOL> ( nempty_samples , <EOL> "<STR_LIT:U+002C>" . join ( [ sample_names [ x ] <EOL> for x , y in enumerate ( empty_samples ) if y ] ) ) ) <EOL> R ( '''<STR_LIT>''' % <EOL> filter_min_counts_per_sample ) <EOL> nobservations , nsamples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> E . debug ( "<STR_LIT>" % R ( '''<STR_LIT>''' ) ) <EOL> nrows , ncolumns = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> outfile . write ( "<STR_LIT>" ) <EOL> outfile . write ( "<STR_LIT>" % nobservations ) <EOL> outfile . write ( "<STR_LIT>" % nsamples ) <EOL> R ( '''<STR_LIT>''' ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R ( '''<STR_LIT>''' % <EOL> locals ( ) ) <EOL> E . info ( "<STR_LIT>" ) <EOL> R ( '''<STR_LIT>''' ) <EOL> for x in range ( <NUM_LIT:0> , <NUM_LIT:20> ) : <EOL> nempty = tuple ( R ( '''<STR_LIT>''' % x ) ) [ <NUM_LIT:0> ] <EOL> outfile . write ( "<STR_LIT>" % <EOL> ( x , nempty , <NUM_LIT> * nempty / nrows ) ) <EOL> E . info ( "<STR_LIT>" % tuple ( R ( '''<STR_LIT>''' ) ) ) <EOL> observations , samples = tuple ( R ( '''<STR_LIT>''' ) ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( observations , samples ) ) <EOL> R ( '''<STR_LIT>''' ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R . png ( outfilename , width = <NUM_LIT> , height = <NUM_LIT> ) <EOL> plotPairs ( ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R . svg ( outfilename ) <EOL> plotCorrelationHeatmap ( method = "<STR_LIT>" ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R . svg ( outfilename ) <EOL> plotPCA ( groups = False ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> R ( '''<STR_LIT>''' ) <EOL> outfilename = output_filename_pattern + "<STR_LIT>" <EOL> E . info ( "<STR_LIT>" % outfilename ) <EOL> R . svg ( outfilename ) <EOL> try : <EOL> R ( '''<STR_LIT>''' ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> def dumpTagData ( filename_tags , filename_design , outfile ) : <EOL> '''<STR_LIT>''' <EOL> if outfile == sys . stdout : <EOL> outfilename = "<STR_LIT>" <EOL> else : <EOL> outfilename = outfile . name <EOL> loadTagData ( filename_tags , filename_design ) <EOL> nobservations , nsamples = filterTagData ( ) <EOL> R ( '''<STR_LIT>''' % locals ( ) ) <EOL> def runTTest ( outfile , <EOL> outfile_prefix , <EOL> fdr = <NUM_LIT:0.1> , <EOL> ref_group = None , <EOL> ref_regex = None ) : <EOL> '''<STR_LIT>''' <EOL> groups , pairs , has_replicates , has_pairs = groupTagData ( ref_group , <EOL> ref_regex ) <EOL> results = [ ] <EOL> for combination in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> control , treatment = combination <EOL> r = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> for test_id , ttest in zip ( r . names , r ) : <EOL> treatment_mean , control_mean = tuple ( ttest . rx2 ( '<STR_LIT>' ) ) <EOL> fold_change = treatment_mean / control_mean <EOL> pvalue = tuple ( ttest . rx2 ( '<STR_LIT>' ) ) [ <NUM_LIT:0> ] <EOL> significant = ( <NUM_LIT:0> , <NUM_LIT:1> ) [ pvalue < fdr ] <EOL> results . append ( GeneExpressionResult . _make ( ( test_id , <EOL> treatment , <EOL> treatment_mean , <EOL> <NUM_LIT:0> , <EOL> control , <EOL> control_mean , <EOL> <NUM_LIT:0> , <EOL> pvalue , <EOL> pvalue , <EOL> numpy . log2 ( fold_change ) , <EOL> fold_change , <EOL> numpy . log2 ( fold_change ) , <EOL> significant , <EOL> "<STR_LIT:OK>" ) ) ) <EOL> writeExpressionResults ( outfile , results ) <EOL> def loadTagDataPandas ( tags_filename , design_filename ) : <EOL> '''<STR_LIT>''' <EOL> E . info ( "<STR_LIT>" % tags_filename ) <EOL> inf = IOTools . openFile ( tags_filename ) <EOL> counts_table = pandas . read_csv ( inf , <EOL> sep = "<STR_LIT:\t>" , <EOL> index_col = <NUM_LIT:0> , <EOL> comment = "<STR_LIT:#>" ) <EOL> inf . close ( ) <EOL> E . info ( "<STR_LIT>" % <EOL> counts_table . shape ) <EOL> E . debug ( "<STR_LIT>" % list ( counts_table . columns ) ) <EOL> inf = IOTools . openFile ( design_filename ) <EOL> design_table = pandas . read_csv ( inf , sep = "<STR_LIT:\t>" , index_col = <NUM_LIT:0> ) <EOL> inf . close ( ) <EOL> E . debug ( "<STR_LIT>" % list ( design_table . index ) ) <EOL> missing = set ( counts_table . columns ) . difference ( design_table . index ) <EOL> if missing : <EOL> E . warn ( "<STR_LIT>" % missing ) <EOL> design_table = design_table [ design_table [ "<STR_LIT>" ] != <NUM_LIT:0> ] <EOL> E . debug ( "<STR_LIT>" % list ( design_table . index ) ) <EOL> counts_table = counts_table [ list ( design_table . index ) ] <EOL> E . info ( "<STR_LIT>" % <EOL> counts_table . shape ) <EOL> return counts_table , design_table <EOL> def filterTagDataPandas ( counts_table , <EOL> design_table , <EOL> filter_min_counts_per_row = <NUM_LIT:1> , <EOL> filter_min_counts_per_sample = <NUM_LIT:10> , <EOL> filter_percentile_rowsums = <NUM_LIT:0> ) : <EOL> '''<STR_LIT>''' <EOL> max_counts_per_row = counts_table . max ( <NUM_LIT:1> ) <EOL> counts_table = counts_table [ <EOL> max_counts_per_row >= filter_min_counts_per_row ] <EOL> observations , samples = counts_table . shape <EOL> E . info ( "<STR_LIT>" % <EOL> ( observations , samples ) ) <EOL> max_counts_per_sample = counts_table . max ( ) <EOL> empty_samples = max_counts_per_sample < filter_min_counts_per_sample <EOL> sample_names = counts_table . columns <EOL> nempty_samples = sum ( empty_samples ) <EOL> if nempty_samples : <EOL> E . warn ( "<STR_LIT>" % <EOL> ( nempty_samples , <EOL> "<STR_LIT:U+002C>" . join ( [ sample_names [ x ] for x , y in <EOL> enumerate ( empty_samples ) if y ] ) ) ) <EOL> raise NotImplementedError ( "<STR_LIT>" ) <EOL> if filter_percentile_rowsums > <NUM_LIT:0> : <EOL> percentile = float ( filter_percentile_rowsums ) / <NUM_LIT> <EOL> sum_counts = counts_table . sum ( <NUM_LIT:1> ) <EOL> take = sum_counts > sum_counts . quantile ( percentile ) <EOL> E . info ( "<STR_LIT>" % <EOL> ( filter_percentile_rowsums , <EOL> sum ( take ) , <EOL> len ( take ) - sum ( take ) ) ) <EOL> counts_table = counts_table [ take ] <EOL> return counts_table <EOL> def identifyVariablesPandas ( design_table ) : <EOL> design_table = design_table [ design_table [ "<STR_LIT>" ] != <NUM_LIT:0> ] <EOL> conds = design_table [ '<STR_LIT>' ] . tolist ( ) <EOL> pairs = design_table [ '<STR_LIT>' ] . tolist ( ) <EOL> if len ( design_table . columns ) > <NUM_LIT:3> : <EOL> factors = design_table . iloc [ : , <NUM_LIT:3> : ] <EOL> else : <EOL> factors = None <EOL> return conds , pairs , factors <EOL> def checkTagGroupsPandas ( design_table , ref_group = None ) : <EOL> '''<STR_LIT>''' <EOL> conds , pairs , factors = identifyVariablesPandas ( design_table ) <EOL> groups = list ( set ( conds ) ) <EOL> if factors : <EOL> E . warn ( "<STR_LIT>" <EOL> "<STR_LIT>" , factors ) <EOL> else : <EOL> pass <EOL> max_per_group = max ( [ conds . count ( x ) for x in groups ] ) <EOL> has_replicates = max_per_group >= <NUM_LIT:2> <EOL> npairs = len ( set ( pairs ) ) <EOL> has_pairs = npairs == <NUM_LIT:2> <EOL> if has_pairs : <EOL> min_per_pair = min ( [ pairs . count ( x ) for x in set ( pairs ) ] ) <EOL> has_pairs = min_per_pair >= <NUM_LIT:2> <EOL> return groups , pairs , conds , factors , has_replicates , has_pairs <EOL> ResultColumns = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT:status>" ] <EOL> ResultColumns_dtype = { "<STR_LIT>" : object , "<STR_LIT>" : object , <EOL> "<STR_LIT>" : float , "<STR_LIT>" : <EOL> float , "<STR_LIT>" : object , "<STR_LIT>" : <EOL> float , "<STR_LIT>" : float , "<STR_LIT>" : float , <EOL> "<STR_LIT>" : float , "<STR_LIT>" : float , "<STR_LIT>" : <EOL> float , "<STR_LIT>" : float , <EOL> "<STR_LIT>" : int , "<STR_LIT:status>" : object } <EOL> def makeEmptyDataFrameDict ( ) : <EOL> return { key : [ ] for key in ResultColumns } <EOL> def runTTestPandas ( counts_table , <EOL> design_table , <EOL> outfile , <EOL> outfile_prefix , <EOL> fdr , <EOL> ref_group = None ) : <EOL> '''<STR_LIT>''' <EOL> stats = importr ( '<STR_LIT>' ) <EOL> ( groups , pairs , conds , factors , has_replicates , <EOL> has_pairs ) = checkTagGroupsPandas ( design_table , ref_group ) <EOL> df_dict = makeEmptyDataFrameDict ( ) <EOL> for combination in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> control , treatment = combination <EOL> n_rows = counts_table . shape [ <NUM_LIT:0> ] <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( control , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( treatment , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( counts_table . index . tolist ( ) ) <EOL> c_keep = [ x == control for x in conds ] <EOL> control_counts = counts_table . iloc [ : , c_keep ] <EOL> t_keep = [ x == treatment for x in conds ] <EOL> treatment_counts = counts_table . iloc [ : , t_keep ] <EOL> c_mean = control_counts . mean ( axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( c_mean ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( control_counts . std ( axis = <NUM_LIT:1> ) ) <EOL> t_mean = treatment_counts . mean ( axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( t_mean ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( treatment_counts . std ( axis = <NUM_LIT:1> ) ) <EOL> t , prob = ttest_ind ( control_counts , treatment_counts , axis = <NUM_LIT:1> ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( prob ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( t_mean / c_mean ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( <EOL> list ( stats . p_adjust ( FloatVector ( df_dict [ "<STR_LIT>" ] ) , method = '<STR_LIT>' ) ) ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( <EOL> [ int ( x < fdr ) for x in df_dict [ "<STR_LIT>" ] ] ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( list ( numpy . log2 ( df_dict [ "<STR_LIT>" ] ) ) ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( list ( numpy . log2 ( df_dict [ "<STR_LIT>" ] ) ) ) <EOL> df_dict [ "<STR_LIT:status>" ] . extend ( ( "<STR_LIT:OK>" , ) * n_rows ) <EOL> results = pandas . DataFrame ( df_dict ) <EOL> results . set_index ( "<STR_LIT>" , inplace = True ) <EOL> results . to_csv ( outfile , sep = "<STR_LIT:\t>" , header = True , index = True ) <EOL> def plotCorrelationHeatmapMatplot ( counts , outfile , method = "<STR_LIT>" , <EOL> cor_method = "<STR_LIT>" ) : <EOL> '''<STR_LIT>''' <EOL> heatmap_cdict_b_to_y = { <EOL> '<STR_LIT>' : ( ( <NUM_LIT:0.0> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT:1.> , <NUM_LIT> , <NUM_LIT> ) ) , <EOL> '<STR_LIT>' : ( ( <NUM_LIT:0.0> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT:1.> , <NUM_LIT> , <NUM_LIT> ) ) , <EOL> '<STR_LIT>' : ( ( <NUM_LIT:0.0> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , ( <NUM_LIT:1.> , <NUM_LIT> , <NUM_LIT> ) ) } <EOL> cm = matplotlib . colors . LinearSegmentedColormap ( <EOL> '<STR_LIT>' , heatmap_cdict_b_to_y , <NUM_LIT> ) <EOL> df = counts . corr ( method = cor_method ) <EOL> plt . pcolor ( np . array ( df ) , cmap = cm ) <EOL> plt . colorbar ( ) <EOL> plt . title ( "<STR_LIT>" % locals ( ) ) <EOL> plt . yticks ( np . arange ( <NUM_LIT:0.5> , len ( df . index ) , <NUM_LIT:1> ) , df . index ) <EOL> plt . xticks ( np . arange ( <NUM_LIT:0.5> , len ( df . columns ) , <NUM_LIT:1> ) , df . columns , rotation = <NUM_LIT> ) <EOL> plt . tight_layout ( ) <EOL> plt . savefig ( outfile ) <EOL> def runEdgeRPandas ( counts , <EOL> design_table , <EOL> outfile , <EOL> outfile_prefix = "<STR_LIT>" , <EOL> fdr = <NUM_LIT:0.1> , <EOL> prefix = "<STR_LIT>" , <EOL> dispersion = None , <EOL> ref_group = None ) : <EOL> '''<STR_LIT>''' <EOL> R ( '''<STR_LIT>''' ) <EOL> ( groups , pairs , conds , factors , has_replicates , <EOL> has_pairs ) = checkTagGroupsPandas ( design_table , ref_group ) <EOL> if not has_replicates and dispersion is None : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> plotCorrelationHeatmapMatplot ( counts , <EOL> '<STR_LIT>' % locals ( ) , <EOL> cor_method = "<STR_LIT>" ) <EOL> E . info ( '<STR_LIT>' % <EOL> ( groups , pairs , has_replicates , has_pairs ) ) <EOL> r_counts = pandas2ri . py2ri ( counts ) <EOL> passDFtoRGlobalEnvironment = R ( '''<STR_LIT>''' ) <EOL> passDFtoRGlobalEnvironment ( r_counts ) <EOL> if has_pairs : <EOL> first = True <EOL> pairs_df = pandas . DataFrame ( ) <EOL> nrows = len ( counts . index ) <EOL> n = <NUM_LIT:0> <EOL> for g1 , g2 in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> keep_a = [ x == g1 for x in conds ] <EOL> counts_a = counts . iloc [ : , keep_a ] <EOL> keep_b = [ x == g2 for x in conds ] <EOL> counts_b = counts . iloc [ : , keep_b ] <EOL> index = range ( n , n + nrows ) <EOL> n += nrows <EOL> a = counts_a . sum ( axis = <NUM_LIT:1> ) <EOL> b = counts_b . sum ( axis = <NUM_LIT:1> ) <EOL> diff = a - b <EOL> diff . sort ( ) <EOL> temp_df = pandas . DataFrame ( { "<STR_LIT>" : np . cumsum ( diff ) . tolist ( ) , <EOL> "<STR_LIT>" : "<STR_LIT>" . join ( [ g1 , g2 ] ) , <EOL> "<STR_LIT:id>" : range ( <NUM_LIT:0> , nrows ) } , <EOL> index = index ) <EOL> pairs_df = pairs_df . append ( temp_df ) <EOL> plot_pairs = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> r_pairs_df = pandas2ri . py2ri ( pairs_df ) <EOL> plot_pairs ( r_pairs_df ) <EOL> first = True <EOL> legend = [ ] <EOL> n = <NUM_LIT:0> <EOL> pairs_in_groups_df = pandas . DataFrame ( ) <EOL> for pair in set ( pairs ) : <EOL> for g1 , g2 in itertools . combinations ( groups , <NUM_LIT:2> ) : <EOL> key = "<STR_LIT>" % ( pair , g1 , g2 ) <EOL> legend . append ( key ) <EOL> keep_a = [ x == pair and y == g1 for x , y in zip ( pairs , conds ) ] <EOL> keep_b = [ x == pair and y == g2 for x , y in zip ( pairs , conds ) ] <EOL> if sum ( keep_a ) > <NUM_LIT:0> and sum ( keep_b ) > <NUM_LIT:0> : <EOL> counts_a = counts . iloc [ : , keep_a ] <EOL> counts_b = counts . iloc [ : , keep_b ] <EOL> index = range ( n , n + nrows ) <EOL> n += nrows <EOL> a = counts_a . sum ( axis = <NUM_LIT:1> ) <EOL> b = counts_b . sum ( axis = <NUM_LIT:1> ) <EOL> diff = a - b <EOL> diff . sort ( ) <EOL> comparison = "<STR_LIT>" % ( pair , g1 , g2 ) <EOL> temp_df = pandas . DataFrame ( { "<STR_LIT>" : np . cumsum ( diff ) . tolist ( ) , <EOL> "<STR_LIT>" : comparison , <EOL> "<STR_LIT:id>" : range ( <NUM_LIT:0> , nrows ) } , <EOL> index = index ) <EOL> pairs_in_groups_df = pairs_in_groups_df . append ( temp_df ) <EOL> plot_pairs_in_groups = R ( '''<STR_LIT>''' % locals ( ) ) <EOL> r_pairs_in_groups_df = pandas2ri . py2ri ( pairs_in_groups_df ) <EOL> plot_pairs_in_groups ( r_pairs_in_groups_df ) <EOL> r_counts = pandas2ri . py2ri ( counts ) <EOL> r_groups = ro . StrVector ( conds ) <EOL> r_pairs = ro . StrVector ( pairs ) <EOL> r_has_pairs = ro . default_py2ri ( has_pairs ) <EOL> r_has_replicates = ro . default_py2ri ( has_replicates ) <EOL> if dispersion is not None : <EOL> r_dispersion = ro . default_py2ri ( dispersion ) <EOL> else : <EOL> r_dispersion = ro . default_py2ri ( False ) <EOL> if ref_group is not None : <EOL> r_ref_group = ro . default_py2ri ( ref_group ) <EOL> else : <EOL> r_ref_group = ro . default_py2ri ( groups [ <NUM_LIT:0> ] ) <EOL> buildDGEList = R ( '''<STR_LIT>''' ) <EOL> r_countsTable = buildDGEList ( r_counts , r_groups , r_ref_group ) <EOL> R . png ( '''<STR_LIT>''' % locals ( ) ) <EOL> try : <EOL> MDSplot = R ( '''<STR_LIT>''' ) <EOL> MDSplot ( r_counts ) <EOL> except rpy2 . rinterface . RRuntimeError , msg : <EOL> E . warn ( "<STR_LIT>" % msg ) <EOL> R [ '<STR_LIT>' ] ( ) <EOL> buildDesign = R ( '''<STR_LIT>''' ) <EOL> r_design = buildDesign ( r_countsTable , r_has_pairs , r_pairs ) <EOL> fitModel = R ( '''<STR_LIT>''' ) <EOL> r_fit = fitModel ( r_countsTable , r_design , r_has_replicates , r_dispersion ) <EOL> E . info ( "<STR_LIT>" ) <EOL> lrtTest = R ( '''<STR_LIT>''' ) <EOL> r_lrt = lrtTest ( r_fit , outfile_prefix ) <EOL> extractTable = R ( '''<STR_LIT>''' ) <EOL> r_lrt_table = extractTable ( r_lrt ) <EOL> outputCPMTable = R ( '''<STR_LIT>''' ) <EOL> outputCPMTable ( r_countsTable , outfile_prefix , r_lrt ) <EOL> if len ( groups ) == <NUM_LIT:2> : <EOL> plotSmear = R ( '''<STR_LIT>''' % "<STR_LIT>" . join ( groups ) ) <EOL> lrt_table = pandas2ri . ri2py ( r_lrt_table ) <EOL> n_rows = lrt_table . shape [ <NUM_LIT:0> ] <EOL> df_dict = makeEmptyDataFrameDict ( ) <EOL> stats = importr ( '<STR_LIT>' ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( groups [ <NUM_LIT:0> ] , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( groups [ <NUM_LIT:1> ] , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( lrt_table . index ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( lrt_table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( lrt_table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( <NUM_LIT:0> , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( ( <NUM_LIT:0> , ) * n_rows ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( lrt_table [ '<STR_LIT>' ] ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( <EOL> list ( stats . p_adjust ( FloatVector ( df_dict [ "<STR_LIT>" ] ) , method = '<STR_LIT>' ) ) ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( <EOL> [ int ( float ( x ) < fdr ) for x in df_dict [ "<STR_LIT>" ] ] ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( list ( numpy . log2 ( lrt_table [ '<STR_LIT>' ] ) ) ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( list ( numpy . log2 ( lrt_table [ '<STR_LIT>' ] ) ) ) <EOL> df_dict [ "<STR_LIT>" ] . extend ( [ math . pow ( <NUM_LIT:2> , float ( x ) ) for x in lrt_table [ '<STR_LIT>' ] ] ) <EOL> df_dict [ "<STR_LIT:status>" ] . extend ( ( "<STR_LIT:OK>" , ) * n_rows ) <EOL> results = pandas . DataFrame ( df_dict ) <EOL> results . set_index ( "<STR_LIT>" , inplace = True ) <EOL> results . to_csv ( outfile , sep = "<STR_LIT:\t>" , header = True , index = True ) <EOL> counts = E . Counter ( ) <EOL> counts . signficant = sum ( results [ '<STR_LIT>' ] ) <EOL> counts . insignficant = ( len ( results [ '<STR_LIT>' ] ) - counts . signficant ) <EOL> counts . all_over = sum ( [ x > <NUM_LIT:0> for x in results [ '<STR_LIT>' ] ] ) <EOL> counts . all_under = sum ( [ x < <NUM_LIT:0> for x in results [ '<STR_LIT>' ] ] ) <EOL> counts . signficant_over = sum ( [ results [ '<STR_LIT>' ] [ x ] == <NUM_LIT:1> and <EOL> results [ '<STR_LIT>' ] [ x ] > <NUM_LIT:1> for <EOL> x in range ( <NUM_LIT:0> , n_rows ) ] ) <EOL> counts . signficant_under = sum ( [ results [ '<STR_LIT>' ] [ x ] == <NUM_LIT:1> and <EOL> results [ '<STR_LIT>' ] [ x ] < <NUM_LIT:1> for <EOL> x in range ( <NUM_LIT:0> , n_rows ) ] ) <EOL> outf = IOTools . openFile ( "<STR_LIT>" % locals ( ) , "<STR_LIT:w>" ) <EOL> outf . write ( "<STR_LIT>" % counts . asTable ( ) ) <EOL> outf . close ( ) <EOL> def outputSpikeIns ( filename_tags , <EOL> outfile , <EOL> output_filename_pattern , <EOL> filename_design = None , <EOL> foldchange_max = <NUM_LIT> , <EOL> expression_max = <NUM_LIT> , <EOL> max_counts_per_bin = <NUM_LIT:100> , <EOL> expression_bin_width = <NUM_LIT:0.5> , <EOL> foldchange_bin_width = <NUM_LIT:0.5> , <EOL> iterations = <NUM_LIT:1> ) : <EOL> E . info ( "<STR_LIT>" % filename_tags ) <EOL> if filename_design is not None : <EOL> counts_table , design_table = loadTagDataPandas ( <EOL> filename_tags , filename_design ) <EOL> counts_table = filterTagDataPandas ( counts_table , design_table ) <EOL> else : <EOL> raise NotImplementedError ( "<STR_LIT>" ) <EOL> nobservations , nsamples = counts_table . shape <EOL> groups = list ( set ( design_table [ "<STR_LIT>" ] ) ) <EOL> if len ( groups ) != <NUM_LIT:2> : <EOL> raise NotImplementedError ( <EOL> "<STR_LIT>" ) <EOL> group1 = counts_table [ design_table [ design_table . group == groups [ <NUM_LIT:0> ] ] . index ] <EOL> group2 = counts_table [ design_table [ design_table . group == groups [ <NUM_LIT:1> ] ] . index ] <EOL> outfile . write ( "<STR_LIT>" % <EOL> ( "<STR_LIT:\t>" . join ( group1 . columns ) , "<STR_LIT:\t>" . join ( group2 . columns ) ) ) <EOL> outf_info = IOTools . openFile ( output_filename_pattern + "<STR_LIT>" , "<STR_LIT:w>" ) <EOL> outf_info . write ( "<STR_LIT>" ) <EOL> group1 = group1 . as_matrix ( ) <EOL> group2 = group2 . as_matrix ( ) <EOL> expression_bins = numpy . arange ( <NUM_LIT:0> , expression_max , expression_bin_width ) <EOL> fold_change_bins = numpy . arange ( - foldchange_max , <EOL> foldchange_max , foldchange_bin_width ) <EOL> E . info ( "<STR_LIT>" % expression_bins ) <EOL> E . info ( "<STR_LIT>" % fold_change_bins ) <EOL> output_counts = numpy . zeros ( <EOL> ( len ( expression_bins ) + <NUM_LIT:1> , len ( fold_change_bins ) + <NUM_LIT:1> ) ) <EOL> for iteration in range ( iterations ) : <EOL> group1 = numpy . random . permutation ( group1 ) <EOL> group2 = numpy . random . permutation ( group2 ) <EOL> group1_mean = group1 . mean ( <NUM_LIT:1> ) <EOL> group2_mean = group2 . mean ( <NUM_LIT:1> ) <EOL> l10average = numpy . log ( ( group1_mean + group2_mean ) / <NUM_LIT> ) <EOL> l2fold = numpy . log2 ( ( group1_mean + <NUM_LIT:1.0> ) / ( group2_mean + <NUM_LIT:1.0> ) ) <EOL> l10average_idx = numpy . digitize ( l10average , expression_bins ) <EOL> l2fold_idx = numpy . digitize ( l2fold , fold_change_bins ) <EOL> interval_id = <NUM_LIT:0> <EOL> for idx , coord in enumerate ( zip ( l10average_idx , l2fold_idx ) ) : <EOL> if output_counts [ coord ] >= max_counts_per_bin : <EOL> continue <EOL> output_counts [ coord ] += <NUM_LIT:1> <EOL> outf_info . write ( "<STR_LIT>" % ( interval_id , <EOL> l10average [ idx ] , <EOL> l2fold [ idx ] ) ) <EOL> outfile . write ( "<STR_LIT>" % <EOL> ( interval_id , <EOL> "<STR_LIT:\t>" . join ( <EOL> map ( str , list ( group1 [ idx ] ) ) ) , <EOL> "<STR_LIT:\t>" . join ( map ( str , list ( group2 [ idx ] ) ) ) ) ) <EOL> interval_id += <NUM_LIT:1> <EOL> outf_info . close ( ) <EOL> df = pandas . DataFrame ( output_counts ) <EOL> df . index = list ( expression_bins ) + [ expression_max ] <EOL> df . columns = list ( fold_change_bins ) + [ foldchange_max ] <EOL> df . to_csv ( IOTools . openFile ( output_filename_pattern + "<STR_LIT>" , "<STR_LIT:w>" ) , <EOL> sep = "<STR_LIT:\t>" ) <EOL> E . info ( "<STR_LIT>" % interval_id ) </s>
