<jupyter_start><jupyter_text>RNA Regulation through Localization- finding key features in subcellular RNA populations*****Purpose:** Restricting localization of RNA translation is one method that can regulate protein function within a cell. If an RNA is translated only in neurites, the protein function is also limited to the neurites unless transported away. This noteboook tries to determine key features of genes with transcripts translated in neurites (cell projections of neurons, axons/dendrites) over those that are sequestered to the cell body (soma). These features could predict how the cell determines where to send RNAs within a cell. **Data:** Source: Ouwenga et al 2019 (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6361723/)     1. List of ribosome bound RNAs enriched in the soma of cortical inhibitory neurons (VGAT expressing).      2. List of ribosome bound RNAs enriched in the neurites of cortical inhibitory neurons (VGAT expressing).**Pre-processing of data:** Extraction of gene/RNA/protein properties from Biomart Ensembl, Mouses Strains 102, Mouse C57BL/6NJ genes(C57BL_6NJ_v1)(http://uswest.ensembl.org/biomart/martview) **Input:** Dataframes of the following for soma enriched (v_soma.csv) and neurite enriched (v_neurite.csv) transcripts    1. Percent GC of the transcript     2. Number of different transcripts for the encoding gene    3. Which strand the RNA is created from        4. Chromosome     5. Whether or not the transcript encodes a transmembrane helix      **Output:**         1. Comparison of each feature between Neurite and Soma populations at multiple levels        Chromosome Level: Chr        Gene Level: Transcript count        RNA Level: PercentGC        Protein Level: TMHelix    2. Classification Tree    3. Random Forest Classification<jupyter_code>#importing packages
import pandas as pd
import numpy as np
import os 

import matplotlib as plt
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns 
plt.style.use('fivethirtyeight')
from IPython.display import Image


import scipy.stats as stats
from scipy.stats import ttest_ind 
from subprocess import check_output

import sklearn
from sklearn import tree
from sklearn.tree import DecisionTreeClassifier 
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

#Directory management 
os.getcwd()
os.chdir("C:\\Users\\ouwen\\Desktop\\Pythonclass")<jupyter_output><empty_output><jupyter_text>Data Import and Cleanup***<jupyter_code>#Defining subcellular location where each RNA population was harvested from 
#Enter Directory below
neuron = plt.imread("Neuron.jpg")
plt.imshow(neuron)

print('In neurons the cell body is called the soma and the dendrites and axons are collectively called neurites')
#Importing dataset produced by biomart from RNAs enriched in the neurites of cortical neurons
NeuriteTable = pd.read_csv("V_neurite.csv")
NeuriteTable.head()
#Importing dataset produced by biomart from RNAs enriched in the soma of cortical neurons
SomaTable = pd.read_csv("v_soma.csv")
SomaTable.head()
###Cleaning up the data 
#Changing the Headers to be easier to work with 
NeuriteTable = NeuriteTable.rename(columns={'Gene % GC content': 'PercentGC', 
                                            'Chromosome/scaffold name': 'Chr', 
                                            'Transmembrane helices': 'TMHelix'})

SomaTable = SomaTable.rename(columns={'Gene % GC content': 'PercentGC', 
                                            'Chromosome/scaffold name': 'Chr', 
                                            'Transmembrane helices': 'TMHelix'})

#Changing the NaN values in TMHelix to False:0 if it does not contain a TM Helix, True:1 if it does
NeuriteTable.TMHelix = NeuriteTable.TMHelix.map({'TMhelix': 1, np.nan: 0})
SomaTable.TMHelix = SomaTable.TMHelix.map({'TMhelix': 1, np.nan: 0})

#Checking one of the files to see the changes
NeuriteTable.head()
#SomaTable.head()
#Checking that all of the NaN values were changed and no other columns contain null data
print("Null Values in Somatic File\n", SomaTable.isnull().sum())
print("Null Values in Neurite File\n", NeuriteTable.isnull().sum())

#Checking the types of each column 
print(NeuriteTable.dtypes)
print("Note that Chromosome loaded correctly as a classifier, not a scalar value int64")
#Summaries of the data
print("Neurite Summary", "\n", NeuriteTable.describe(), "\n\n", 
      "Soma Summary", "\n", SomaTable.describe())
#Merging the data to one dataframe

#Giving each file a classier column Neurite: 1, Soma: 0
NeuriteTable['Class'] = 1
SomaTable['Class'] = 0

#merging
df= NeuriteTable.append(SomaTable)
df.head()<jupyter_output><empty_output><jupyter_text>*** Comparing Features between Soma and Neurite populations    Chromosome Level: Chr    Gene Level: Transcript count    RNA Level: PercentGC    Protein Level: TMHelix<jupyter_code>###CHROMOSOME LEVEL: Chromosome Location
#Which chromosomes encode the genes for these RNAs?

#Sorting the chromosomes to be in numerical order as objects
NChrList = list(NeuriteTable['Chr'])
NChrList.sort(key= lambda x: int(x) if x != 'X' else 23)
SChrList = list(SomaTable['Chr'])
SChrList.sort(key= lambda x: int(x) if x != 'X' else 23)

#plotting
plt.figure(figsize=(8,6))
sns.histplot(SChrList, label="Soma")
sns.histplot(NChrList, label="Neurite", color= "purple")
plt.title('Chromosome Location of Populations')
plt.xlabel('Chromosome')
plt.ylabel('Count')
plt.legend(prop={'size': 12})

#stats and summary
print("Given that chromosomes are named in descending size\n",
      "the appearance higher counts for chromosomes 1-11 is not unexpected\n",
      "An ANOVA could tell me if any are overrepresented, but it would be too underpowered to trust")
###GENE LEVEL: Transcript count
#Trancript Number for the genes that encode the RNAs located in neurites 
#An increased transcript number could mean location regulation through motif exclusion/inclusion 

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 9))

