corrections2.py

# This script will make further corrections on the "paper_figs_clc.py" script.

# First, we will import the appropriate modules.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from statannotations.Annotator import Annotator as An
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from scikit_posthocs import posthoc_tukey
import scipy.stats as stats
from scipy.stats import f_oneway
from io import StringIO
import warnings
import os.path
import statsmodels.api as sm
from statsmodels.formula.api import ols
import statsmodels.api as sm
from statsmodels.formula.api import ols

# Define paths.
All = './All/'
No_1 = './No_1Mo/'
PK = './Only_PKO/'


# First we will read in the file containing all CLC data as a dataframe 'df'
# and create a count column so that things can be counted and stats can be
# done..
df = pd.read_excel('Concat_with_all.xlsx')
df['count'] =0

df = df[df['Type'] != 'Replacement']

# The Frequency column will be renamed to Heteroplasmy and Reference Position
# to Mitochondrial Position as these are what these columns represent.
df = df.rename(columns ={'Frequency' : 'Heteroplasmy', 'Reference Position' :
    'Mitochondrial Position'})

# Now we will create subsets - first without 1 month old animals, and then
# without W402A animals.
no1mo = df[df['Strain'] != 'WT\n1Mo']
no1mo = no1mo[no1mo['Strain'] != 'PKO\n1Mo']
w402a = no1mo[no1mo['Strain'] != 'Polg-W402A']

# We will replace the 'PKO\n1Mo' to 'Polg-PKO\n1Mo' as this is a more accurate
# description of the strain.
df = df.replace({'PKO\n1Mo' : 'Polg-PKO\n1Mo'})

# Now we will subset sample_type. Use the '.copy()' because otherwise there is
# a warning downstream when the count column is created because without the
# .copy(), it's a view versus a copy - may check docs.
Syn = df[df['Sample_type'] == 'Syn'].copy()
WBH = df[df['Sample_type'] == 'Homogenate'].copy()
no1moSyn = no1mo[no1mo['Sample_type'] == 'Syn'].copy()
no1moWBH = no1mo[no1mo['Sample_type'] == 'Homogenate'].copy()
w402aSyn = w402a[w402a['Sample_type'] == 'Syn'].copy()
w402aWBH = w402a[w402a['Sample_type'] == 'Homogenate'].copy()

# We are going to create hue orders for the different dataframes. These hue
# orders will be used in the graphs.
Main_hue = ['WT', 'WT\n1Mo', 'Polg', 'Polg-PKO', 'Polg-PKO\n1Mo', 'Polg-W402A']
no1mo_hue = ['WT', 'Polg', 'Polg-PKO', 'Polg-W402A']
w402a_hue = ['WT', 'Polg', 'Polg-PKO']

# For each of the dataframes above, group by 'Strain', 'Name', and aggregate as
# count. This will add the counts into the 'count' column. In addition, group
# by 'Strain, 'Name', and 'Type', and aggregate as count. Finally, the index
# needs to be reset so that the variables are available in the downstream
# graphs.
Syngrp = Syn.groupby(['Strain', 'Name'], observed=True).count().reset_index()
Syngrp1 = Syn.groupby(['Strain', 'Name', 'Type'], observed=True).count().reset_index()
WBHgrp = WBH.groupby(['Strain', 'Name'], observed=True).count().reset_index()
WBHgrp1 = WBH.groupby(['Strain', 'Name','Type'], observed=True).count().reset_index()
no1Syngrp = no1moSyn.groupby(['Name', 'Strain'], observed=True).count().reset_index()
no1Syngrp1 = no1moSyn.groupby(['Strain', 'Name', 'Type'], observed=True).count().reset_index()
no1WBHgrp = no1moWBH.groupby(['Strain', 'Name'], observed=True).count().reset_index()
no1WBHgrp1 = no1moWBH.groupby(['Strain', 'Name', 'Type'], observed=True).count().reset_index()
w402aSyngrp = w402aSyn.groupby(['Strain', 'Name']).count().reset_index()
w402aSyngrp1 = w402aSyn.groupby(['Name', 'Strain', 'Type'], observed=True).count().reset_index()
w402awbhgrp = w402aWBH.groupby(['Name', 'Strain'], observed=True).count().reset_index()
w402awbhgrp1 = w402aWBH.groupby(['Strain', 'Name', 'Type'], observed=True).count().reset_index()

