<jupyter_start><jupyter_text>import scVI models<jupyter_code>from scvi.inference import UnsupervisedTrainer, SemiSupervisedTrainer
from scvi.models.scanvi import SCANVI
from scvi.models.vae import VAE
import numpy as np
gene_dataset = ZeiselMoleArchData()
gene_dataset_fine = ZeiselMoleArchData(coarse=False)
print("The number of classes in the coarse data is %i\n" % len(np.unique(gene_dataset.labels)))
print("The number of classes in the fine data is %i\n" % len(np.unique(gene_dataset_fine.labels)))
if np.sum([gene_dataset_fine.labels_groups[x] for x in gene_dataset_fine.labels.ravel()] == gene_dataset.labels.ravel()) == len(gene_dataset_fine):
    print("the mapping of fine to coarse labels is Correct\n") 
else: print("check that labels_groups is right")
     

import pandas as pd
genes1 = pd.read_table('../Zeisel/Zeisel.1.hvg_info.csv', delimiter=',')
geneid1 = np.asarray([x.replace('gene_', '') for x in genes1[genes1.keys()[0]]]).astype('int')
genenames1 = genes1['genename']
assert np.sum(np.asarray(genenames1) == gene_dataset.gene_names) == len(gene_dataset.gene_names)
genes = gene_dataset.gene_names[geneid1[:1000]]
from scvi.harmonization.utils_chenling import subsetByGenenames
gene_dataset = subsetByGenenames(gene_dataset,genes)
gene_dataset_fine = subsetByGenenames(gene_dataset_fine,genes)
vae = VAE(gene_dataset.nb_genes, n_batch=gene_dataset.n_batches, n_labels=gene_dataset.n_labels,
          n_hidden=128, n_latent=10, n_layers=2, dispersion='gene')
trainer = UnsupervisedTrainer(vae, gene_dataset, train_size=1.0)
# trainer.model.load_state_dict(torch.load('../Zeisel/vae.pkl'))
# trainer.model.eval()
# trainer.train(n_epochs=100)
# torch.save(trainer.model.state_dict(),'../Zeisel/vae.pkl')
full = trainer.create_posterior(trainer.model, gene_dataset, indices=np.arange(len(gene_dataset)))<jupyter_output><empty_output><jupyter_text>SCANVI coarse<jupyter_code>scanvi = SCANVI(gene_dataset.nb_genes, gene_dataset.n_batches, gene_dataset.n_labels, n_layers=2)
scanvi.load_state_dict(full.model.state_dict(), strict=False)
trainer_scanvi = SemiSupervisedTrainer(scanvi, gene_dataset, classification_ratio=50,
                                       n_epochs_classifier=1, lr_classification=5 * 1e-3)
# trainer_scanvi.train(n_epochs=5)
# torch.save(trainer_scanvi.model.state_dict(),'../Zeisel/scanvi.coarse.pkl')
trainer_scanvi.model.load_state_dict(torch.load('../Zeisel/scanvi.coarse.pkl'))
trainer_scanvi.model.eval()
full_coarse = trainer_scanvi.create_posterior(trainer_scanvi.model, gene_dataset, indices=np.arange(len(gene_dataset)))<jupyter_output><empty_output><jupyter_text>SCANVI Fine<jupyter_code>K=5
# ,'dropout_rate':0.2
scanvi = SCANVI(gene_dataset_fine.nb_genes, gene_dataset_fine.n_batches, gene_dataset_fine.n_labels, n_layers=2,
                classifier_parameters = {'dropout_rate':0.2, 'n_hidden':256, 'n_layers':2})
scanvi.load_state_dict(full.model.state_dict(), strict=False)
trainer_scanvi = SemiSupervisedTrainer(scanvi, gene_dataset_fine, classification_ratio=100,
                                       n_epochs_classifier=1, lr_classification=5 * 1e-3,n_labelled_samples_per_class=K)
trainer_scanvi.train(n_epochs=5)
torch.save(trainer_scanvi.model.state_dict(),'../Zeisel/scanvi.%i.pkl'%K)
trainer_scanvi.model.load_state_dict(torch.load('../Zeisel/scanvi.%i.pkl'%K))
trainer_scanvi.model.eval()
full_fine = trainer_scanvi.create_posterior(trainer_scanvi.model, gene_dataset_fine, indices=np.arange(len(gene_dataset_fine)))
unlabelled_idx= trainer_scanvi.unlabelled_set.indices
labelled_idx= trainer_scanvi.labelled_set.indices
_,SCANVI_fine_pred = full_fine.sequential().compute_predictions()
len(np.unique(gene_dataset_fine.labels[unlabelled_idx]))
len(np.unique(SCANVI_fine_pred[unlabelled_idx]))
scanvi_fine_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],SCANVI_fine_pred[unlabelled_idx])
scanvi_fine_acc[1]<jupyter_output><empty_output><jupyter_text>SCANVI Hierarchical<jupyter_code>scanvi = SCANVI(gene_dataset_fine.nb_genes, gene_dataset_fine.n_batches, gene_dataset_fine.n_labels, n_layers=2,
                classifier_parameters = {'dropout_rate':0.2, 'n_hidden':256, 'n_layers':2},
               labels_groups=gene_dataset_fine.labels_groups, use_labels_groups=True)

