CRMs_annot.py

# This module will graph and annotate the CRMs data.

# First do imports
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 os.path

# Define paths.
No_1 = './No_1Mo/'
PK = './Only_PKO/'
syn = './prism/syn/'

# WE will read in the file containing all CRMs as a dataframe. Rename the
# sample column to keep it consistent with other data.
df = pd.read_excel('CRMs.xlsx')
df = df.rename(columns={'sample' : 'Name'})

# Change the strain names from PKO and W402A to Polg-PKO and Polg-W402A as that
# is more accurate and is consistent with other data.
df = df.replace({'PKO' : 'Polg-PKO', 'W402A' : 'Polg-W402A'})

# Create a count column so that the CRMs can be counted.
df['count'] = 0
# Create a subset without W402A and create hue orders for graphing.
w402a = df[df['Strain'] != 'Polg-W402A']
Main_hue = ['Polg', 'Polg-PKO', 'Polg-W402A']
w402a_hue = ['Polg', 'Polg-PKO']

grouped = df.groupby(['Strain', 'Name'], observed=True).count().reset_index()
groupedw = w402a.groupby(['Strain', 'Name'], observed=True).count().reset_index()

# Get all unique strains for the three different conditions into the variables
# 'Strain', and 'Strain1'.
Strain = df.Strain.unique()
Strain1 = w402a.Strain.unique()

