CNSAnalysis/CNS_commandline.py
def CNS_RetentionvsGeneticDistance(CNS_bed, reference_species, tree_file, species_file):
    """Reference should be CNS bed species"""
    import plotly.graph_objs as go
    import plotly.offline as py
    import matplotlib.pyplot as plt
    from collections import defaultdict
    from scipy import stats
    tree_mat = tree2matrix(tree_file)
    #genetic_distances = tree_mat[reference_species]
    species = tree_mat.index.values
    print tree_mat
    print species
    #FIXME binarize CNS matrix
    included_species = os.popen('cat '+species_file).read().splitlines()
    if not included_species:
        included_species = species
    refs = reference_species.split(',')
    c = ['hsl(' + str(h) + ',50%' + ',50%)' for h in np.linspace(0, 360, len(set(refs)) + 2)]
    plots = []
    retention_rates = {ref:defaultdict(list) for ref in refs}
    x_f, y_f = [],[]
    for i,ref in enumerate(refs):
        CNS_bed = CNS_bed[:CNS_bed.rfind('/')+1]+ref+'_AllCNSElements.bed'
        CNS_matrix = CNS_matrix_generation(CNS_bed).drop_duplicates()
        print CNS_matrix
        CNS_matrix = CNS_matrix.iloc[:,2:]
        print CNS_matrix
        total_CNS = CNS_matrix.shape[0]
        plt.figure()
        x = []
        y = []
        name = []
        for sp in species:
            if sp != reference_species:
                distance = tree_mat[sp][ref] #FIXME debug
                retention_rate = float(len(CNS_matrix[sp][CNS_matrix[sp] > 0].as_matrix()))/total_CNS
                plt.scatter(retention_rate,distance,s=3,label=sp)
                x_f.append(distance)
                y_f.append(retention_rate)
                if sp in included_species:
                    retention_rates[ref][sp] = retention_rate
                    y.append(retention_rate)
                    x.append(distance)
                    name.append(sp)
        plots.append(go.Scatter(x=x, y=y, name=ref,
                     mode='markers',
                     marker=dict(size=8,color=c[i]), text=name))
        print(x,y,name)
        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
        plt.xlabel('CNS Retention Rate')
        plt.ylabel('Genetic Distance')
        plt.title(reference_species+' CNS Retention Rate vs Genetic Distance, %s'%ref)
        plt.savefig('CNS_retention_vs_distance_%s.pdf'%ref,dpi=300)
        #FIXME fix local PCA
    py.plot(go.Figure(data=plots,layout=go.Layout(title='Genetic Distance to Reference vs CNS Retention',yaxis=dict(title='Percent of CNS Elements Retained With Respect to Reference'),xaxis=dict(title='Genetic Distance from Reference'))),filename='CNS_retention_vs_distance.html',auto_open=False)
    plots = []
    for ref in refs:
        plots.append(go.Bar(x=retention_rates[ref].keys(),y=retention_rates[ref].values(),name=ref))# FIXME finish, add tests for independence between genetic distance and retention, check if polyploid reduction in retention
    py.plot(go.Figure(data=plots,layout=go.Layout({'barmode':'group'})),filename='CNS_retention.html',auto_open=False)
    x_f = np.nan_to_num(np.array(x_f))
    y_f = np.nan_to_num(np.array(y_f))
    print x_f,y_f
    p = np.polyfit(x_f,y_f,1)
    f = lambda x: np.polyval(p,x)
    expected = f(x_f)
    chi_sq = sum((y_f - expected)**2/expected)
    print('X^2',chi_sq,'p',1-stats.chi2.cdf(chi_sq,len(x_f)-2),'r',stats.pearsonr(x_f,y_f))
    # FIXME Chi-square doesn't matter if you do not have any uncertainty about the y-data point
Reference should be CNS bed species
Reference should be CNS bed species
Reference should be CNS bed species
def CNS_RetentionvsGeneticDistance(CNS_bed, reference_species, tree_file, species_file):
    import plotly.graph_objs as go
    import plotly.offline as py
    import matplotlib.pyplot as plt
    from collections import defaultdict
    from scipy import stats
    tree_mat = tree2matrix(tree_file)
    species = tree_mat.index.values
    print tree_mat
    print species
    included_species = os.popen('cat '+species_file).read().splitlines()
    if not included_species:
        included_species = species
    refs = reference_species.split(',')
    c = ['hsl(' + str(h) + ',50%' + ',50%)' for h in np.linspace(0, 360, len(set(refs)) + 2)]
    plots = []
    retention_rates = {ref:defaultdict(list) for ref in refs}
    x_f, y_f = [],[]
    for i,ref in enumerate(refs):
        CNS_bed = CNS_bed[:CNS_bed.rfind('/')+1]+ref+'_AllCNSElements.bed'
        CNS_matrix = CNS_matrix_generation(CNS_bed).drop_duplicates()
        print CNS_matrix
        CNS_matrix = CNS_matrix.iloc[:,2:]
        print CNS_matrix
        total_CNS = CNS_matrix.shape[0]
        plt.figure()
        x = []
        y = []
        name = []
        for sp in species:
            if sp != reference_species:
                distance = tree_mat[sp][ref] 
                retention_rate = float(len(CNS_matrix[sp][CNS_matrix[sp] > 0].as_matrix()))/total_CNS
                plt.scatter(retention_rate,distance,s=3,label=sp)
                x_f.append(distance)
                y_f.append(retention_rate)
                if sp in included_species:
                    retention_rates[ref][sp] = retention_rate
                    y.append(retention_rate)
                    x.append(distance)
                    name.append(sp)
        plots.append(go.Scatter(x=x, y=y, name=ref,
                     mode='markers',
                     marker=dict(size=8,color=c[i]), text=name))
        print(x,y,name)
        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
        plt.xlabel('CNS Retention Rate')
        plt.ylabel('Genetic Distance')
        plt.title(reference_species+' CNS Retention Rate vs Genetic Distance, %s'%ref)
        plt.savefig('CNS_retention_vs_distance_%s.pdf'%ref,dpi=300)
    py.plot(go.Figure(data=plots,layout=go.Layout(title='Genetic Distance to Reference vs CNS Retention',yaxis=dict(title='Percent of CNS Elements Retained With Respect to Reference'),xaxis=dict(title='Genetic Distance from Reference'))),filename='CNS_retention_vs_distance.html',auto_open=False)
    plots = []
    for ref in refs:
        plots.append(go.Bar(x=retention_rates[ref].keys(),y=retention_rates[ref].values(),name=ref))
    py.plot(go.Figure(data=plots,layout=go.Layout({'barmode':'group'})),filename='CNS_retention.html',auto_open=False)
    x_f = np.nan_to_num(np.array(x_f))
    y_f = np.nan_to_num(np.array(y_f))
    print x_f,y_f
    p = np.polyfit(x_f,y_f,1)
    f = lambda x: np.polyval(p,x)
    expected = f(x_f)
    chi_sq = sum((y_f - expected)**2/expected)
    print('X^2',chi_sq,'p',1-stats.chi2.cdf(chi_sq,len(x_f)-2),'r',stats.pearsonr(x_f,y_f))