scanvi.load_state_dict(full.model.state_dict(), strict=False)


trainer_scanvi = SemiSupervisedTrainer(scanvi, gene_dataset_fine, classification_ratio=100,
                                       n_epochs_classifier=1, lr_classification=5 * 1e-3,n_labelled_samples_per_class=K)

trainer_scanvi.train(n_epochs=5)
torch.save(trainer_scanvi.model.state_dict(),'../Zeisel/scanvi.hier.%i.pkl'%K)
trainer_scanvi.model.load_state_dict(torch.load('../Zeisel/scanvi.hier.%i.pkl'%K))
trainer_scanvi.model.eval()
full_hier = trainer_scanvi.create_posterior(trainer_scanvi.model, gene_dataset_fine, indices=np.arange(len(gene_dataset_fine)))
_,SCANVI_hier_pred = full_hier.sequential().compute_predictions()
scanvi_hier_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],SCANVI_hier_pred[unlabelled_idx])

len(np.unique(SCANVI_hier_pred[unlabelled_idx]))
scanvi_hier_acc[1]<jupyter_output><empty_output><jupyter_text>scVI and Seurat latent space<jupyter_code>vae_posterior = trainer.create_posterior(trainer.model, gene_dataset, indices=np.arange(len(gene_dataset)))
seurat_latent = np.genfromtxt('../Zeisel/Zeisel.1.CCA.txt')
latent, batch_indices, labels = vae_posterior.sequential().get_latent()
batch_indices = batch_indices.ravel()<jupyter_output><empty_output><jupyter_text>Coarse to Coarse<jupyter_code>import copy
gene_dataset_coarse = copy.deepcopy(gene_dataset_fine)
gene_dataset_coarse.labels = gene_dataset_fine.labels_groups[gene_dataset_fine.labels.ravel()].reshape(len(gene_dataset_fine),1)
gene_dataset_coarse.cell_types = gene_dataset_fine.groups
gene_dataset_coarse.n_labels = gene_dataset_fine.groups.shape[0]
scanvi = SCANVI(gene_dataset.nb_genes, gene_dataset.n_batches, gene_dataset.n_labels, n_layers=2)
scanvi.load_state_dict(full.model.state_dict(), strict=False)
trainer_scanvi = SemiSupervisedTrainer(scanvi, gene_dataset, classification_ratio=50,
                                       n_epochs_classifier=1, lr_classification=5 * 1e-3)
labelled_idx = trainer_scanvi.labelled_set.indices
unlabelled_idx= trainer_scanvi.unlabelled_set.indices
from sklearn.neighbors import KNeighborsClassifier
def PartialPrediction(latent,labelled_idx,unlabelled_idx,labels):
        latent_labelled = latent[labelled_idx, :]
        latent_unlabelled = latent[unlabelled_idx, :]
        labels_labelled = labels[labelled_idx]
        labels_unlabelled = labels[unlabelled_idx]
        neigh = KNeighborsClassifier(n_neighbors=10)
        neigh = neigh.fit(latent_labelled, labels_labelled)
        labels_pred = neigh.predict(latent_unlabelled)
        return labels_pred

coarse_latent, _, _ = full_coarse.sequential().get_latent()
coarse_pred=(PartialPrediction(coarse_latent, labelled_idx,unlabelled_idx,gene_dataset.labels.ravel()))

# seurat coarse
seurat_coarse_pred=(PartialPrediction(seurat_latent, labelled_idx,unlabelled_idx,gene_dataset.labels.ravel()))
# scVI coarse
scVI_coarse_pred=(PartialPrediction(latent, labelled_idx,unlabelled_idx,gene_dataset.labels.ravel()))
coarse_pred=(PartialPrediction(coarse_latent, labelled_idx,unlabelled_idx,gene_dataset.labels.ravel()))
y,SCANVI_coarse_pred = full_coarse.sequential().compute_predictions()
seurat_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx], seurat_coarse_pred)
scVI_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],scVI_coarse_pred)
coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],coarse_pred)
scanvi_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],SCANVI_coarse_pred[unlabelled_idx])
seurat_coarse_acc[1]
scVI_coarse_acc[1]
scanvi_coarse_acc[1]
coarse_acc[1]<jupyter_output><empty_output><jupyter_text>Fine to Fine<jupyter_code>scanvi = SCANVI(gene_dataset_fine.nb_genes, gene_dataset_fine.n_batches, gene_dataset_fine.n_labels, n_layers=2)
scanvi.load_state_dict(full.model.state_dict(), strict=False)
trainer_scanvi = SemiSupervisedTrainer(scanvi, gene_dataset_fine, classification_ratio=50,
                                       n_epochs_classifier=1, lr_classification=5 * 1e-3,n_labelled_samples_per_class=K)