sns.histplot(NeuriteTable['Transcript count'], ax=ax1)
ax1.set_ylim([0, (len(df['Transcript count'])/5)])
ax1.set_xlim([1, df['Transcript count'].max()])
ax1.set_ylabel('Count')
ax1.set_title('Neurite Genes')

sns.histplot(SomaTable['Transcript count'], ax=ax2)
ax2.set_ylim([0, (len(df['Transcript count'])/5)])
ax2.set_xlim([1, df['Transcript count'].max()])
ax2.set_ylabel('')
ax2.set_title('Soma Genes')

sns.violinplot(x='Class', y="Transcript count", data=df, ax=ax3)
ax3.set_title('Distrubution of Counts')
ax3.set_xticklabels(['Soma', 'Neurite'])
ax3.set_xlabel('Population')

f.suptitle('Transcript Counts and Distribution Across Populations')

#stats
#This isn't normally distributed, but it will be generally work without further stats packages 
print(stats.ttest_ind(NeuriteTable['Transcript count'],
                SomaTable['Transcript count']))
#summary
print("The average number of transcripts of encoding genes is significantly greater for the neurite population")
             
###RNA SEQUENCE LEVE: GC Content
#Percent of the transcript that is GC 
#Higher GC content can create stable structures and binding 
f, (ax1, ax2) = plt.subplots(2, 1)

sns.boxplot(NeuriteTable['PercentGC'], ax=ax1)
ax1.set_ylabel('Neurite')
sns.boxplot(SomaTable['PercentGC'], ax=ax2)
ax2.set_ylabel('Soma')

f.suptitle('Percentage GC in Transcripts')

#stats
print(stats.ttest_ind(NeuriteTable['PercentGC'],
                SomaTable['PercentGC']))
#summary
print("Neurite RNAs have a significatly higher average GC content")



###PROTEIN LEVEL: Transmembrane Helix 
#A TMHelix is common in proteins involved with neurite functions

#Percentage of Neurite RNAs encoding a protein that includes a Transmembrane Helix
print('Neutite: ' + str(sum(NeuriteTable['TMHelix'])/len(NeuriteTable['TMHelix'])))

#Percentage of Soma RNAs encoding a Transmembrane Helix
print('Soma: ' + str(sum(SomaTable['TMHelix'])/len(SomaTable['TMHelix'])))

#stats
oddsratio, pvalue = stats.fisher_exact([[sum(NeuriteTable['TMHelix']), len(NeuriteTable['TMHelix']) - sum(NeuriteTable['TMHelix'])], [sum(NeuriteTable['TMHelix']), len(SomaTable['TMHelix']) - sum(SomaTable['TMHelix'])]])
print("pvalue:", pvalue)


#summary
print("Neurite RNAs are more likely to encode for a transmembrane helix")<jupyter_output>Neutite: 0.17634854771784234
Soma: 0.04407443682664055
pvalue: 3.2471475719141496e-21
Neurite RNAs are more likely to encode for a transmembrane helix<jupyter_text>*** Classification TreeIndividual features of the RNA populations differ, but can we use all of the features to attempt to classify them?<jupyter_code>#Splitting the data to testing and training data

testsize=0.25

X_train, X_test, Y_train, Y_test = train_test_split(df[['Strand', 'PercentGC', 'Transcript count', 'TMHelix']], 
                                                    df['Class'], random_state=0, test_size=testsize)
X_train.head()
#testing to see which depth gives us the best classification tree

max_depth_range = range(1,25)
accuracy = []
for depth in max_depth_range:
    
    clf = DecisionTreeClassifier(max_depth = depth, random_state = 0, criterion= "entropy")
    clf.fit(X_train, Y_train)
    score = clf.score(X_test, Y_test)
    accuracy.append(score)

print(accuracy)
print("\nMax accuracy:", max(accuracy), )
print("Ideal depth:", (accuracy.index(max(accuracy))+1))

#Plotting accuracy over depth to visualize the best depth to use
fig, ax = plt.subplots(figsize=[8,5])
ax.set_ylabel('Accuracy')
ax.set_xlabel('Depth')
ax.set_title('Accuracy over Depth')
ax.set_xticks(max_depth_range)
ax.grid(True)
plt.plot(max_depth_range, accuracy)
plt.show()
#Training the classification modeler using the depth determined above
clf = DecisionTreeClassifier(max_depth = (accuracy.index(max(accuracy))+1) , random_state = 0, criterion= "entropy")
clf.fit(X_train, Y_train)
score = clf.score(X_test, Y_test)
print(score)
#Running the test set against the classifier
predictions = clf.predict(X_test)

