<jupyter_start><jupyter_text>Investigating different sources of error- This notebook investigates the different sources of error in age inference. Besides the "bimodality" that we discuss at length in the manuscript, we also looked at whether protein length, number of domains, and evolutionary rate correlated with node error, our main error statistic- We also investigate whether protein length and evolutionary rate affect the age call, with shorter, more highly evolving proteins coming out younger than they should due to limited sensitivity in the inital search used to identify homologs<jupyter_code>import numpy as np
import pandas as pd
from scipy import stats
from matplotlib import pyplot as plt
from itertools import combinations
%matplotlib inline
## Suppress copious FutureWarnings from pandas

import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Read in the different files with error-source information<jupyter_code># Protein lengths are from Uniprot proteome. Number of domains was inferred with hmmer. Similarities here are derived 
# from pairwise alignments (needleman-wunsch) between one-to-one human and mouse orthologs, defined by PhylomeDB. 
# We look atHuman-Yeast pairs below, but chose to use Mouse here because there are more one-to-ones.

lengths = pd.read_csv("lengthMapping.csv",index_col=0,names=["protLength"]).sort_index()
numDomains = pd.read_csv("numDomains.csv",index_col=0).sort_index()
mouse_human = pd.read_csv("h-m_similarities.csv",index_col=0,names=["Hs-Mm_sim"]).sort_index()
nodeStats = pd.read_csv("../../nodeStats/nodeStats_HUMAN.csv",index_col=0).sort_index()[["Bimodality","NodeError"]]
lossTaxa = pd.read_csv("../../Errors/lossStats_HUMAN.csv",index_col=0).sort_index()[["mean"]]
lossTaxa.columns = ["meanLossTaxa"]
fracOverSplit = pd.read_csv("../../Errors/HUMAN_LDO_summary.csv",index_col=0).sort_index()
fracOverSplit.columns = ["fracOverSplit"]<jupyter_output><empty_output><jupyter_text>Make concatenated dataframe<jupyter_code>all_stats = pd.concat([lengths,numDomains,mouse_human,nodeStats,lossTaxa,fracOverSplit],axis=1)
all_stats.dropna(inplace=True)
all_stats = all_stats.convert_objects(convert_numeric=True)
all_stats.columns = ["Protein Length","# Domains","Hs-Mm_evoSim","Bimodality",
                     "Node Error","Avg. Loss Taxa","% Oversplit"]
all_stats.head()<jupyter_output><empty_output><jupyter_text>Scatterplot and correlation analysis> Only bimodality correlates strongly with the node error statistic<jupyter_code>fig, axes = plt.subplots(3,3,sharey=False)
index = 0
for col in all_stats.columns:
    axis = axes.flat[index]
    all_stats.plot(ax=axis,kind='scatter',x=col,y='Node Error',color='grey',title=col)
    axis.set_ylabel('')
    axis.set_xticklabels([])
    index += 1
    
#fig.savefig("errorScatterPlot.png",dpi=200)
# Get Spearman rank correlations between all variables

correlations = all_stats.corr("spearman")
correlations<jupyter_output><empty_output><jupyter_text>Plot absolute Spearman rank correlations with Node Error<jupyter_code>absCorrs = correlations.applymap(np.abs)

absCorrs.drop("Node Error",axis=1).ix["Node Error",].plot(
    kind='bar',color='grey',ylim=(0,1),title="Correlation with Node Error")
plt.ylabel("Absolute Spearman corr. coef.")

#plt.savefig("correlation_barPlot.svg")
# correlations.to_csv("errorCorrelations.csv")<jupyter_output><empty_output><jupyter_text>Human-yeast evolutionary rates and node errors- Now we'll take a quick look at yeast-human similarities. We used mouse above because there are many one-to-one orthologs, but the tradeoff is a small dynamic range. Human-yeast similarities have a stronger correlation with node error than human-mouse.<jupyter_code>yeast_human = pd.read_csv("h-y_similarities.csv",index_col=0,names=["Hs-Sc_sim"]).sort_index()
error_yeastSims = pd.concat([nodeStats,yeast_human],axis=1).dropna()[["Hs-Sc_sim","NodeError"]]
error_yeastSims = error_yeastSims.convert_objects(convert_numeric=True)
print "Number of alignments: %d" % len(error_yeastSims)
error_yeastSims.corr()
error_yeastSims.plot(kind='scatter',y="NodeError",x="Hs-Sc_sim",color='grey',title="Yeast-Human Distances Scatter")

