<jupyter_start><jupyter_text>Prepare master table<jupyter_code>import os
import pandas as pd
import numpy as np
import pickle
import lap_v2_py3 as lap_v2

reprocess_new_basis = False

#Folder with data:
folder = '../Data/'


if reprocess_new_basis:
    #Load in the conversion table
    conv = pd.read_csv(folder+'ann.csv',usecols=[1,2],index_col=0,squeeze=True)
    #Throw out ambiguous reads
    conv = conv[~conv.index.duplicated(keep=False)]
    #Load in the new basis data
    LU = pd.read_excel(folder+'InVivo.xlsx',sheetname='LU_POS vs FG',index_col=0,usecols=[0,5,6,7]).reindex(index=conv.keys()).mean(axis=1)
    FG = pd.read_excel(folder+'InVivo.xlsx',sheetname='LU_POS vs FG',index_col=0,usecols=[0,2,3,4]).reindex(index=conv.keys()).mean(axis=1)
    TH = pd.read_excel(folder+'InVivo.xlsx',sheetname='TH_POS vs FG',index_col=0,usecols=[0,5,6,7]).reindex(index=conv.keys()).mean(axis=1)
    BR = pd.read_excel(folder+'InVivo.xlsx',sheetname='BR_POS vs ECT',index_col=0,usecols=[0,2,3,4]).reindex(index=conv.keys()).mean(axis=1)
    ECT = pd.read_excel(folder+'InVivo.xlsx',sheetname='BR_POS vs ECT',index_col=0,usecols=[0,5,6,7]).reindex(index=conv.keys()).mean(axis=1)
    newdict = {'E.9 Lung Prim Nkx+':LU,'E.13 Thyroid':TH,'E.8.25 Foregut Endoderm':FG,'E.9 Forebrain':BR,'E.8.25 Ectoderm':ECT}

    #Reindex using Entrez ID's, add name, and average duplicates
    for name in newdict.keys():
        newdict[name].index=conv
        newdict[name].dropna(inplace=True)
        newdict[name].name = name
        temp = newdict[name].copy()
        newdict[name] = newdict[name][~newdict[name].index.duplicated()]
        for item in newdict[name].index:
            newdict[name].loc[item] = temp.loc[item].mean()
        del temp
    
    f = open(folder+'NewBasis.dat','wb')
    pickle.dump(newdict,f)
    f.close()
    
else:
    f = open(folder+'NewBasis.dat','rb')
    newdict = pickle.load(f)
    f.close()

#%% Load in the basis data
basis = pd.read_csv(folder+'Mouse_Basis_Data.txt',delimiter='\t',index_col=0,usecols=[0]+list(range(3,64)))

thresh = 1
# Append new basis data
basis_new = basis.copy()
newdict_log = {}
for name in newdict.keys():
    newdict_log[name] = np.log2(newdict[name]+1)
    basis_new = basis_new.join(newdict_log[name][newdict_log[name]>thresh],how='inner')
basis_new.dropna(inplace=True)

#Load Keri's data
keri_index = pd.read_table('../Alex_Analysis/Final_Data/entrez_id.txt',index_col=None,header=None).squeeze().values
keri_label = pd.read_table('../Alex_Analysis/Final_Data/keri_cell_lbl.txt',index_col=None,header=None).squeeze()
keri = pd.read_table('../Alex_Analysis/Final_Data/keri_ranknorm_data_corr.txt',index_col=None,header=None)
keri.index=keri_index
keri = keri.rename(columns=keri_label)

# Load original project data 
data = pd.read_excel(folder+'InVitroNew.xlsx',index_col=0,usecols=[0]+list(range(38,50)))

#Load new data, reindex using Entrez ID's, and average duplicates
#Load in the conversion table
conv = pd.read_table(folder+'AnnotationTable_mogene10sttranscriptcluster.txt',index_col=1,usecols=[2,3],squeeze=True)
#Throw out ambiguous reads
conv = conv[~conv.index.duplicated(keep=False)]
data_new = pd.read_excel(folder+'rlog_7-2018.xlsx',index_col=0).reindex(index=conv.keys())
data_new.index=conv
data_new.dropna(inplace=True)
data_new = data_new[(data_new.T != 0).any()]
temp = data_new.copy()
data_new = data_new[~data_new.index.duplicated()]
#Right now, not averaging duplicates, because it is slow and doesn't matter
#for label in data_new.keys():
#    for item in data_new.index:
#        data_new.loc[item,label] = temp[label].loc[item].mean()
del temp

#Load old test data, reindex using Entrez ID's, and average duplicates
conv_probe = pd.read_table(folder+'AnnotationTable_mogene10sttranscriptcluster.txt',index_col=0,usecols=[0,2],squeeze=True)
data_old = pd.read_table(folder+'InVitroOld.txt',index_col=0,skipfooter=1,engine='python').reindex(index=conv_probe.keys())
data_old.index = conv_probe
data_old_dedup = data_old[~data_old.index.duplicated()].copy()
#for item in data_old_dedup.index:
#    data_old_dedup.loc[item] = data_old.loc[item].mean()

