<jupyter_start><jupyter_text>Combine Results by Running the RF again on the pooled gene list<jupyter_code>classifier_adata = adata[adata.obs['Cell Type'].isin(celltype_list)]
classifier_adata = classifier_adata[:,gene_list_combined_df.index]



classifier_adata = classifier_adata[(classifier_adata.obs['corona_plus'] | classifier_adata.obs['corona_minus'] )]
classifier_adata.obs['DP_vs_DN'] = 'bla'
classifier_adata.obs.loc[classifier_adata.obs['corona_plus']==1,'DP_vs_DN'] = 'DP'
classifier_adata.obs.loc[classifier_adata.obs['corona_plus']==0,'DP_vs_DN'] = 'DN'

if 0 : 
    DP_genes_list = get_deg_edv(label = 'DP',deg_adata = classifier_adata,obs='DP_vs_DN' ,
                   n_obs = (classifier_adata.obs['DP_vs_DN'] == 'DP').sum() ,
                   return_genes = 10000 , log2fc_min=0).names.values

    classifier_adata = classifier_adata[:,np.setdiff1d(
                                            np.setdiff1d(DP_genes_list,[gene1,gene2]) , 
                                            mito_ribo_genes)] #deg_genes_list
### Run the Classifier for each celltype and save gene list

X = classifier_adata.X
y = classifier_adata.obs['corona_plus']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, stratify=y, random_state=0 )


clf = RandomForestClassifier(n_estimators=100, criterion='gini',random_state=0 ,
                             n_jobs = 50 , class_weight='balanced_subsample') #max_depth=4
clf.fit(X_train, y_train)
print(clf.score(X_test, y_test))
#clf.predict_proba(X_test )

importances = clf.feature_importances_
indices = np.argsort(importances)[::-1]
ordered_importances = importances[indices]
ordered_genes = classifier_adata.var.index[indices]
std = np.std([tree.feature_importances_ for tree in clf.estimators_],
             axis=0)
pd.DataFrame(ordered_genes).to_csv(result_prefix+'gene_list.csv')

feature_importance_df = pd.DataFrame( data=ordered_importances , columns = {'Feature Importance'},
                                index = ordered_genes )<jupyter_output>Trying to set attribute `.obs` of view, copying.<jupyter_text>Generate and Save Results and Plots<jupyter_code>go_top_n_genes =500
signature_genes = list(ordered_genes[:go_top_n_genes])
go_df = sc.queries.enrich(signature_genes, org="hsapiens", gprofiler_kwargs=dict(no_evidences=False))
go_df['neg_log10_p_value'] = -np.log10(go_df['p_value'])

go_df.to_csv(result_prefix+'go_enrichment.csv')
n_genes = 100
feature_importance_df = pd.DataFrame( data=ordered_importances[:n_genes] , columns = {'Feature Importance'},
                                    index = ordered_genes[0:n_genes] )

set_fig_params()
fig = plt.figure( figsize =(7.2*2 , 7.2*2) , dpi=200)

ax = sns.barplot(x='Feature Importance',y=ordered_genes[:n_genes],data=feature_importance_df ,  
                 edgecolor='black',linewidth = 0.25 , palette='RdPu_r')

ax.grid(False)

#plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
#plt.xticks([])
plt.title('ACE2+TMPRSS2+ | Top Enriched Genes')
plt.savefig(result_prefix+'ACE2+TMPRSS2+_top_enriched_genes.pdf')
#plt.show()
#plt.close(fig)


set_fig_params()


for source in np.unique(go_df['source']) : 
    fig = plt.figure( figsize =(7.2*8 , 7.2*2) , dpi=200)


    ax = sns.barplot(x='neg_log10_p_value',y='name', data= go_df[go_df['source']==source] ,  
                     edgecolor='black',linewidth = 0.25 , palette='RdPu_r')

    #plt.axvline(-np.log10(0.05) , c = 'k' , linestyle = 'dashed')
    #ax.grid(False)

    #plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    plt.xlabel('$-log_{10}(pval)$')
    plt.title('ACE2+TMPRSS2+ | Significant '+source+' Terms | '+' EDV')
    plt.savefig(result_prefix+'ACE2+TMPRSS2+_top_enriched_goterms_'+source+'.pdf')
    #plt.show()



sc.tl.score_genes(full_adata, signature_genes, score_name='ACE2+TMPRSS2+ signature score') 