# Get all unique strains for the three different conditions into the variables
# 'Strain', 'Strain1', and 'Strain2'.
Strain = Syngrp.Strain.unique()
Strain1 = no1Syngrp.Strain.unique()
Strain2 = w402aSyngrp.Strain.unique()

# Create an empty dataframes 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain:
    df3.append(Syngrp[Syngrp['Strain'] == s]['count'])
F = f_oneway(*df3)    

# Cast the f_oneway statistics into a string so that it can then be saved as a
# dataframe. Use StringIO to implement a file-like class on the string ('F') so
# that it can be read in as a CSV and hence become a dataframe.
F = str(F)
F = StringIO(F)
F = pd.read_csv(F, sep=";")

# Perform groupwise comparisons using tukey HSD.
tukey = pairwise_tukeyhsd(endog=Syngrp['count'],
        groups=Syngrp['Strain'],
        alpha=0.05)

# Save the tukey data as a dataframe and append the F stats from the F
# dataframe into it. Finally, export the dataframe to excel as
# 'Syn_all_stats.xlsx'
tukey = pd.DataFrame(data=tukey._results_table.data[1:],
        columns=tukey._results_table.data[0])
tukey = pd.concat([tukey, F])
tukey.to_csv(os.path.join(All,'Syn_all_stats.csv'),index=False)

# In order to annonate the resulting graphs with stars for statistical
# significance, the tukey results need to be put into a simple matrix (just
# showing the p values, not the other parameters such as meandiff, lower and
# upper etc).
tukey_df = posthoc_tukey(Syngrp, val_col='count', group_col='Strain')

# The matrix needs to be converted to a non-redundant list of comparisons with
# the p-value. This is done by removing the lower half and diagonal of the
# matrix and turning the matrix format into a long dataframe using melt(). The
# code and resulting dataframe are shown below.
remove = np.tril(np.ones(tukey_df.shape), k=0).astype('bool')
tukey_df[remove] = np.nan
molten_df = tukey_df.melt(ignore_index=False).reset_index().dropna()

# Get all unique strains for the three different conditions into the variables
# 'Strain', 'Strain1', and 'Strain2'.
Strain = Syngrp.Strain.unique()
Strain1 = no1Syngrp.Strain.unique()
Strain2 = w402aSyngrp.Strain.unique()

# Create an empty dataframes 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain:
    df3.append(Syngrp[Syngrp['Strain'] == s]['count'])
F = f_oneway(*df3)    

# Cast the f_oneway statistics into a string so that it can then be saved as a
# dataframe. Use StringIO to implement a file-like class on the string ('F') so
# that it can be read in as a CSV and hence become a dataframe.
F = str(F)
F = StringIO(F)
F = pd.read_csv(F, sep=";")

# Perform groupwise comparisons using tukey HSD.
tukey = pairwise_tukeyhsd(endog=Syngrp['count'],
        groups=Syngrp['Strain'],
        alpha=0.05)

# Save the tukey data as a dataframe and append the F stats from the F
# dataframe into it. Finally, export the dataframe to excel as
# 'Syn_all_stats.xlsx'
tukey = pd.DataFrame(data=tukey._results_table.data[1:],
        columns=tukey._results_table.data[0])
tukey = pd.concat([tukey, F])
tukey.to_csv(os.path.join(All,'Syn_all_stats.csv'),index=False)

# In order to annonate the resulting graphs with stars for statistical
# significance, the tukey results need to be put into a simple matrix (just
# showing the p values, not the other parameters such as meandiff, lower and
# upper etc).
tukey_df = posthoc_tukey(Syngrp, val_col='count', group_col='Strain')

# The matrix needs to be converted to a non-redundant list of comparisons with
# the p-value. This is done by removing the lower half and diagonal of the
# matrix and turning the matrix format into a long dataframe using melt(). The
# code and resulting dataframe are shown below.
remove = np.tril(np.ones(tukey_df.shape), k=0).astype('bool')
tukey_df[remove] = np.nan
molten_df = tukey_df.melt(ignore_index=False).reset_index().dropna()

# x, y, and order are defined.
x = 'Strain'
y = 'count'
z = Main_hue
sns.set_style('whitegrid', {'axes.grid' : False})
fig, axes = plt.subplots(1, 2,figsize=(17,4.5))
fig.suptitle('Synaptosomes', weight='bold')
ax = sns.violinplot(ax=axes[0],data=Syngrp,x=x, y=y, order=z,color='white')
                       # facecolor=(1,1,1,0),edgecolor='.2')
ax.set_title('All Mutations')