# Combine into master table
master_table = basis_new.join(data,how='inner')
master_table = master_table.join(data_new,how='inner')
master_table = master_table.join(data_old_dedup,how='inner')
master_table = master_table.join(keri,how='inner')

# Rank-norm data
N = np.shape(master_table)[0]
for item in master_table:
    master_table[item] = lap_v2.rank_norm(np.asarray(master_table[item]), dist='normal',norm=N)

data_keys = data.keys()
data_keys = data_keys.append(data_old.keys())
data_keys = data_keys.append(data_new.keys())
data_keys = data_keys.append(keri.keys())

f = open(folder+'master_table.dat','wb')
pickle.dump([master_table,data_keys],f)
f.close()<jupyter_output><empty_output><jupyter_text>Make LAP function for pandas<jupyter_code>import os
os.chdir('/Users/robertmarsland/Dropbox/BU/Laertis-Analysis/Alex_Analysis/Code/')
import pandas as pd
import numpy as np
import pickle
import lap_v2
import matplotlib.pyplot as plt
import matplotlib
matplotlib.style.use('ggplot')
%matplotlib inline

#Folder with data:
folder = '../../Data/'

f = open(folder+'master_table.dat','rb')
master_table, data_keys = pickle.load(f)
f.close()

master_table.dropna(inplace=True)

#Remove Nkx2-1
master_table.drop(21869,axis=0,inplace=True)

def pd_lap(names, exceptions = None, data = master_table):
    basis = data.drop(data_keys,1)
    
    if exceptions != None:
        reduced_basis = basis.drop(exceptions,1)
    else:
        reduced_basis = basis
        
    [a, A, eta] = lap_v2.lap(reduced_basis.values, data[names].values, full_output=True)

    projections = pd.DataFrame(data = a, index = reduced_basis.keys(), columns = names)
    
    return projections<jupyter_output><empty_output><jupyter_text>Average groups of cell types that are too similar<jupyter_code>def average_groups(data,name_dict):
    names_all = [item for sublist in name_dict.values() for item in sublist]
    reduced_data = data.drop(names_all,1)
    
    for name in name_dict:
        reduced_data = reduced_data.join(pd.DataFrame(data = data[name_dict[name]].mean(1), columns = [name]))

    return reduced_data
namedict = {'epithelial':['medullary thymic epithelial', 'terminal bronchiolar epithelial cells', 
                      'endothelial cell - gonad', 'endothelial cell - blood vessel',
                      'endothelial cell - lymphatic vessel'],
            'monocytes':['monocytes - classical', 'Granulocyte-Monocyte Progenitor (GMP)'],
            'adipose':['adipose - brown','adipose - white'],
            'killer cells':['natural killer cells','T Cell'],
            'intestine':['intestine - small','intestine - large'],
            'neuron':['neuron - glutamatergic','nsc'],
           'fibroblast':['fibroblast - thymus','fibroblast - skin','fibroblast - cardiac']}
table_fixed = average_groups(master_table,namedict)<jupyter_output><empty_output><jupyter_text>Apply to data<jupyter_code>thyroid_protocol = pd_lap(['BMP4+FGF2 mCherry+ 1','BMP4+FGF2 mCherry+ 2','BMP4+FGF2 mCherry+ 3','BMP4+FGF2 mCherry- 1','BMP4+FGF2 mCherry- 2','BMP4+FGF2 mCherry- 3'],
                          data = master_table)
lung_protocol = pd_lap(['BMP4+Wnt3A mCherry+ 1','BMP4+Wnt3A mCherry+ 2','BMP4+Wnt3A mCherry+ 3','BMP4+Wnt3A mCherry- 1','BMP4+Wnt3A mCherry- 2','BMP4+Wnt3A mCherry- 3'],
                       data = master_table)
#old_protocol = pd_lap(['GSM861609','GSM861610','GSM861611','GSM861612','GSM861613','GSM861615'],
#                      data = master_table)
new_protocols = pd_lap(['3D_EpCAM_2', '3D_EpCAM_3', '3D_EpCAM_4', '3D_mCherry-EpCAM_2','3D_mCherry-EpCAM_3', '3D_mCherry-EpCAM_4', 'Gelatin_mCherry-EpCAM_2','Gelatin_mCherry-EpCAM_3', 'Gelatin_mCherry-EpCAM_4'],
                      data = master_table)
keri_protocol = pd_lap(['D1-1','D1-2','D7-1','D7-2','D7-3','D14-1','D14-2'],data = master_table)
thyroid_protocol = pd_lap(['BMP4+FGF2 mCherry+ 1','BMP4+FGF2 mCherry+ 2','BMP4+FGF2 mCherry+ 3','BMP4+FGF2 mCherry- 1','BMP4+FGF2 mCherry- 2','BMP4+FGF2 mCherry- 3'],
              exceptions = ['Common Lymphoid Progenitor (CLP)','thyroid','lung','HSC','sperm','oocyte','mef'],
                         data = table_fixed)