# Create an empty dataframe 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain:
        df3.append(grouped[grouped['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=grouped['count'],
                groups=grouped['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 as a csv  as
# 'Syn_all_stats.csv'
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(No_1,'CRMs_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(grouped, 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, 1,figsize=(5,5))
fig.suptitle('Number of CRMs', weight='bold')
ax = sns.violinplot(ax=axes,data=grouped,x=x, y=y, order=z,
                                facecolor=(1,1,1,0),edgecolor='.2')
grouped.to_excel(os.path.join(syn,'CRMs.xlsx'),index=False)

# 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, data=grouped, x=x, y=y,color='.2', order=z)
ax.set_ylim(top=max(grouped['count']) + 10)
ax.set_ylabel('Average Count')

# 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=grouped, x=x, y=y, order=z)
an.configure(text_format='star', loc='inside')
an.set_pvalues_and_annotate(p_values)

# Save the figure.
plt.savefig(os.path.join(No_1,'CRMs_all.png'), bbox_inches='tight')

# Create an empty dataframe 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain1:
        df3.append(groupedw[groupedw['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=groupedw['count'],
                groups=groupedw['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 as a csv  as
# 'Syn_all_stats.csv'
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(PK,'CRMs_w402a_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(groupedw, 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 = w402a_hue
sns.set_style('whitegrid', {'axes.grid' : False})
fig, axes = plt.subplots(1, 1,figsize=(5,5))
fig.suptitle('Number of CRMs', weight='bold')
ax = sns.barplot(ax=axes,data=groupedw,x=x, y=y, order=z,
                                facecolor=(1,1,1,0),edgecolor='.2')

# 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, data=groupedw, x=x, y=y,color='.2', order=z)
ax.set_ylim(top=max(groupedw['count']) + 10)
ax.set_ylabel('Average Count')

# 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=groupedw, x=x, y=y, order=z)
an.configure(text_format='star', loc='inside')
an.set_pvalues_and_annotate(p_values)

# Save the figure.
plt.savefig(os.path.join(PK,'CRMs_w402a.png'), bbox_inches='tight')

# We will read in the file again because we do not want the count coulmn as we
# will be looking at mean size of CRMs. Rename the
# sample column to keep it consistent with other data.
df = pd.read_excel('CRMs.xlsx')
df = df.rename(columns={'sample' : 'Name'})

# Change the strain names from PKO and W402A to Polg-PKO and Polg-W402A as that
# is more accurate and is consistent with other data.
df = df.replace({'PKO' : 'Polg-PKO', 'W402A' : 'Polg-W402A'})
# Create a subset without W402A and create hue orders for graphing.
w402a = df[df['Strain'] != 'Polg-W402A']
Main_hue = ['Polg', 'Polg-PKO', 'Polg-W402A']
w402a_hue = ['Polg', 'Polg-PKO']

grouped = df.groupby(['Strain', 'Name'])['final.size'].mean().reset_index()
groupedw = df.groupby(['Strain', 'Name'])['final.size'].mean().reset_index()

grouped.to_excel(os.path.join(syn,'CRM_size.xlsx'),index=False)

# Get all unique strains for the three different conditions into the variables
# 'Strain', and 'Strain1'.
Strain = df.Strain.unique()
Strain1 = w402a.Strain.unique()

# Create an empty dataframe 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain:
        df3.append(grouped[grouped['Strain'] == s]['final.size'])
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=grouped['final.size'],
                groups=grouped['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 as a csv  as
# 'Syn_all_stats.csv'
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(No_1,'CRMs_all_size_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(grouped, val_col='final.size', 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 = 'final.size'
z = Main_hue
sns.set_style('whitegrid', {'axes.grid' : False})
fig, axes = plt.subplots(1, 1,figsize=(5,5))
fig.suptitle('Size of CRMs', weight='bold')
ax = sns.boxplot(ax=axes,data=grouped,x=x, y=y, order=z,color='grey')
                              #  facecolor='grey',edgecolor='.2')#

# 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, data=grouped, x=x, y=y,color='.2', order=z)
ax.set_ylim(top=max(grouped['final.size']) + 10)
ax.set_ylabel('Average Size (NTs) of CRMs')

# 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=grouped, x=x, y=y, order=z)
#an.configure(text_format='star', loc='inside')
#an.set_pvalues_and_annotate(p_values)

# Save the figure.
plt.savefig(os.path.join(No_1,'CRMs_all_size.png'), bbox_inches='tight')

# Create an empty dataframe 'df3', and append Strain, Strain1,
# and Strain2 into the dataframe for the apporiate groups..
df3 = []
for s in Strain1:
        df3.append(groupedw[groupedw['Strain'] == s]['final.size'])
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=groupedw['final.size'],
                groups=groupedw['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 as a csv  as
# 'Syn_all_stats.csv'
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(PK,'CRMs_w402a_size_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(groupedw, val_col='final.size', 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 = 'final.size'
z = w402a_hue
sns.set_style('whitegrid', {'axes.grid' : False})
fig, axes = plt.subplots(1, 1,figsize=(5,5))
fig.suptitle('Size of CRMs', weight='bold')
ax = sns.barplot(ax=axes,data=groupedw,x=x, y=y, order=z,
                                facecolor=('grey'),edgecolor='.2')

# 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, data=groupedw, x=x, y=y,color='.2', order=z)
ax.set_ylim(top=max(groupedw['final.size']) + 10)
ax.set_ylabel('Average Size (NTs) of CRMs')

# 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=groupedw, x=x, y=y, order=z)
#an.configure(text_format='star', loc='inside')
#an.set_pvalues_and_annotate(p_values)

# Save the figure.
plt.savefig(os.path.join(PK,'CRMs_size_w402a.png'), bbox_inches='tight')

df = df.rename(columns={'final.end' : 'Mitochondrial Position', 'heteroplasmy' :
    'Heteroplasmy'})
w402a = df[df['Strain'] != 'Polg-W402A']
# A facetgrid is called. The aspect will be 2 to make it wider than taller. The
# setting 'legend_out is set to Tarue to allow room for the legend. 
X = sns.FacetGrid(data=df, col='Strain', col_wrap=1, despine=True,
                         hue='Strain', hue_order=Main_hue,
                                                  palette='deep',col_order=Main_hue,aspect=2)
# The graph title is afjusted to be on top of the graph, but not intruding on
# it, and it is emboldened.
X.fig.subplots_adjust(top=0.8)
X.fig.suptitle('End Position of CRM', fontsize=28,weight='bold')

# A histplot will be inserted into the facetgrid.  X is the Reference Position
# Y is SIFT-SCORE. We will graph the end of the deletion versus heteroplasmy.
# The start there is not much of a difference
X.map(sns.histplot, 'Mitochondrial Position',
                 'Heteroplasmy')

# Since facetgrids puts the title of each graph with '=', it is changes so
# that only the name without '=' is posted on the title of each graph.
X.set_titles(col_template='{col_name}',weight='bold')


# Save the figure as an SVG, with 600 dpi.
plt.savefig(os.path.join(No_1,'End_CRMs_all_histogram.svg'),dpi=600)

# A facetgrid is called. The aspect will be 2 to make it wider than taller. The
# setting 'legend_out is set to Tarue to allow room for the legend. 
X = sns.FacetGrid(data=w402a, col='Strain', col_wrap=1, despine=True,
                         hue='Strain', hue_order=w402a_hue,
                                                  palette='deep',col_order=w402a_hue,aspect=2)
# The graph title is afjusted to be on top of the graph, but not intruding on
# it, and it is emboldened.
X.fig.subplots_adjust(top=0.8)
X.fig.suptitle('End Position of CRM', fontsize=28,weight='bold')

# A histplot will be inserted into the facetgrid.  X is the Reference Position
# Y is SIFT-SCORE. We will graph the end of the deletion versus heteroplasmy.
# The start there is not much of a difference
X.map(sns.histplot, 'Mitochondrial Position',
                 'Heteroplasmy')

# Since facetgrids puts the title of each graph with '=', it is changes so
# that only the name without '=' is posted on the title of each graph.
X.set_titles(col_template='{col_name}',weight='bold')


# Save the figure as an SVG, with 600 dpi.
plt.savefig(os.path.join(PK,'End_CRMs_w402a_histogram.svg'),dpi=600)