# Add a swarmplot to visualize individual datapoints on the barplot. Color the
# swarmplot black so that all points are visible.
ax = sns.swarmplot(ax=axes[0], data=Syngrp, x=x, y=y,color='.2', order=z)
ax.set_ylim(top=max(Syngrp['count']) + 10)
ax.set_ylabel('Number of Mutations')

# In order to annonate the graph where there are significant differences, the
# dataframe 'molten_df', which contains the p-values, will be filtered so that
# only significant p values (<=0.05) are in there. Note, if all notations are
# desirable, skip the filtering step.
molten_df = molten_df.loc[molten_df['value'] <= 0.05]

# The pairs for multiple comparisions is defined as all strains in the p value
# table.
pairs = [(i[1]['index'], i[1]['variable']) for i in molten_df.iterrows()]

# A list of p values is generated from the molten_df dataframe. The annonator
# is then defined and configured to annonate the graph with stars using the p
# values from the list.
p_values = [i[1]['value'] for i in molten_df.iterrows()]
an = An(ax, pairs, data=Syngrp, x=x, y=y, order=z)
an.configure(text_format='star', loc='inside')
an.set_pvalues_and_annotate(p_values)

# Make the graph for mutations separated by type.
ax2 = sns.swarmplot(ax=axes[1], data=Syngrp1, x=x, y=y, hue='Type',
        order=z)
s = dict(x=x, y=y, data=Syngrp1)

#perform two-way ANOVA
model = ols('count ~ C(Strain) + C(Type) + C(Strain):C(Type)', data=Syngrp1).fit()
sm.stats.anova_lm(model, typ=2)
print(type(sm.stats.anova_lm(model, typ=2)))

# Perform groupwise comparisons using tukey HSD.
tukey1 = pairwise_tukeyhsd(endog=Syngrp1['count'],
        groups=Syngrp1['Strain'],
        alpha=0.05)
# Save the tukey data as a dataframe and append the F stats from the F
# dataframe into it. Finally, export the dataframe to excel as
# 'Syn_all_stats.xlsx'
tukey1 = pd.DataFrame(data=tukey1._results_table.data[1:],
        columns=tukey1._results_table.data[0])
#tukey1 = pd.concat([tukey, F])
tukey1.to_csv(os.path.join(All,'Syn_all_stats_type.csv'),index=False)

# In order to annonate the resulting graphs with stars for statistical
# significance, the tukey results need to be put into a simple matrix (just
# showing the p values, not the other parameters such as meandiff, lower and
# upper etc).
tukey1_df = posthoc_tukey(Syngrp1, val_col='count', group_col='Strain')

# The matrix needs to be converted to a non-redundant list of comparisons with
# the p-value. This is done by removing the lower half and diagonal of the
# matrix and turning the matrix format into a long dataframe using melt(). The
# code and resulting dataframe are shown below.
remove = np.tril(np.ones(tukey1_df.shape), k=0).astype('bool')
tukey1_df[remove] = np.nan
molten_df = tukey1_df.melt(ignore_index=False).reset_index().dropna()
# In order to annonate the graph where there are significant differences, the
# dataframe 'molten_df', which contains the p-values, will be filtered so that
# only significant p values (<=0.05) are in there. Note, if all notations are
# desirable, skip the filtering step.
molten_df = molten_df.loc[molten_df['value'] <= 0.05]

# The pairs for multiple comparisions is defined as all strains in the p value
# table.
pairs = [(i[1]['index'], i[1]['variable']) for i in molten_df.iterrows()]

# A list of p values is generated from the molten_df dataframe. The annonator
# is then defined and configured to annonate the graph with stars using the p
# values from the list.
p_values = [i[1]['value'] for i in molten_df.iterrows()]
an = An(ax2, pairs, data=Syngrp1, x=x, y=y, order=z)
an.configure(text_format='star', loc='inside')
an.set_pvalues_and_annotate(p_values)
sns.pointplot(**s,join=False, ci=0,capsize=.7,scale=0, order=z, hue='Type')

# The following was added because otherwise the legend comes up twice.
handles, labels=ax2.get_legend_handles_labels()
ax2.legend(handles[:3],labels[:3])
ax2.set_title('Mutations by Type')
ax2.set_ylabel('Number of Mutations')
sns.move_legend(ax2, 'upper left', bbox_to_anchor=(1, 1.025))
plt.savefig(os.path.join(All,'Syn_all_Mut_Num.png'), bbox_inches='tight')