if result_prefix ==  'internal_human_lung/': 
    scanpy.set_figure_params(scanpy=True, dpi=80, dpi_save=200, figsize =(7.2 , 7.2) )

    sc.pl.umap(full_adata , color =  ['ACE2+TMPRSS2+ signature score' , 'Cell Type' ] , size=(7.2,7.2) ,
              save='_signature_score_celltypes.pdf' , show =0)



if result_prefix ==  'lung_european/': 
    scanpy.set_figure_params(scanpy=True, dpi=80, dpi_save=200, figsize =(7.2 , 7.2) )

    sc.pl.umap(full_adata , color =  ['ACE2+TMPRSS2+ signature score' , 'ann_level_2' ] , size=(7.2,7.2) ,
              save='_signature_score_celltypes.pdf', show=0)


signature_classifier_adata = full_adata[classifier_adata.obs_names] 

Double_Positive = signature_classifier_adata[signature_classifier_adata.obs['corona_plus']].obs['ACE2+TMPRSS2+ signature score']
Double_Negative = signature_classifier_adata[signature_classifier_adata.obs['corona_minus']].obs['ACE2+TMPRSS2+ signature score']

ranksum_pval = scipy.stats.ranksums(Double_Positive,Double_Negative)
ks2samp_pval = scipy.stats.ks_2samp(Double_Positive,Double_Negative)






set_fig_params()

fig = plt.figure( figsize =(2*7.2, 2*7.2/1.6 ) , dpi=200)


if ((result_prefix ==  'internal_human_lung/') | (result_prefix ==  'lung_european/') |
   (result_prefix ==  'internal_fibroblast_only/') ): 
    sns.kdeplot((Double_Positive) , shade=1 , color=sns.color_palette("RdBu_r", 10)[-1] ,label='ACE2+TMPRSS2+ Cells' )
    sns.kdeplot((Double_Negative) , shade=1, color=sns.color_palette("RdBu_r", 10)[0], label='ACE2-TMPRSS2- Cells'  )
    plt.xlabel('Signature Score')


if ((result_prefix ==  'gut/') | (result_prefix ==  'nasal/')): 
    sns.kdeplot(np.log2(Double_Positive) , shade=1 , color=sns.color_palette("RdBu_r", 10)[-1] ,label='ACE2+TMPRSS2+ Cells' )
    sns.kdeplot(np.log2(Double_Negative) , shade=1, color=sns.color_palette("RdBu_r", 10)[0], label='ACE2-TMPRSS2- Cells'  )
    plt.xlabel('$log_{2}(Signature Score)$')



ax.grid(False)


plt.ylabel('Density')
plt.title('Wilcoxon Rank-sum test | Statistic = {:.3f} | p-val = {:.3f} \n \
K-S Two-sample test | Statistic = {:.3f} | p-val = {:.3f}'.format(ranksum_pval[0],ranksum_pval[1],
                                                                  ks2samp_pval[0],ks2samp_pval[1]))
plt.savefig(result_prefix+'ACE2+TMPRSS2+_kde_pvals.pdf')


plt.close('all')<jupyter_output>/ahg/regevdata/users/edv/software/anaconda3/envs/st/lib/python3.7/site-packages/pandas/core/series.py:853: RuntimeWarning: invalid value encountered in log2
  result = getattr(ufunc, method)(*inputs, **kwargs)
/ahg/regevdata/users/edv/software/anaconda3/envs/st/lib/python3.7/site-packages/statsmodels/nonparametric/kde.py:447: RuntimeWarning: invalid value encountered in greater
  X = X[np.logical_and(X > clip[0], X < clip[1])] # won't work for two columns.
/ahg/regevdata/users/edv/software/anaconda3/envs/st/lib/python3.7/site-packages/statsmodels/nonparametric/kde.py:447: RuntimeWarning: invalid value encountered in less
  X = X[np.logical_and(X > clip[0], X < clip[1])] # won't work for two columns.
{'fraction_non_alphanumeric': 0.15152832101984645, 'fraction_numerical': 0.029697233087063596, 'mean_word_length': 4.518784972022382, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '401125', 'license': 'MIT', 'n_tokens_mistral': 2797, 'n_tokens_neox': 2604, 'n_words': 365, 'path': 'Figure4_tissue_programs/analysis.ipynb', 'repo_name': 'JasmineTee/Covid_meta_analysis', 'size': '6903'}