<jupyter_start><jupyter_text>Analysis of Variance of Read CountsTo support evaluation of classifiers, we need to know the variance (standard deviations) of read counts that are normalized in several ways: by library size, by gene size, and w.r.t. time 0. Two variances are of interest. The first is the variability of replications for the same gene and time. The second is the variability of mean values of replicas for a gene.<jupyter_code>import init
from common import constants as cn
from common.data_provider import DataProvider
from common_python.util import dataframe
from common_python.plots import util_plots
from common_python.plots.plotter import Plotter
from common_python.classifier.random_hypergrid_harness import RandomHypergridHarness

import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.distributions.empirical_distribution as stats

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Data<jupyter_code>provider = DataProvider()
provider.do()
# The data used in this analysis is a list of dataframe "matrix" that are the log2 values of the
# libary and gene adjusted readcounts divided by time 0.
# Each data frame has columns of genes and rows of time.
dfs_raw = provider.dfs_adjusted_read_count_wrt0
# Ignore the first timepoint since it's one because it's the reference.
dfs_nonref = [df.T[1:] for df in dfs_raw]
# Eliminate nan values
dfs_data = []
for df in dfs_nonref:
    dfs_data.append(df.applymap(lambda v: 0 if v == 0 else np.log(v)))
dfs_data[0].head()
dfs_data[1].head()
dfs_data[2].head()
len(dfs_data[0].columns)<jupyter_output><empty_output><jupyter_text>Standard Deviation Between Timepoints By Gene<jupyter_code>df_mean = sum(dfs_data)/len(dfs_data)
df_mean.head()
ser_stdb = df_mean.std() # Standard deviation between timepoints
ser_stdb.head()
ax = ser_stdb.hist(bins=100, range=(0, 2), density=False)

plotter = Plotter()
ecdf = stats.ECDF(ser_stdb)
plotter.ax.plot(ecdf.x, ecdf.y)
plotter.do(xlim=[0, 1.5], ylim=[0, 1.0], 
           title="STD Between Timepoints for %d genes" % len(df_mean.columns), 
           xlabel="std", ylabel="probability")<jupyter_output><empty_output><jupyter_text>Standard Deviation Within Gene, TimepointThis is the standard deviation across replications for a gene and timepoint. We use the average value across timepoints for each gene.<jupyter_code>df_stdw = dataframe.std(dfs_data)  # Dataframe of standard deviations within a timepoint
ser_stdw = df_stdw.mean()
#ser_stdw = pd.Series(np.reshape(df_stdw.to_numpy(), df_stdw.size))  # Convert to a series
ser_stdw.head()
plotter = Plotter()
ecdf = stats.ECDF(ser_stdw)
plotter.ax.plot(ecdf.x, ecdf.y)
plotter.do(xlim=[0, 1], ylim=[0, 1.0], title="STD Within Timepoints for %d genes" % len(df_stdw.columns), 
           xlabel="std", ylabel="probability")<jupyter_output><empty_output><jupyter_text>Analysis of Top GenesThese are genes identified by SVM classifier as having the most effect on the classification vector as indicated by the plot below.![](pngs/top10genes.png)![](pngs/top15genes.png)<jupyter_code>STDB = "stdb"
STDW = "stdw"
top_genes10 = ["Rv2382c", "Rv2329c", "Rv3083", "Rv2958c", "Rv3613c", "Rv2385", "Rv3841", "Rv3342", "Rv0040c"]
top_genes15 = ["Rv0886", "Rv2382c", "Rv3083", "Rv3613c", "Rv2329c", "Rv3136", "Rv2958c", "Rv3841", "Rv2377c",
            "Rv2380c", "Rv3342", "Rv1479", "Rv1700", "Rv2385"]
gene_dct = {k: {STDB: None, STDW: None} for k in top_genes15}
for gene in top_genes15:
    dct = gene_dct[gene]
    dct[STDB] = ser_stdb[gene]
    dct[STDW] = ser_stdw[gene]
df = pd.DataFrame(gene_dct)
df.transpose()
dff = df.transpose()
dff[STDB] / dff[STDW]
ecdf = stats.ECDF(dff[STDB])
plotter = Plotter()
plotter.ax.plot(ecdf.x, ecdf.y)
plotter.do(xlim=[0, 1], ylim=[0, 1.0], title="CDF of STD Between for 15 top Genes",
          xlabel="STDB", ylabel="Probability")
ecdf = stats.ECDF(dff[STDW])
plotter = Plotter()
plotter.ax.plot(ecdf.x, ecdf.y)
plotter.do(xlim=[0, 1], ylim=[0, 1.0], title="CDF of STD Within for 15 top Genes",
          xlabel="STDW", ylabel="Probability")
len(set(top_genes10).intersection(top_genes15))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11959654178674352, 'fraction_numerical': 0.035782901056676274, 'mean_word_length': 5.611111111111111, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '447656', 'license': 'MIT', 'n_tokens_mistral': 1651, 'n_tokens_neox': 1514, 'n_words': 415, 'path': 'notebooks/Variance of Read Counts.ipynb', 'repo_name': 'uwescience/expression_state', 'size': '4164'}