def CNS_RetentionvsGeneticDistance ( CNS_bed , reference_species , tree_file , species_file ) : import plotly . graph_objs as go import plotly . offline as py import matplotlib . pyplot as plt from collections import defaultdict from scipy import stats tree_mat = tree2matrix ( tree_file ) species = tree_mat . index . values print tree_mat print species included_species = os . popen ( 'cat ' + species_file ) . read ( ) . splitlines ( ) if not included_species : included_species = species refs = reference_species . split ( ',' ) c = [ 'hsl(' + str ( h ) + ',50%' + ',50%)' for h in np . linspace ( 0 , 360 , len ( set ( refs ) ) + 2 ) ] plots = [ ] retention_rates = { ref : defaultdict ( list ) for ref in refs } x_f , y_f = [ ] , [ ] for i , ref in enumerate ( refs ) : CNS_bed = CNS_bed [ : CNS_bed . rfind ( '/' ) + 1 ] + ref + '_AllCNSElements.bed' CNS_matrix = CNS_matrix_generation ( CNS_bed ) . drop_duplicates ( ) print CNS_matrix CNS_matrix = CNS_matrix . iloc [ : , 2 : ] print CNS_matrix total_CNS = CNS_matrix . shape [ 0 ] plt . figure ( ) x = [ ] y = [ ] name = [ ] for sp in species : if sp != reference_species : distance = tree_mat [ sp ] [ ref ] retention_rate = float ( len ( CNS_matrix [ sp ] [ CNS_matrix [ sp ] > 0 ] . as_matrix ( ) ) ) / total_CNS plt . scatter ( retention_rate , distance , s = 3 , label = sp ) x_f . append ( distance ) y_f . append ( retention_rate ) if sp in included_species : retention_rates [ ref ] [ sp ] = retention_rate y . append ( retention_rate ) x . append ( distance ) name . append ( sp ) plots . append ( go . Scatter ( x = x , y = y , name = ref , mode = 'markers' , marker = dict ( size = 8 , color = c [ i ] ) , text = name ) ) print ( x , y , name ) plt . legend ( bbox_to_anchor = ( 1.05 , 1 ) , loc = 2 , borderaxespad = 0. ) plt . xlabel ( 'CNS Retention Rate' ) plt . ylabel ( 'Genetic Distance' ) plt . title ( reference_species + ' CNS Retention Rate vs Genetic Distance, %s' % ref ) plt . savefig ( 'CNS_retention_vs_distance_%s.pdf' % ref , dpi = 300 ) py . plot ( go . Figure ( data = plots , layout = go . Layout ( title = 'Genetic Distance to Reference vs CNS Retention' , yaxis = dict ( title = 'Percent of CNS Elements Retained With Respect to Reference' ) , xaxis = dict ( title = 'Genetic Distance from Reference' ) ) ) , filename = 'CNS_retention_vs_distance.html' , auto_open = False ) plots = [ ] for ref in refs : plots . append ( go . Bar ( x = retention_rates [ ref ] . keys ( ) , y = retention_rates [ ref ] . values ( ) , name = ref ) ) py . plot ( go . Figure ( data = plots , layout = go . Layout ( { 'barmode' : 'group' } ) ) , filename = 'CNS_retention.html' , auto_open = False ) x_f = np . nan_to_num ( np . array ( x_f ) ) y_f = np . nan_to_num ( np . array ( y_f ) ) print x_f , y_f p = np . polyfit ( x_f , y_f , 1 ) f = lambda x : np . polyval ( p , x ) expected = f ( x_f ) chi_sq = sum ( ( y_f - expected ) ** 2 / expected ) print ( 'X^2' , chi_sq , 'p' , 1 - stats . chi2 . cdf ( chi_sq , len ( x_f ) - 2 ) , 'r' , stats . pearsonr ( x_f , y_f ) )
Reference should be CNS bed species
Reference should be CNS bed species
"""Reference should be CNS bed species""" #genetic_distances = tree_mat[reference_species] #FIXME binarize CNS matrix #FIXME debug #FIXME fix local PCA # FIXME finish, add tests for independence between genetic distance and retention, check if polyploid reduction in retention # FIXME Chi-square doesn't matter if you do not have any uncertainty about the y-data point
CNS_bed
reference_species
tree_file
species_file