<jupyter_start><jupyter_text>Means sigmas<jupyter_code>df_mv = pd.read_csv("meanVariances.csv", index_col = [0])
#type_of_gene='protein-coding'
#df_g = pd.read_csv('/Users/filippo/Developer/tesi/genes.txt', header=0, index_col=0)
#df_mv=df_mv[df_mv.index.isin(df_g[df_g['type_of_gene']=='protein-coding'].index)]
#df_mv = df_mv.loc[df_mv['type_of_gene']==type_of_gene]
if not df_mv.columns.isin(['occurrence']).any():
    df_mv_occ=pd.read_csv("O.dat", header=None)
    df_mv.insert(3, 'occurrence', df_mv_occ.values)
df_mv.dropna(axis=0,how='any',inplace=True)
df_mv.head()
#df_mv.round(2).to_csv("meanVariances.csv",index=True,header=True)
#df_mv.fillna(value=0.,inplace=True)
means = df_mv['mean'].values
variances = df_mv['variance'].values
occurrences = np.array(df_mv['occurrence'].values, dtype=float)
abundances = pd.read_csv("A.dat", header=None).values
len(df_mv)<jupyter_output><empty_output><jupyter_text>plot  **var** versus **mean**<jupyter_code>fig=plt.figure(figsize=(15,8))
ax=fig.subplots()
plotvarmen(means, variances, ax=ax, normalisation_str=normalisation_str)
fig.savefig("varmean_loglog.png")
fig=plt.figure(figsize=(15,8))
ax=fig.subplots()
plotcv2mean(means, variances, ax=ax, normalisation_str=normalisation_str)
fig.savefig("cvmean_loglog.png")
import scanpy as sc
adata = sc.AnnData(df.transpose(),obs=df_files.reindex(index=df.transpose().index), var=df_mv)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=1000, n_bins=100)
x = adata.var['mean'].values
y = adata.var['variance'].values
cv2 = [yi/xi/xi if xi!=0 else 0 for xi, yi in zip(x,y)]
y=np.array(cv2)
HVx = adata.var[adata.var['highly_variable']==True]['mean'].values
HVy = adata.var[adata.var['highly_variable']==True]['variance'].values
HVcv2 = [yi/xi/xi if xi!=0 else 0 for xi, yi in zip(HVx,HVy)]
HVy=np.array(HVcv2)
fig=plt.figure(figsize=(15,8))
plt.scatter(x,y)
plt.scatter(HVx,HVy)
plt.ylim((y[y!=0].min()*0.95,y.max()*1.5))
plt.xlim((x[x!=0].min()*0.95,x.max()*1.5))
plt.xscale('log')
plt.yscale('log')
plt.show()
x = adata.var['means'].values
y = adata.var['dispersions'].values
HVx = adata.var[adata.var['highly_variable']==True]['means'].values
HVy = adata.var[adata.var['highly_variable']==True]['dispersions'].values
fig=plt.figure(figsize=(20,10))
plt.scatter(x,y, label='genes')
plt.scatter(HVx,HVy, label='highly variable genes')
plt.ylim((np.nanmin(y)-0.1,np.nanmax(y)))
plt.xlim((-0.1,x.max()))
plt.xlabel("Log - Mean", fontsize=18)
plt.ylabel("Log - Dispersion", fontsize=18)
plt.legend(fontsize=20)
plt.show()
fig.savefig("DispersionMean.pdf")<jupyter_output><empty_output><jupyter_text>overexpressed<jupyter_code>how_many_sigmas = 1
distance = 10
fig=plt.figure(figsize=(15,8))
ax=fig.subplots()
plotoversigmacv2(means, variances, ax=ax, normalisation_str=normalisation_str, how_many_sigmas=how_many_sigmas)
fig.savefig("cvmean_loglog_%dsigma.png"%how_many_sigmas)
def get_mean_cv2(mean, cv2, knee=1., distance=distance):
    if mean < knee:
        return(mean, cv2, -1, -1, cv2 > distance+1./mean)
    else:
        return(mean, cv2, -1, -1, cv2 > 1e-1+distance)
