morar¶

import os
import pandas as pd
import numpy as np
import morar
import matplotlib.pyplot as plt
%matplotlib inline

# load in example single cell data
data = pd.read_csv("/home/scott/Dropbox/single_cell_example_data.csv")
n_rows, n_cols = data.shape
print("{} rows \n{} columns".format(n_rows, n_cols))
data.head()

data.featurecols[:10] # first 10

data.featuredata.head()

data.metacols

data.agg(on="Metadata_platename")

data.agg(on="Metadata_platename", method="mean")

# function to plot nuclei eccentricity per plate
def norm_plot(dat, title):
    names, vals = [], []
    grouping = dat.groupby("Metadata_platename")
    for name, group in grouping:
        tmp = group.Nuclei_AreaShape_Eccentricity.values.tolist()
        vals.append(tmp)
        names.append(name)

    plt.boxplot(vals)
    plt.xticks(range(1, len(names)+1), names)
    plt.xlabel("plate name")
    plt.title(title, fontweight="bold")
    


plt.figure(figsize=[8,6])

plt.subplot(221)
norm_plot(data, "Un-normalised data")
plt.ylabel("nuclei eccentricity")

plt.subplot(222)

# normalise data per plate (cell-line)
data_norm = data.normalise(plate_id="Metadata_platename")

norm_plot(data_norm, "Normalised data")

plt.tight_layout()
plt.show()

# function to make plots
def scale_plot(dat, title):
    """make a plot of cell-area and nuclei-eccentricity"""
    plt.figure(figsize=[8, 4])
    plt.subplot(221)
    plt.title("Cell Area")
    plt.hist(dat["Cells_AreaShape_Area"], bins=100, normed=1)

    plt.subplot(222)
    plt.title("Nuclei Eccentricity")
    plt.hist(dat["Nuclei_AreaShape_Eccentricity"], bins=100, normed=1)
    plt.tight_layout()
    plt.subplots_adjust(top=0.85)
    plt.suptitle(title, fontweight="bold")
    plt.show()
    

    
scale_plot(data, "Unscaled data")

data_scaled = data.scale_features()

scale_plot(data_scaled, "Scaled data")

area_mean = data_scaled.Cells_AreaShape_Area.mean()
area_std = data_scaled.Cells_AreaShape_Area.std()
print("mean = {:.2f}\n\
standard deviation = {:.2f}".format(area_mean, area_std))

pca_df.head()

# plot the calculated principal components
plt.figure(figsize=[12, 9])
plt.suptitle("Plot of first 2 principal components", fontweight="bold")
plt.subplot(221)
plt.scatter(pca_df.PC1, pca_df.PC2, s=3, alpha=0.3)
plt.xlim([-10, 25]); plt.ylim([-10, 20])
plt.subplot(222)
for name, group in pca_df.groupby("Metadata_platename"):
    plt.scatter(group.PC1, group.PC2, label=name, s=3, alpha=0.5)
plt.xlim([-10, 25]); plt.ylim([-10, 20])
plt.legend()
plt.subplots_adjust(top=0.95)
plt.show()

plt.figure(figsize=[10, 6])
plt.subplot(221)
plt.plot(pca_var, ".")
plt.vlines(x=range(len(pca_var)), ymin=0, ymax=pca_var)
plt.ylabel("proportion of variance")
plt.xlabel("principal component")
plt.title("Proportion of variance explained \nby each principal component")

plt.subplot(222)
cumulative_var = np.cumsum(pca_var) / np.sum(pca_var)
plt.plot(cumulative_var)
plt.title("Cumulative variation explained by \nsuccessive principal components")
plt.xlabel("principal component")
plt.ylabel("cumulative proportion of variance")
plt.tight_layout()
plt.show()

wanted_proportion = 0.8
# returns the first index of where `cumulative_var` exceeds or equals the wanted_proportion
# as False == 1
np.argmax(cumulative_var >= wanted_proportion) + 1 # as zero indexed

pca_df2, pca_var2 = data_scaled.pca(n_components=0.8)
pca_df2.columns

from morar import feature_selection

# find one of pairs of features that have a correlation > 0.95
features_to_rm = feature_selection.find_correlation(data_scaled, threshold=0.95)

features_to_rm