#plt.savefig("Yeast-Human_dists_scatter.png",dpi=600)<jupyter_output><empty_output><jupyter_text>Effect of protein lengths and evolutionary distances on age-call.- Does the limited sensitivity of homology search algorithms affect the age-call? Prediction of Moyers and Zhang (2015,2016) is that quickly evolving proteins and short proteins will be called at a younger age incorrectly.- Takeaway from the below is that protein length and evolutionary rate may affect some of the proteins in the "Mammal" age class. So we cannot rule out that some of these are false negatives due to limited sensitivity in the initial homology search.  First we'll take a look at the distribution of human-mouse one-to-one alignment similarities. Most are highly similar.<jupyter_code>mouse_human.hist(bins=50,color='grey')<jupyter_output><empty_output><jupyter_text>Distribution of protein lengths. Most are short<jupyter_code>lengths.hist(bins=50,color='grey')
plt.yscale("log")<jupyter_output><empty_output><jupyter_text>Add in the mode-ages<jupyter_code>conAges = pd.read_csv("../../Main/main_HUMAN.csv",index_col=0)["modeAge"]
conAges.head()
statsAges = all_stats.join(conAges).sort("modeAge")
statsAges.head()
ages = ["Cellular_organisms","Euk+Bac","Euk_Archaea","Eukaryota",
        "Opisthokonta","Eumetazoa","Vertebrata","Mammalia"]

ageLengths = {}

for i in ages:
    df = statsAges[statsAges["modeAge"] == i]
    assert len(df) > 50 # just checkin'
    ageLengths[i] = df["Protein Length"].mean()<jupyter_output><empty_output><jupyter_text>Average protein lengths in each age category. No strong correlation with age, but Mammals may be enriched for small proteins.<jupyter_code>ax = pd.Series(ageLengths,index=ages).plot(color='black')
plt.xticks(rotation=70)
ax.set_ylabel("Average Protein Length")
ax.set_title("Average Protein Lengths in each Age Category")

#plt.savefig("lengthsXage.svg")<jupyter_output><empty_output><jupyter_text>Average similarity in each age class. Mammals definitely enriched for more quickly evolving proteins. Some of these may be older than "Mammal" and are false-negatives.<jupyter_code>ageSims = {}

for i in ages:
    df = statsAges[statsAges["modeAge"] == i]
    assert len(df) > 50 # just checkin'
    ageSims[i] = df["Hs-Mm_evoSim"].mean()
    
ax = pd.Series(ageSims,index=ages).plot(color='black')
plt.xticks(rotation=70)
ax.set_ylabel("Average Hs-Mm Similarity")
ax.set_title("Hs-Mm Evolutionary Similarity in each Age Category")

#plt.savefig("Hs-Mm_similarityXage.svg")<jupyter_output><empty_output><jupyter_text>Double-check that Mammal age class has significantly lower similarities in its proteins than Vertebrate class<jupyter_code># Create series for Mann-Whitney U-test

ageSimsSets = {}

for i in ages:
    df = statsAges[statsAges["modeAge"] == i]
    assert len(df) > 50 # just checkin'
    ageSimsSets[i] = df["Hs-Mm_evoSim"]

ageSimsSets["Mammalia"].head()
## Mann-Whitney U-test

f,p = stats.mannwhitneyu(ageSimsSets["Vertebrata"],ageSimsSets["Mammalia"])
print "p-val Mammals and Vertebrate classes have same mean: %s" % p
ageCatSimsDF = pd.DataFrame(ageSimsSets,columns=ages)<jupyter_output><empty_output><jupyter_text>Box plot of similarities<jupyter_code>ax = ageCatSimsDF.boxplot(sym=".")
plt.ylim(20,110)
plt.xticks(rotation=70)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11692635504061792, 'fraction_numerical': 0.006925023305366893, 'mean_word_length': 6.872117400419287, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '69524', 'license': 'MIT', 'n_tokens_mistral': 2579, 'n_tokens_neox': 2405, 'n_words': 598, 'path': 'Notebooks/errorCorrelation_plotting.ipynb', 'repo_name': 'bliebeskind/Gene-Ages', 'size': '7509'}