def get_mean_cv2(mean, cv2, bin_means=[], bin_sigmas=[],bin_edges=[],how_many_sigmas=3):
    bin_i = 0
    for i in range(len(bin_edges[:-1])):
        if mean<=bin_edges[i+1] and mean > bin_edges[i]:
            bin_i = i
            break
    return(mean, cv2, bin_means[bin_i], bin_sigmas[bin_i], cv2>(bin_means[bin_i]+how_many_sigmas*bin_sigmas[bin_i]))
over, over_plot = getovergenes(df_mv,get_mean_cv2, how_many_sigmas=how_many_sigmas, method='sigma')
fig=plt.figure(figsize=(15,8))
ax = fig.subplots()
plotoverpoints(means, variances, over_plot, ax=ax, how_many_sigmas=how_many_sigmas)
fig.savefig("cvmean_loglog_oversigma.png")
#discretize_df_columns(df.loc[over, df.columns[np.unique(np.random.randint(len(df.columns), size=2000))]]).to_csv("main_Table_Altmann.csv",index=True, header=True)
df[df.index.isin(over)].dropna().to_csv("mainTable_over.csv",index=True, header=True)<jupyter_output><empty_output><jupyter_text>mean versus occurrence<jupyter_code>fig=plt.figure(figsize=(8,5))
plt.scatter(occurrences*nfiles, means, c='b', alpha=0.8, label='genes')
if 'counts' in normalisation_str:
    plt.plot(np.linspace(1,nfiles), np.linspace(1,nfiles)/(nfiles), lw=4, label='bound', c='cyan', ls='--')
