lncrna_of_scanpy_test_with_imputation.py

# -*- coding: utf-8 -*-
"""lncRNA of ScanPy_Test with Imputation

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1EX7Bhp9RiZ9c1cBaYF8Uzd-WQi9ioPDq

# Installs
"""

!pip install scanpy

!pip3 install leidenalg

!pip install MulticoreTSNE

"""## Imputation Tool"""

!pip install deepimpute

"""# 3K PBMC Import"""

import scanpy as sc
import numpy as np
import pandas as pd

sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')

results_file = 'write/pbmc3k.h5ad'  # the file that will store the analysis results

!mkdir data
!wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz
!cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz
!mkdir write

pbmc = sc.read_10x_mtx('data/filtered_gene_bc_matrices/hg19/', var_names= 'gene_symbols',
                       cache = True)

"""# Filter lncRNA"""

lncF = 'https://raw.githubusercontent.com/josoga2/dataset-repos/main/lncTest/genes.tsv'
lnc = pd.read_csv(lncF, sep='\t', header=None)
lnc.columns = ['gene','name']
lnc.head(1)

pbmc.to_df().transpose().to_csv('PBMC.tsv', sep='\t')
PBMC = pd.read_csv('PBMC.tsv', sep='\t')

pbmcLNC = pd.merge(PBMC, lnc, left_on='Unnamed: 0', right_on='name', how='inner')

PBMC.head(5)

PBMC.shape

pbmcLNC = pbmcLNC.set_index('Unnamed: 0')
pbmcLNC = pbmcLNC.drop(['gene','name'], axis=1)
pbmcLNC.shape

"""# Dropout Imputation"""

from deepimpute.multinet import MultiNet
model = MultiNet()

#transpose Data
pbmcLNC = pbmcLNC.transpose()
print('Working on {} cells and {} genes'.format(*pbmcLNC.shape))

True in pbmcLNC.columns.duplicated()

def get_duplicate_cols(df: pd.DataFrame) -> pd.Series:
    return pd.Series(df.columns).value_counts()[lambda x: x>1]

get_duplicate_cols(pbmcLNC)

pbmcLNC = pbmcLNC.drop(['Y_RNA', 'Y_RNA', 'U1', '7SK'], axis=1)

pbmcLNC.shape

model.fit(pbmcLNC)

imputedData = model.predict(pbmcLNC)

imputedData.head(1)

import matplotlib.pyplot as plt
import numpy as np

limits = [0,100]

fig,ax = plt.subplots()

jitter = np.random.normal(0,1,pbmcLNC.size) # Add some jittering to better see the point density
ax.scatter(pbmcLNC.values.flatten()+jitter,imputedData.values.flatten(),s=2)
ax.plot(limits,limits,'r-.',linewidth=2)
ax.set_xlim(limits)
ax.set_ylim(limits)

plt.show()





adata = sc.AnnData(X= imputedData)

# not needed after data imputation
#adata = adata.T

adata.var_names_make_unique()

adata

"""# Preprocessing

Selection of highly variable genes per cell
"""

sc.pl.highest_expr_genes(adata, n_top=10)

sc.pp.filter_cells(adata, min_genes=10)

sc.pp.filter_genes(adata, min_cells=3)

adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

adata

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')

#Removing unwanted data of interest (mitochondiral genes and bad cells)

adata = adata[adata.obs.n_genes_by_counts < 20, :]
adata = adata[adata.obs.pct_counts_mt < 0.01, :]

"""# Select Highly Variable Data"""

#normalize total count data btw 1E0 - 1E4
sc.pp.normalize_total(adata, target_sum=10000)
#log of counts
sc.pp.log1p(adata)

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

import seaborn as sns

sns.set('paper', color_codes='CST3')

sc.pl.highly_variable_genes(adata)

adata

adata.raw = adata

"""# Visualization & Clustering"""

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

sc.tl.pca(adata, svd_solver='arpack')

sc.pl.pca(adata, color=['NEAT1','MALAT1','MIR142'])

sc.tl.umap(adata)

sc.pl.umap(adata, color=['NEAT1','MALAT1','MIR142'], use_raw=False)

sc.tl.leiden(adata)

sc.pl.umap(adata, color=['leiden'])

adata.write(results_file)

"""# Using Highly Variable Marker genes to profile individual cells"""

sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=30, sharey=False)

pbmc.write(results_file)

# USING tSNE

sc.tl.tsne(pbmc, perplexity=30, learning_rate=1000, random_state=0)

pbmc.write(results_file)

pbmc

sc.pl.tsne(pbmc, color='n_genes')

#I don't understand how this data was extracted 

marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
                'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']

adata = sc.read(results_file)

pd.DataFrame(adata.uns['rank_genes_groups']['names']).head()

result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
    {group + '_' + key[:1]: result[key][group]
    for group in groups for key in ['names', 'pvals']}).head()

sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=10)

sc.pl.rank_genes_groups_tracksplot(adata, groupby='leiden')

new_cluster_names = [
    'CD4 T', 'CD14 Monocytes',
    'B', 'CD8 T',
    'NK', 'FCGR3A Monocytes','UNK',
    'Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)

sc.pl.dotplot(adata, marker_genes, groupby='leiden')

sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90)