labelled_idx = trainer_scanvi.labelled_set.indices
unlabelled_idx= trainer_scanvi.unlabelled_set.indices
scanvi_latent,_,_ = full_fine.sequential().get_latent()
scanvi_hier_latent,_,_ = full_hier.sequential().get_latent()
len(labelled_idx)
seurat_fine_pred=(PartialPrediction(seurat_latent, labelled_idx,unlabelled_idx,gene_dataset_fine.labels.ravel()))
scVI_fine_pred=(PartialPrediction(latent, labelled_idx,unlabelled_idx,gene_dataset_fine.labels.ravel()))
fine_pred=(PartialPrediction(scanvi_latent, labelled_idx,unlabelled_idx,gene_dataset_fine.labels.ravel()))
hier_pred=(PartialPrediction(scanvi_hier_latent, labelled_idx,unlabelled_idx,gene_dataset_fine.labels.ravel()))
scVI_fine_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],scVI_fine_pred)
seurat_fine_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],seurat_fine_pred)

# _,SCANVI_fine_pred = full_fine.sequential().compute_predictions()
# _,SCANVI_hier_pred = full_hier.sequential().compute_predictions()
scanvi_fine_acc[1]
scanvi_hier_acc[1]
scVI_fine_acc[1]
seurat_fine_acc[1]
fine_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],fine_pred)
fine_acc[1]
hier_acc = compute_accuracy_tuple(gene_dataset_fine.labels.ravel()[unlabelled_idx],hier_pred)
hier_acc[1]
print("%.3f & %.3f &%.3f &%.3f & %.3f & %.3f \\\\"%(seurat_fine_acc[0], scVI_fine_acc[0],scanvi_fine_acc[0],fine_acc[0],scanvi_hier_acc[0],hier_acc[0]))
print("%.3f & %.3f &%.3f &%.3f & %.3f & %.3f \\\\"%(seurat_fine_acc[1], scVI_fine_acc[1],scanvi_fine_acc[1],fine_acc[1],scanvi_hier_acc[1],hier_acc[1]))<jupyter_output>0.183 & 0.244 &0.367 &0.268 & 0.365 & 0.276 \\<jupyter_text>Fine to Coarse<jupyter_code>len(SCANVI_hier_pred)
len(fine_pred)
seurat_fine_pred_coarse = gene_dataset_fine.labels_groups[seurat_fine_pred]
scVI_fine_pred_coarse = gene_dataset_fine.labels_groups[scVI_fine_pred]
SCANVI_fine_pred_coarse = gene_dataset_fine.labels_groups[SCANVI_fine_pred]
SCANVI_hier_pred_coarse = gene_dataset_fine.labels_groups[SCANVI_hier_pred]
fine_pred_coarse = gene_dataset_fine.labels_groups[fine_pred]
hier_pred_coarse = gene_dataset_fine.labels_groups[hier_pred]
seurat_fine_coarse_acc = compute_accuracy_tuple(seurat_fine_pred_coarse,gene_dataset.labels.ravel()[unlabelled_idx])
scVI_fine_coarse_acc = compute_accuracy_tuple(scVI_fine_pred_coarse,gene_dataset.labels.ravel()[unlabelled_idx])
scanvi_fine_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],SCANVI_fine_pred_coarse[unlabelled_idx])
scanvi_hier_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],SCANVI_hier_pred_coarse[unlabelled_idx])
fine_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],fine_pred_coarse)
hier_coarse_acc = compute_accuracy_tuple(gene_dataset.labels.ravel()[unlabelled_idx],hier_pred_coarse)

seurat_fine_coarse_acc
scVI_fine_coarse_acc
scanvi_fine_coarse_acc
scanvi_hier_coarse_acc
fine_coarse_acc
hier_coarse_acc<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11414847161572052, 'fraction_numerical': 0.012663755458515284, 'mean_word_length': 8.728971962616823, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '273148', 'license': 'MIT', 'n_tokens_mistral': 4557, 'n_tokens_neox': 4200, 'n_words': 413, 'path': 'notebooks/Hierarchical-final.ipynb', 'repo_name': 'chenlingantelope/HarmonizationSCANVI', 'size': '11450'}