#Looking at the predictions
clf.predict(X_test[1:10])
#Plotting the tree
tree.plot_tree(clf) 
#Confusion Matrix
cm = metrics.confusion_matrix(Y_test, predictions)
plt.figure(figsize=(4,4))
sns.heatmap(cm, annot=True, fmt=".0f", linewidths=.5, square= True);
plt.ylabel("Actual");
plt.xlabel("Predicted");
plt.title('Accuracy Score: {}'.format(score), size = 12);

print("There are a lot of inaccurate classifications, this may be overfitted")
#null hypothesis test
MajorClass = max(Y_test.value_counts())
null = MajorClass / len(X_test)
print(null)
           
if null < (max(accuracy)):
    print("We reject the null!")
else:
    print("This could likely just be random")
importances_tree = pd.DataFrame({'feature':X_train.columns, 'importance':np.round(clf.feature_importances_,3)})
importances_tree = importances_tree.sort_values('importance', ascending=False).set_index('feature')

print(importances_tree)

print("\nGiven that there are so few have a TM Helix(a structure with high GC content) and strand was not significant by itself, \nit makes sense that the splits would rely more heavily on Percent GC and Transcript count to do prediction. \nSince features are correlated, perhaps a random forest would do a better job at predicting classifications")<jupyter_output>importance
feature                     
PercentGC              0.512
TMHelix                0.276
Transcript count       0.213
Strand                 0.000

Given that there are so few have a TM Helix(a structure with high GC content) and strand was not significant by itself, 
it makes sense that the splits would rely more heavily on Percent GC and Transcript count to do prediction. 
Since features are correlated, perhaps a random forest would do a better job at predicting classifications<jupyter_text>*** Random Forest<jupyter_code>#Declaring which columns to use for classifiers
feature_cols = ['Strand', 'PercentGC', 'Transcript count', 'TMHelix']
X = df[feature_cols]
y = df.Class
#splitting data
#adding rf for random forest to not reuse variables from the classification tree
#Using the same random state to better compare between trees
Xrf_train, Xrf_test, yrf_train, yrf_test = train_test_split(X,y, random_state = 0)
#determing the best depth 
max_depth_range = range(1,25)
accuracy = []

for depth in max_depth_range:
    
    clf_rf = RandomForestClassifier(n_estimators=150, max_depth= depth, random_state= 0)
    clf_rf.fit(Xrf_train, yrf_train)
    score = clf_rf.score(Xrf_test, yrf_test)
    accuracy.append(score)

print(accuracy)
print("\nMax accuracy:", max(accuracy), )
print("Ideal depth:", (accuracy.index(max(accuracy))+1))
#Visualizing accuracy over depth 
fig, ax = plt.subplots(figsize=[10,5])
ax.set_ylabel('Accuracy')
ax.set_xlabel('Depth')
ax.set_title('Accuracy over Depth')
ax.set_xticks(max_depth_range)
ax.grid(True)
plt.plot(max_depth_range, accuracy)
plt.show()
#Train the random forest 
clf_rf = RandomForestClassifier(n_estimators=150, max_depth= (accuracy.index(max(accuracy))+1) , random_state= 0)
clf_rf.fit(Xrf_train, yrf_train)
predictions_forest = clf_rf.predict(Xrf_test)
score = clf_rf.score(Xrf_test, yrf_test)
print(score)
# Compute feature importances.
importances_forest = pd.DataFrame({'feature':feature_cols, 'importance':clf_rf.feature_importances_}).sort_values(by='importance')
importances_forest = importances_forest.sort_values('importance', ascending=False).set_index('feature')

print("Importances from a simple decision tree:\n", importances_tree, "\n\n")
print("Importances from a random forest model:\n", importances_forest)

#Confusion Matrix

cmrf = metrics.confusion_matrix(yrf_test, predictions_forest)
plt.figure(figsize=(4,4))
sns.heatmap(cmrf, annot=True, fmt=".0f", linewidths=.5, square= True);
plt.ylabel("Actual");
plt.xlabel("Predicted");
plt.title('Accuracy Score: {}'.format(score), size = 12);

print("There are EVEN MORE inaccurate classifications- we must be overfitting")<jupyter_output>There are EVEN MORE inaccurate classifications- we must be overfitting
{'fraction_non_alphanumeric': 0.10390247320114902, 'fraction_numerical': 0.0145029075877531, 'mean_word_length': 5.1953125, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 18, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '293521', 'license': 'MIT', 'n_tokens_mistral': 4796, 'n_tokens_neox': 4391, 'n_words': 1406, 'path': 'Ouwenga_RNA Regulation through Localization.ipynb', 'repo_name': 'rebecca0nline/SynapticRNA_Regulation', 'size': '14273'}