bin_means, bin_edges, _ = st.binned_statistic(occurrences*nfiles, means, statistic='mean', bins=np.logspace(-3,6))
x = (bin_edges[1:]+bin_edges[:-1])/2
plt.scatter(x,bin_means, marker='x', c='r', label='binned average')
plt.ylabel("$<%s>$"%normalisation_str, fontsize=16)
#plt.xlabel("$\Sigma_j\Theta(FPKM-0.1)\Theta(1e5-FPKM)$", fontsize=16)
plt.xlabel("$\Sigma_j\Theta(%s)$"%normalisation_str, fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.ylim(means[means.nonzero()].min()/5,np.power(10,np.log10(means.max())+1))
plt.xlim(5e-1,nfiles+800)
plt.legend(fontsize=20)
plt.show()
fig.savefig("meanDiff_loglog.png")
fig=plt.figure(figsize=(8,5))
plt.scatter(means, occurrences*nfiles, c='b', alpha=0.6, label='data')
bin_means, bin_edges, _ = st.binned_statistic(means, occurrences*nfiles, statistic='mean', bins=np.logspace(-3,6))
x = (bin_edges[1:]+bin_edges[:-1])/2
plt.scatter(x,bin_means, marker='x', c='r', label='binned average')
plt.xlabel("$<%s>$"%normalisation_str, fontsize=16)
#plt.xlabel("$\Sigma_j\Theta(FPKM-0.1)\Theta(1e5-FPKM)$", fontsize=16)
plt.ylabel("$\Sigma_j\Theta(%s)$"%normalisation_str, fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.xlim(means[means.nonzero()].min()/5,np.power(10,np.log10(means.max())+1))
plt.ylim(5e-1,nfiles+800)
plt.legend(fontsize=20)
plt.show()
fig.savefig("diffMean_loglog.png")<jupyter_output><empty_output><jupyter_text>Abundance vs occurrence<jupyter_code>fig=plt.figure()
plt.scatter(occurrences, abundances.T[0])
plt.xscale('log')
plt.yscale('log')
plt.xlim((1e-5,1))
plt.ylim((np.min(abundances.T[0][abundances.T[0].nonzero()]),np.max(abundances.T[0])))
plt.show()
plt.scatter(np.divide(abundances,occurrences), means)<jupyter_output><empty_output><jupyter_text>Distributions<jupyter_code>fig = plt.figure()
data = means
mu = np.median(data)
s = np.std(data)
log_bins = np.logspace(-5,7)
plt.hist(data, density = False, histtype='step', bins=log_bins)
plt.title("means", fontsize=16)
plt.xlabel("$<%s>$"%normalisation_str, fontsize=16)
plt.ylabel("#", fontsize=16)
plt.yscale('log')
plt.xscale('log')
plt.xlim(5e-5,5e5)
plt.show()
fig.savefig("mean_distr.pdf")
fig = plt.figure(figsize=(7,7))
data = -np.log10(means[means.nonzero()])
minimum = data.min()
data -= minimum
mu = np.median(data)
s = np.std(data)
fit_params = st.lognorm.fit(data)
n, c, _ = plt.hist(data, density = True, histtype='step')
plt.plot(np.linspace(-1,10),st.lognorm.pdf(np.linspace(-1,10),*fit_params), label='lognormal (s:%.2f,loc:%.2f,scale:%.2f)'%(fit_params))
plt.title("means", fontsize=16)
plt.xlabel("$-log10(<%s>)+%.2f$"%(normalisation_str,-minimum), fontsize=16)
plt.ylabel("pdf", fontsize=16)
plt.yscale('log')
plt.ylim(1e-4,1)
plt.xlim(-1,10)
plt.legend(fontsize=16)
plt.show()
fig.savefig("mean_distr_scaled.pdf")
bins = 80
bins=np.logspace(np.log10(5e-1),np.log10(variances.max()))
fig = plt.figure()
n, c, _ = plt.hist(variances, density = False, histtype='step', bins=bins)

plt.title("vars", fontsize=16)
plt.xlabel("$<\sigma_{%s}^2>$"%normalisation_str, fontsize=16)
plt.ylabel("#", fontsize=16)
plt.yscale('log')
plt.xscale('log')
plt.show()
fig.savefig("var_distr.pdf")
fig = plt.figure(figsize=(7,7))
data = -np.log10(variances[variances.nonzero()])
minimum = data.min()
data -= minimum
mu = np.median(data)
s = np.std(data)
fit_params = st.lognorm.fit(data)
n, c, _ = plt.hist(data, density = True, histtype='step')
plt.plot(np.linspace(-1,data.max()),st.lognorm.pdf(np.linspace(-1,data.max()),*fit_params), label='lognormal (s:%.2f,loc:%.2f,scale:%.2f)'%(fit_params))
plt.title("means", fontsize=16)
plt.xlabel("$-log10(<%s>)+%.2f$"%(normalisation_str,-minimum), fontsize=16)
plt.ylabel("pdf", fontsize=16)
plt.yscale('log')
plt.ylim(1e-4,1)
plt.xlim(data.min()-1,data.max()+1)
plt.legend(fontsize=16)
plt.show()
fig.savefig("var_distr_scaled.pdf")<jupyter_output><empty_output><jupyter_text>null<jupyter_code>df_null = pd.read_csv(("%s/nullTable.csv"%working_dir), header=None, index_col=[0])
df_null.head()<jupyter_output><empty_output><jupyter_text>meanVariances<jupyter_code>df_mv_null = pd.read_csv("meanVariances_null.csv", usecols=[1,2])
df_occ_null = pd.read_csv("O_null.dat", header=None)
df_mv_null.insert(2,'occurrence', np.array(df_occ_null.values,dtype=float))
#df_mv_null.to_csv("meanVariances_null.csv", index=False, header=True)
df_mv_null.head()
means_null = np.round(df_mv_null['mean'].values,1)
variances_null = np.round(df_mv_null['variance'].values,1)
occurrences_null = np.round(np.array(df_mv_null['occurrence'].values, dtype=float),1)
len(df_mv_null)
x = means
y = variances
x_lin = np.logspace(np.log10(x[x.nonzero()].min()),np.log10(x[x.nonzero()].max()), dtype=float,num=50)
high_means = df_mv[df_mv['occurrence']==1]['mean'].values
high_var = df_mv[df_mv['occurrence']==1]['variance'].values
fig=plt.figure(figsize=(12,7))

plt.scatter(x, y, label = 'genes', marker='o', alpha=0.8, linewidths=0.1)

#plt.scatter(high_means, high_var, label = '$O_i=1$', marker='o',c='cyan', alpha=0.4, lw=0.1)

log_bins_for_x = np.logspace(-5, np.log10(np.max(x)), num=50)
bin_means, bin_edges, binnumber = st.binned_statistic(x, y, statistic='mean', bins=log_bins_for_x)
bin_centres = (bin_edges[:-1]+bin_edges[1:])/2
plt.scatter((bin_edges[:-1]+bin_edges[1:])/2., bin_means, color='r', marker='x', lw=2, label='binned average')

plt.plot(x_lin[-40:],np.power(x_lin[-40:],2), 'g-', lw=3.5, label='$<%s>^2$ (Taylor)'%normalisation_str)
plt.plot(x_lin[:20],x_lin[:20], 'r-', lw=3.5, label='$<%s>$ (Poisson)'%normalisation_str)


#popt, pcov = curve_fit(lambda x,a,b : a*np.power(x,b), bin_centres, bin_means, bounds=([1,1],[35,5]))
#plt.plot(bin_centres, popt[0]*np.power(bin_centres, popt[1]), color='y', lw=3, label='fit')
#print(popt[0],popt[1])

#bin_sigmas,  bin_sigmas_edges, binsigmanumber = stats.binned_statistic(x, y, statistic=np.std, bins=log_bins_for_x)
#plt.plot((bin_edges[:-1] + bin_edges[1:])/2, bin_means+bin_sigmas*3, lw=3, color='yellow', label='binned average + $3\sigma$')

plt.scatter(means_null, variances_null, label='sampling')

plt.xlabel("$<%s>$"%normalisation_str, fontsize=16)
plt.ylabel("$\sigma^2_{%s}$"%normalisation_str, fontsize=16)
plt.legend(fontsize=20)
plt.xscale('log')
plt.yscale('log')
plt.xlim(x[x.nonzero()].min()/100,np.power(10,np.log10(x.max())+1))
plt.ylim((y[y.nonzero()].min()/100,np.power(10,np.log10(y.max())+1)))
plt.show()
fig.savefig("varmean_3sigma.png")
cv2 = np.array([variances[i]/(np.power(mean,2)) for i,mean in enumerate(means) if mean>0])
cv2_null = [variances_null[i]/(np.power(mean,2)) for i,mean in enumerate(means_null) if mean>0]
fig=plt.figure(figsize=(15,8))
plt.scatter(means[means.nonzero()], cv2, c='b', label ='genes')
plt.scatter(means_null[means_null.nonzero()], cv2_null, lw=3, c='orange', label='sampling')
plt.plot(x_lin[:30],1./x_lin[:30], 'r-', lw=1.5, label='Poisson')
plt.plot(x_lin[20:],[1e-1 for _ in x_lin[20:]], 'g-', lw=1.5, label='Taylor')
#plt.plot(x_lin,[nfiles-1 for _ in x_lin], color='cyan', lw=3, ls='--', label='bound')

bin_means, bin_edges, binnumber = st.binned_statistic(means[means.nonzero()], cv2, statistic='mean', bins=log_bins_for_x)
bin_centres = (bin_edges[:-1]+bin_edges[1:])/2
plt.scatter(np.array((bin_edges[:-1]+bin_edges[1:])/2.)[:-6], bin_means[:-6], color='r', marker='x', lw=2, label='binned average')

plt.xlabel("$<%s>$"%normalisation_str, fontsize=20)
plt.ylabel("$cv^2$", fontsize=20)
plt.xscale('log')
plt.yscale('log')
plt.xlim(means[means.nonzero()].min()/5,np.power(10,np.log10(means.max())+1))
plt.ylim((cv2[cv2.nonzero()].min()/10,np.power(10,np.log10(cv2.max())+1)))
plt.legend(fontsize=20)
plt.show()
fig.savefig("cvmean_loglog_sampling.png")
x = means
y = variances

# INIT FIGURE #################################################################

fig = plt.figure(figsize=(12, 6))
ax = fig.subplots()


# AX #########################################################################

xmin = np.log10(1e-3)
xmax = np.log10(x.max())
ymin = np.log10(1e-6)
ymax = np.log10(y.max())

nbins=80

xbins = np.logspace(xmin, xmax, nbins) # <- make a range from 10**xmin to 10**xmax
ybins = np.logspace(ymin, ymax, nbins) # <- make a range from 10**ymin to 10**ymax
counts, _, _, _ = ax.hist2d(x, y, bins=(xbins, ybins));

pcm = ax.pcolormesh(xbins, ybins, counts.T)
plt.colorbar(pcm)
#fig.colorbar(pcm, ax=ax)  # this works too


ax.set_xscale("log")               # <- Activate log scale on X axis
ax.set_yscale("log")               # <- Activate log scale on Y axis

ax.set_xlim(xmin=xbins[0])
ax.set_xlim(xmax=xbins[-1])
ax.set_ylim(ymin=ybins[0])
ax.set_ylim(ymax=ybins[-1])

ax.set_title("")
ax.set_xlabel("$<%s>$"%normalisation_str, fontsize=16)
ax.set_ylabel("$\sigma^2_{%s}$"%normalisation_str, fontsize=16)

# SHOW AND SAVE FILE ##########################################################

plt.tight_layout()
plt.show()
fig.savefig("varmean_density.png")
fig=plt.figure(figsize=(8,5))
plt.scatter(occurrences*nfiles, means, c='b', alpha=0.6, label='genes')

log_bins_for_x = np.logspace(-5, np.log10(np.max(x)), num=50)
bin_means, bin_edges, binnumber = st.binned_statistic(occurrences*nfiles, means, statistic='mean', bins=log_bins_for_x)
bin_centres = (bin_edges[:-1]+bin_edges[1:])/2
plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='r', lw=5, label='binned average')

plt.scatter(occurrences_null, means_null, c='orange', alpha=0.6, label='sampling')

plt.plot(np.linspace(1,nfiles), np.linspace(1,nfiles)/(nfiles), lw=3, label='$\\frac{O_i}{Nsamples}$', c='cyan', ls='--')
plt.ylabel("$<%s>$"%normalisation_str, fontsize=16)
#plt.xlabel("$\Sigma_j\Theta(FPKM-0.1)\Theta(1e5-FPKM)$", fontsize=16)
plt.xlabel("$\Sigma_j\Theta(%s)$"%normalisation_str, fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.ylim(5e-5,np.power(10,np.log10(means.max())+1))
plt.xlim(5e-1,nfiles+800)
plt.legend(fontsize=16, loc='upper left')
plt.show()
fig.savefig("meanDiff_binned_sampling.png")
x = occurrences
y = means

# INIT FIGURE #################################################################

fig = plt.figure(figsize=(12, 6))
ax = fig.subplots()


# AX #########################################################################

xmin = np.log10(9e-1)
xmax = np.log10(x.max())
ymin = np.log10(5e-4)
ymax = np.log10(y.max())

nbins=80

xbins = np.logspace(xmin, xmax, nbins) # <- make a range from 10**xmin to 10**xmax
ybins = np.logspace(ymin, ymax, nbins) # <- make a range from 10**ymin to 10**ymax
counts, _, _, _ = ax.hist2d(x, y, bins=(xbins, ybins));

pcm = ax.pcolormesh(xbins, ybins, counts.T)
plt.colorbar(pcm)
#fig.colorbar(pcm, ax=ax)  # this works too


ax.set_xscale("log")               # <- Activate log scale on X axis
ax.set_yscale("log")               # <- Activate log scale on Y axis

ax.set_xlim(xmin=xbins[0])
ax.set_xlim(xmax=xbins[-1])
ax.set_ylim(ymin=ybins[0])
ax.set_ylim(ymax=ybins[-1])

ax.set_title("")
ax.set_ylabel("$<%s>$"%normalisation_str, fontsize=16)
ax.set_xlabel("$\Sigma_j\Theta(%s)$"%normalisation_str, fontsize=16)

# SHOW AND SAVE FILE ##########################################################

plt.tight_layout()
plt.show()
fig.savefig("meanDiff_density.png")<jupyter_output><empty_output><jupyter_text>Length<jupyter_code>q_many = pd.read_csv("genes.txt", index_col=[0], header=[0])
q_many = q_many[q_many['type_of_gene']=='protein-coding']
lengths = q_many['exons']
df_mv.insert(3,'length',lengths)
df_mv.head()
from scipy.stats import binned_statistic
fig=plt.figure(figsize=(15,7))
lengths = df_mv['length'].values
means = df_mv['mean'].values
bin_means, bin_edges, _ = binned_statistic(lengths, means, statistic='mean', bins=np.logspace(1,7))
plt.scatter(lengths,means)
plt.scatter((bin_edges[1:]+bin_edges[:-1])/2., bin_means, marker='x')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('length (bp)', fontsize=16)
plt.ylabel('mean (counts)', fontsize=16)
plt.xlim((np.nanmin(lengths)/10,np.nanmax(lengths)*10))
plt.ylim((np.nanmin(means[means.nonzero()])/10,np.nanmax(means)*10))
plt.show()
fig.savefig("meanLength.pdf")
df_mv_counts=pd.read_csv("results/counts/meanVariances.csv", index_col=[0]).dropna(axis=0, how='any')
df_mv_counts.head()
df_mv_fpkm=pd.read_csv("results/proteincoding/meanVariances.csv", index_col=[0]).dropna(axis=0, how='any')
df_mv_fpkm.head()
df_mv_rpk=pd.read_csv("results/rpk/meanVariances.csv", index_col=[0]).dropna(axis=0, how='any')
df_mv_rpk.head()
diffnormmeans = []
for g in df_mv.index:
    try:
        diffnormmeans.append([g,df_mv.at[g,'length'],df_mv_counts.at[g,'mean'],df_mv_fpkm.at[g,'mean'], df_mv_rpk.at[g,'mean']])
    except:
        print("skipping %s"%g)
fig=plt.figure(figsize=(24,8))
ax=fig.subplots(1,3)
dataset = np.array(diffnormmeans).T
lengths = np.array(dataset[1],dtype=float)
means = np.array(dataset[2],dtype=float)
bin_means, bin_edges, _ = binned_statistic(lengths, means, statistic='mean', bins=np.logspace(1,7))
ax[0].scatter(lengths,means, label='genes')
ax[0].scatter((bin_edges[1:]+bin_edges[:-1])/2., bin_means, marker='x', label='binned average')
ax[0].set_yscale('log')
ax[0].set_xscale('log')
ax[0].set_xlabel('lenght (bp)', fontsize=16)
ax[0].set_ylabel('mean (counts)', fontsize=16)
ax[0].set_xlim((np.nanmin(lengths)/10,np.nanmax(lengths)*10))
ax[0].set_ylim((means[means.nonzero()].min()/10,means.max()*10))
means = np.array(dataset[3],dtype=float)
bin_means, bin_edges, _ = binned_statistic(lengths, means, statistic='mean', bins=np.logspace(1,7))
ax[1].scatter(lengths,means, label='genes')
ax[1].scatter((bin_edges[1:]+bin_edges[:-1])/2., bin_means, marker='x',label='binned average')
ax[1].plot(lenghts,np.power(lenghts,-0.5)*1e3, label='$L[kilobp]^{-0.5}$', lw=2, color='red')
ax[1].set_yscale('log')
ax[1].set_xscale('log')
ax[1].set_xlabel('lenght (bp)', fontsize=16)
ax[1].set_ylabel('mean (FPKM)', fontsize=16)
ax[1].set_xlim((np.nanmin(lengths)/10,np.nanmax(lengths)*10))
ax[1].set_ylim((means[means.nonzero()].min()/10,means.max()*10))
ax[1].legend(fontsize=16)
means = np.array(dataset[4],dtype=float)
bin_means, bin_edges, _ = binned_statistic(lengths, means, statistic='mean', bins=np.logspace(1,7))
ax[2].scatter(lengths,means, label='genes')
ax[2].plot(lengths,np.power(lengths,-1)*1e6, label='$L[kilobp]^{-1}$', lw=2, color='red')
ax[2].scatter((bin_edges[1:]+bin_edges[:-1])/2., bin_means, marker='x',label='binned average')
ax[2].set_yscale('log')
ax[2].set_xscale('log')
ax[2].set_xlabel('length (bp)', fontsize=16)
ax[2].set_ylabel('mean (RPK)', fontsize=16)
ax[2].set_xlim((np.nanmin(lengths)/10,np.nanmax(lengths)*10))
ax[2].set_ylim((means[means.nonzero()].min()/10,means.max()*10))
ax[2].legend(fontsize=16)
plt.show()
fig.savefig("meanLenght_fpkm_rpk.pdf")
scatterdense(np.array(dataset[2],dtype=float), np.array(dataset[4],dtype=float))
plt.plot(np.array(dataset[2],dtype=float),1e-3*np.power(np.array(dataset[2],dtype=float),1), c='r')
plt.xscale('log')
plt.yscale('log')
plt.xlim(np.nanmin(np.array(dataset[2],dtype=float)), np.nanmax(np.array(dataset[2],dtype=float)))
plt.ylim(np.nanmin(np.array(dataset[4],dtype=float)), np.nanmax(np.array(dataset[4],dtype=float)))
plt.show()<jupyter_output><empty_output><jupyter_text>P_i<jupyter_code>#query_g = df_mv[df_mv['mean']>1e-4].index.values
query_g = df_mv.index.values
len(query_g)
N = df.loc[query_g,:].sum(axis=0)
df_null.index=df.index
N_null = df_null.loc[query_g,:].sum(axis=0)
df_p = df.loc[query_g,:].div(N.values, axis=1)
df_p_null = df_null.loc[query_g,:].div(N_null.values, axis=1)
means = df_p.apply(np.average, axis=1).values
variances = df_p.apply(np.var, axis=1).values
o = df_p.apply(lambda x: float(len(np.nonzero(x)[0]))/len(x),axis=1).values
means_null = df_p_null.apply(np.average, axis=1).values
variances_null = df_p_null.apply(np.var, axis=1).values
o_null = df_p_null.apply(lambda x: float(len(np.nonzero(x)[0]))/len(x),axis=1).values
x = means
y = variances
log_bins_for_x = np.logspace(np.log10(x[x.nonzero()].min()),np.log10(x[x.nonzero()].max()), dtype=float,num=30)
x_null = means_null
y_null = variances_null
log_bins_for_x_null = np.logspace(np.log10(x_null[x_null.nonzero()].min()),np.log10(x_null[x_null.nonzero()].max()), dtype=float,num=30)
x_high = df_p.loc[df_mv[df_mv['occurrence']==1].index,:].mean(axis=1).values
y_high = df_p.loc[df_mv[df_mv['occurrence']==1].index,:].var(axis=1).values
fig=plt.figure(figsize=(12,7))

plt.scatter(x, y, label = 'genes', marker='o', alpha=0.8, linewidths=0.1)
plt.scatter(x_high, y_high, label = '$O_i=1$', marker='o', alpha=0.4, color='cyan', linewidths=0.1)
#plt.scatter(x_null, y_null, label = 'sampling', marker='o', alpha=0.8, linewidths=0.1)

bin_means, bin_edges, binnumber = stats.binned_statistic(x, y, statistic='mean', bins=log_bins_for_x)
bin_centres = (bin_edges[:-1]+bin_edges[1:])/2
plt.hlines(bin_means, bin_edges[:-1], bin_edges[1:], colors='r', lw=5, label='binned average')

plt.plot(log_bins_for_x[-40:],np.power(log_bins_for_x[-40:],2), 'g-', lw=3.5, label='$<p_i>^2$ (Taylor)')
plt.plot(log_bins_for_x[:],log_bins_for_x[:], 'r-', lw=3.5, label='$<p_i>$ (Poisson)')
plt.plot(log_bins_for_x[:],log_bins_for_x[:]*1e-8, ls='-', color='darkred', lw=3.5, label='$10^{-8}*<p_i>$ (Poisson)')



plt.xlabel("$<p_i>$", fontsize=16)
plt.ylabel("$\sigma^2_{p_i}$", fontsize=16)
plt.legend(fontsize=18)
plt.xscale('log')
plt.yscale('log')
plt.xlim(x[x.nonzero()].min()/100,np.power(10,np.log10(x.max())+1))
plt.ylim((y[y.nonzero()].min()/100,np.power(10,np.log10(y.max())+1)))
plt.show()
fig.savefig("pi_varmean_3sigma.png")
fig = plt.figure()
data = means
mu = np.median(data)
s = np.std(data)
log_bins = np.logspace(-12,-2)
plt.hist(data, density = False, histtype='step', bins=log_bins)
plt.title("means", fontsize=16)
plt.xlabel("$<pi>$", fontsize=16)
plt.ylabel("#", fontsize=16)
plt.yscale('log')
plt.xscale('log')
plt.xlim(data[data.nonzero()].min(),data.max())
plt.show()
fig.savefig("pi_mean_distr.pdf")<jupyter_output><empty_output><jupyter_text>over mean<jupyter_code>from matplotlib.patches import Rectangle
o_min = 3e1
o_max = nfiles
m_min = 5e4
m_max = 1e6
fig=plt.figure(figsize=(8,5))
plt.scatter(occurrences, means, c='b', alpha=0.6, label='data')
plt.plot(np.linspace(1,nfiles), np.linspace(1,nfiles)/(nfiles*10), lw=2, label='', c='r')

width = o_max - o_min
height = m_max-m_min
plt.gca().add_patch(Rectangle((o_min,m_min), width=width, height=height, fill=False))

plt.ylabel("$<%s>$"%normalisation_str, fontsize=16)
#plt.xlabel("$\Sigma_j\Theta(FPKM-0.1)\Theta(1e5-FPKM)$", fontsize=16)
plt.xlabel("$\Sigma_j\Theta(%s)$"%normalisation_str, fontsize=16)
plt.xscale('log')
plt.yscale('log')
plt.ylim(5e-5,np.power(10,np.log10(means.max())+1))
plt.xlim(5e-1,nfiles+800)
plt.show()
plt.savefig("highmean.png")
up = []
for g in df_mv.index:
    subdf = df_mv.loc[g,:]
    mean = subdf['mean']
    occ = subdf['occurrence']
    if mean>m_min and mean < m_max and occ*nfiles > o_min and occ*nfiles< o_max:
        up.append(g)
len(up)
for g in up:
    print(g)
dat = df.loc['ENSG00000042832'].values
np.var(dat)/(np.average(dat)**2)
df_mv[df_mv['mean']>0].sort_values(by='variance', axis=0, ascending=False)<jupyter_output><empty_output><jupyter_text>set by occurrence<jupyter_code>with open("o1.txt",'a') as f:
    for g in df_mv[df_mv['occurrence']>4990].index:
        f.write("%s\n"%g)<jupyter_output><empty_output><jupyter_text>data size Heaps check<jupyter_code>col = df.loc[:,df.keys()[1]].values
np.sum(col)
len(col[col.nonzero()])
x = []
y = []
for i in range(1, nfiles):
    col = df.loc[:,df.keys()[i]].values
    x.append(np.sum(col))
    y.append(len(col.nonzero()[0]))
plt.scatter(x,y)
i=794
x=[]
y=[]
col = df.loc[:,df.keys()[i]].values
x.append(np.sum(col))
y.append(len(col.nonzero()[0]))
x
y
col[8142:8150]<jupyter_output><empty_output><jupyter_text>Imagesave<jupyter_code>bits = np.array([df.loc[g,:].values for g in df.index])
from PIL import Image
img = Image.fromarray(bits, mode='1')
img.save("mat.bmp")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2318804586282545, 'fraction_numerical': 0.03468686617823585, 'mean_word_length': 7.801457194899818, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 87, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '34500', 'license': 'MIT', 'n_tokens_mistral': 11141, 'n_tokens_neox': 10396, 'n_words': 1272, 'path': 'Table_Analyzer.ipynb', 'repo_name': 'fvalle1/master_thesis', 'size': '24159'}