lung_protocol = pd_lap(['BMP4+Wnt3A mCherry+ 1','BMP4+Wnt3A mCherry+ 2','BMP4+Wnt3A mCherry+ 3','BMP4+Wnt3A mCherry- 1','BMP4+Wnt3A mCherry- 2','BMP4+Wnt3A mCherry- 3'],
              exceptions = ['Common Lymphoid Progenitor (CLP)','thyroid','lung','HSC','sperm','oocyte','mef'],
                      data = table_fixed)
#old_protocol = pd_lap(['GSM861609','GSM861610','GSM861611','GSM861612','GSM861613','GSM861615'],
#                     exceptions = ['Common Lymphoid Progenitor (CLP)','thyroid','lung','HSC','sperm','oocyte','mef'],
#                     data = table_fixed)
new_protocols = pd_lap(['3D_EpCAM_2', '3D_EpCAM_3', '3D_EpCAM_4', '3D_mCherry-EpCAM_2','3D_mCherry-EpCAM_3', '3D_mCherry-EpCAM_4', 'Gelatin_mCherry-EpCAM_2','Gelatin_mCherry-EpCAM_3', 'Gelatin_mCherry-EpCAM_4'],
                        exceptions = ['Common Lymphoid Progenitor (CLP)','thyroid','lung','HSC','sperm','oocyte','mef'],
                        data = table_fixed)
keri_protocol = pd_lap(['D1-1','D1-2','D7-1','D7-2','D7-3','D14-1','D14-2'],
                        exceptions = ['Common Lymphoid Progenitor (CLP)','thyroid','lung','HSC','sperm','oocyte','mef'],
                        data = table_fixed)
fig,axs=plt.subplots(1,4,figsize=(24,32))
new_protocols.plot.barh(ax=axs[0])
thyroid_protocol.plot.barh(ax=axs[1])
lung_protocol.plot.barh(ax=axs[2])
keri_protocol.plot.barh(ax=axs[3])
axs[0].set_xlim((-0.2,0.7))
axs[1].set_xlim((-0.2,0.7))
axs[2].set_xlim((-0.2,0.7))
axs[3].set_xlim((-0.2,0.7))
axs[1].set_yticklabels(())
axs[2].set_yticklabels(())
axs[3].set_yticklabels(())
axs[0].set_title('New Data')
axs[1].set_title('Old Thyroid Protocol')
axs[2].set_title('Old Lung Protocol')
axs[3].set_title('Keri Data')
plt.show()<jupyter_output><empty_output><jupyter_text>Project noise  Multiplicative noise Make noise amplitude proportional to square root of mean abundance. Assume rank-norm correctly extracts original abundance, with mean subtracted off. Add the mean back in (wild guess at 5 times standard deviation of expression levels) before taking square root. These plots show the result when constant of proportionality is set to 0.1, so that noise amplitude at high expression level (4 standard deviations out) is equal to 1/5 of the standard deviation.<jupyter_code>N = np.shape(master_table)[0]
noise = np.random.randn(N,3)#*np.sqrt(master_table[['BMP4+FGF2 mCherry+ 1','BMP4+FGF2 mCherry+ 2','BMP4+FGF2 mCherry+ 3']].values+5.)
[a, A, eta] = lap_v2.lap(basis_new.values, noise, full_output=True)
projections = pd.DataFrame(data = a, index = basis_new.keys(), columns = ['Noise1','Noise2','Noise3'])
projections.plot.barh(figsize=(6, 24))
plt.show()<jupyter_output><empty_output><jupyter_text>Project everything onto everything<jupyter_code>item = basis_new.keys()[0]
reduced_basis = basis_new.drop(item,1)
[a, A, eta] = lap_v2.lap(reduced_basis.values, master_table[item].values, full_output=True)
projections = pd.DataFrame(data = a, index = reduced_basis.keys(), columns = [item])

for item in basis_new.keys()[1:]:
    reduced_basis = basis_new.drop(item,1)
        
    [a, A, eta] = lap_v2.lap(reduced_basis.values, master_table[item].values, full_output=True)
    projections_new = pd.DataFrame(data = a, index = reduced_basis.keys(), columns = [item])
    projections = projections.join(projections_new)

projections = projections.drop('oocyte',1)
projections[projections.isnull()] = 0
plt.pcolor(projections)
plt.colorbar()
plt.show()
plt.hist(projections)
pairs=[[]]
select = projections.values > 0.5
for j in range(len(projections.keys())):
    for k in range(len(projections.keys())):
        if select[j,k]:
            print projections.keys()[k] + ' on ' + projections.keys()[j]
'adipose':['adipose - brown','adipose - white'],
'killer cells':['natural killer cells','T Cell'],
'intestine':['intestine - small','intestine - large'],
'neuron':['neuron - glutamatergic','nsc']
projections[['monocytes - classical','fibroblast - skin','neuron - glutamatergic']].plot.barh(figsize=(6, 24))
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup