Scanpy code for BT145.

scanpy_clp_3_3_2020.py

# docker run -it -v "/home/cplaisier/Dropbox (ASU):/files" cplaisier/scrna_seq_velocity
# docker run -it -v "/home/swilferd:/files/scRNA_seq_Mehta" cplaisier/scrna_seq_velocity
# docker run -it -v "/home/swilferd:/files" cplaisier/scrna_seq_velocity
# pip3 intall mygene

#python3 entry:
import numpy as np
import pandas as pd
import scanpy as sc
import mygene as mg
from sklearn.metrics import adjusted_rand_score, normalized_mutual_info_score, adjusted_mutual_info_score
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')

fileVAR = {'C':{'n_genes':2500, 'percent_mito':0.25}, '3D': {}, '1W':{}, '2W':{}}
## C, 3D, 1W, 2W -- all complete 01/26/2020

## Load GB3_'fileVAR'
for i in fileVAR:
    results_file = './write/GB3_' + i + '.h5ad'
    adata = sc.read_10x_mtx(
            './GB3_' + i + '/filtered_feature_bc_matrix/',
            var_names = 'gene_symbols',
            cache=True)
    adata.var_names_make_unique()
    adata
    
    ## Preprocessing
    # Filter
    sc.pp.filter_cells(adata, min_genes=200)
    sc.pp.filter_genes(adata, min_cells=3)
    
    mito_genes = adata.var_names.str.startswith('MT-')
    adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
    adata.obs['n_counts'] = adata.X.sum(axis=1).A1
    adata.obs['n_counts'] = adata.X.sum(axis=1).A1
    adata.obs['n_genes'] = [i.count_nonzero() for i in adata.X]
    
    # Violin plots
    #sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'], jitter=0.4, multi_panel=True, save= '_GB3_' + i + '_scanpy.pdf')
    
    # Scatter plots
    #sc.pl.scatter(adata, x='n_counts', y='percent_mito', save= '_GB3_' + i + '_n_count_by_percMito.pdf')
    #sc.pl.scatter(adata, x='n_counts', y='n_genes', save= '_GB3_' + i + '_n_count_by_n_genes.pdf')
    
    # Filter
    adata = adata[adata.obs.n_genes > 2500, :]
    adata = adata[adata.obs.percent_mito < 0.25, :]
    print(adata.shape)

## Load BT145_'fileVAR'
BT145_adatas = {}
for i in fileVAR:
    results_file = './write/BT145_' + i + '.h5ad'
    adata = sc.read_10x_mtx(
            './BT145_' + i + '/filtered_feature_bc_matrix/',
            var_names = 'gene_symbols',
            cache=True)
    adata.var_names_make_unique()
    adata
    
    ## Preprocessing
    # Filter
    sc.pp.filter_cells(adata, min_genes=200)
    sc.pp.filter_genes(adata, min_cells=3)
    
    mito_genes = adata.var_names.str.startswith('MT-')
    adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
    adata.obs['n_counts'] = adata.X.sum(axis=1).A1
    adata.obs['n_counts'] = adata.X.sum(axis=1).A1
    adata.obs['n_genes'] = [i.count_nonzero() for i in adata.X]
    
    # Violin plots
    #sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'], jitter=0.4, multi_panel=True, save= '_' + i + '_scanpy.pdf')
    
    # Scatter plots
    #sc.pl.scatter(adata, x='n_counts', y='percent_mito', save= '_' + i + '_n_count_by_percMito.pdf')
    #sc.pl.scatter(adata, x='n_counts', y='n_genes', save= '_' + i + '_n_count_by_n_genes.pdf')
    
    # Filter
    adata = adata[adata.obs.n_genes > 2500, :]
    adata = adata[adata.obs.percent_mito < 0.25, :]
    print(adata.shape)
    adata.write(results_file)
    
    BT145_adatas[i] = adata


##########################################
### BT145_C 3867 cells and 18631 genes ###
##########################################

# Load meta data for BT145_C
meta_C = pd.read_csv('BT145_C/BT145_C_meta.data.csv', header=0, index_col=0)
meta_C.index = [i+'-1' for i in meta_C.index]

adata = BT145_adatas['C']
adata.obs['ccAF'] = meta_C.loc[adata.obs.index,'clusts_WT']

# Normalize
sc.pp.normalize_total(adata, target_sum=1e6)

# Log normalize
sc.pp.log1p(adata)

# Save raw data
adata.raw = adata

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=4000)
sc.pl.highly_variable_genes(adata, save='_C_hvg.pdf')
adata = adata[:, adata.var.highly_variable]

# Regress out unwated noise and scale data
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])
sc.pp.scale(adata, max_value=10)

# PCA
sc.tl.pca(adata, svd_solver='arpack')
#sc.pl.pca(adata, color='CST3')
#sc.pl.pca_variance_ratio(adata, log=True)
adata.write(results_file)
adata

# Compute the neighborhood graph
sc.pp.neighbors(adata) #, n_neighbors=10, n_pcs=40)

# UMAP
sc.tl.umap(adata)
#sc.pl.umap(adata, color=['leiden'], save='_BT145_C_clusters_resolution' + str(res1) + '.pdf')

# Clustering and finding marker genes
resolutions = [0.2, 0.3, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.525, 0.55, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0]
adjRandScores = []
amiScore = []
nmiScore = []
for res1 in resolutions:
    sc.tl.leiden(adata, resolution = res1)
    
    # Writing confusion matrix to CSV
    adjRandScores.append(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF']))
    amiScore.append(adjusted_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    nmiScore.append(normalized_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    #print('res = '+str(res1)+'; Adj Rand Score = '+str(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF'])))
    with open('BT145_C/BT145_C_confusion_matrix_'+str(res1)+'.csv', 'w') as outFile:
        pd.crosstab(adata.obs['leiden'],adata.obs['ccAF']).to_csv(outFile)
    
    
    # Plot UMAP and cluster
    sc.pl.umap(adata, color=['leiden'], save='_BT145_C_clusters_resolution' + str(res1) + '.pdf')
    
    # Identifying Marker Genes
    #sc.tl.dendrogram(adata, 'leiden')
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', corr_method='benjamini-hochberg', n_genes=10000)
    #sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, save='_BT145_C.pdf')
    #sc.tl.filter_rank_genes_groups(adata, min_fold_change=2)
    # visualize results
    #sc.pl.rank_genes_groups(adata, key='rank_genes_groups_filtered')
    # visualize results using dotplot
    #sc.tl.dendrogram(adata, 'leiden')
    #sc.pl.rank_genes_groups_dotplot(adata, key='rank_genes_groups_filtered', save='_BT145_C.pdf')
    
    # Writing Marker Genes to CSV
    with open('BT145_C/BT145_C_markergenes_dataframe_'+str(res1)+'.csv', 'w') as outFile:
        for x in range(max([int(i) for i in list(adata.obs['leiden'])])+1):
            #cluster = [x]*100
            #df = sc.get.rank_genes_groups_df(adata, str(x))
            df = sc.get.rank_genes_groups_df(adata, str(x), pval_cutoff=0.05, log2fc_min=1).sort_values(by='logfoldchanges', ascending=False)
            df['cluster'] = [x]*df.shape[0]
            if x == 0:
                df.to_csv(outFile, header = True)
            else:
                df.to_csv(outFile, header = False)


fig = plt.figure()
ax = plt.axes()
ax.plot(resolutions, adjRandScores)
ax.plot(resolutions, nmiScore, color='red')
ax.plot(resolutions, amiScore, color='green')
plt.xlabel('Leiden resolution')
plt.ylabel('Adjusted Rand Score')
plt.savefig('BT145_C/BT145_C_adjRandScores.pdf')

## Chose 0.475 based on the adjusted rand index
sc.tl.leiden(adata, resolution = 0.475)

from matplotlib.colors import LinearSegmentedColormap
sc.pl.umap(adata, color=['TIMP3', 'GFAP', 'FAM198B', 'PLAT', 'LONRF2', 'PTPRE', 'CRYAB', 'TIMP4'], color_map= LinearSegmentedColormap.from_list('abc',['gainsboro','mistyrose','firebrick'],N=100), save='_C_GFAP.pdf') #, vmin=2.5,vmax=6

                             


################
### BT145_3D ###
################

# Load meta data for BT145_3D
meta_3D = pd.read_csv('BT145_3D/BT145_3D_meta.data.csv', header=0, index_col=0)
meta_3D.index = [i+'-1' for i in meta_3D.index]

adata = BT145_adatas['3D']
adata.obs['ccAF'] = meta_3D.loc[adata.obs.index,'clusts_WT']

# Normalize
sc.pp.normalize_total(adata, target_sum=1e6)

# Log normalize
sc.pp.log1p(adata)

# Save raw data
adata.raw = adata

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=4000)
sc.pl.highly_variable_genes(adata, save='_3D_hvg.pdf')
adata = adata[:, adata.var.highly_variable]

# Regress out unwated noise and scale data
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])
sc.pp.scale(adata, max_value=10)

# PCA
sc.tl.pca(adata, svd_solver='arpack')
#sc.pl.pca(adata, color='CST3')
#sc.pl.pca_variance_ratio(adata, log=True)
adata.write(results_file)
adata

# Compute the neighborhood graph
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

# UMAP
sc.tl.umap(adata)

# Clustering and finding marker genes
resolutions = [0.2, 0.3, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.525, 0.55, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0]
adjRandScores = []
amiScore = []
nmiScore = []
for res1 in resolutions:
    sc.tl.leiden(adata, resolution = res1)
    
    # Writing confusion matrix to CSV
    adjRandScores.append(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF']))
    amiScore.append(adjusted_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    nmiScore.append(normalized_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    #print('res = '+str(res1)+'; Adj Rand Score = '+str(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF'])))
    with open('BT145_C/BT145_3D_confusion_matrix_'+str(res1)+'.csv', 'w') as outFile:
        pd.crosstab(adata.obs['leiden'],adata.obs['ccAF']).to_csv(outFile)
    
    
    # Plot UMAP and cluster
    sc.pl.umap(adata, color=['leiden'], save='_BT145_3D_clusters_resolution' + str(res1) + '.pdf')
    
    # Identifying Marker Genes
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', corr_method='benjamini-hochberg', n_genes=10000)
    #sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, save='_BT145_3D.pdf')
    #sc.tl.filter_rank_genes_groups(adata, min_fold_change=2)
    # visualize results
    #sc.pl.rank_genes_groups(adata, key='rank_genes_groups_filtered')
    # visualize results using dotplot
    #sc.tl.dendrogram(adata, 'leiden')
    #sc.pl.rank_genes_groups_dotplot(adata, key='rank_genes_groups_filtered', save='_BT145_3D.pdf')
    
    # Writing Marker Genes to CSV
    with open('BT145_3D/BT145_3D_markergenes_dataframe_'+str(res1)+'.csv', 'w') as outFile:
        for x in range(max([int(i) for i in list(adata.obs['leiden'])])+1):
            #cluster = [x]*100
            #df = sc.get.rank_genes_groups_df(adata, str(x))
            df = sc.get.rank_genes_groups_df(adata, str(x), pval_cutoff=0.05, log2fc_min=1).sort_values(by='logfoldchanges', ascending=False)
            df['cluster'] = [x]*df.shape[0]
            if x == 0:
                df.to_csv(outFile, header = True)
            else:
                df.to_csv(outFile, header = False)


fig = plt.figure()
ax = plt.axes()
ax.plot(resolutions, adjRandScores)
ax.plot(resolutions, nmiScore, color='red')
ax.plot(resolutions, amiScore, color='green')
plt.xlabel('Leiden resolution')
plt.ylabel('Adjusted Rand Score')
plt.savefig('BT145_3D/BT145_3D_adjRandScores.pdf')

## Chose 0.4 based on the adjusted rand index
#sc.tl.leiden(adata, resolution = 0.4)
#sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', corr_method='benjamini-hochberg', n_genes=10000)
#sc.tl.dendrogram(adata, 'leiden')
#sc.pl.rank_genes_groups_dotplot(adata, save='_3D_dotplot.pdf')


from matplotlib.colors import LinearSegmentedColormap
sc.pl.umap(adata, color=['TIMP3', 'GFAP', 'FAM198B', 'PLAT', 'LONRF2', 'PTPRE', 'CRYAB', 'TIMP4'], color_map= LinearSegmentedColormap.from_list('abc',['gainsboro','mistyrose','firebrick'],N=100), save='_3D_GFAP.pdf') #, vmin=2.5,vmax=6


################
### BT145_1W ###
################

# Load meta data for BT145_1W
meta_1W = pd.read_csv('BT145_1W/BT145_1W_meta.data.csv', header=0, index_col=0)
meta_1W.index = [i+'-1' for i in meta_1W.index]

adata = BT145_adatas['1W']
adata.obs['ccAF'] = meta_1W.loc[adata.obs.index,'clusts_WT']

# Normalize
sc.pp.normalize_total(adata, target_sum=1e6)

# Log normalize
sc.pp.log1p(adata)

# Save raw data
adata.raw = adata

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=4000)
sc.pl.highly_variable_genes(adata, save='_1W_hvg.pdf')
adata = adata[:, adata.var.highly_variable]

# Regress out unwated noise and scale data
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])
sc.pp.scale(adata, max_value=10)

# PCA
sc.tl.pca(adata, svd_solver='arpack')
#sc.pl.pca(adata, color='CST3')
#sc.pl.pca_variance_ratio(adata, log=True)
adata.write(results_file)
adata

# Compute the neighborhood graph
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

# UMAP
sc.tl.umap(adata)

# Clustering and finding marker genes
resolutions = [0.2, 0.3, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.525, 0.55, 0.575, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0]
adjRandScores = []
amiScore = []
nmiScore = []
for res1 in resolutions:
    sc.tl.leiden(adata, resolution = res1)
    
    # Writing confusion matrix to CSV
    adjRandScores.append(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF']))
    amiScore.append(adjusted_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    nmiScore.append(normalized_mutual_info_score(adata.obs['leiden'],adata.obs['ccAF']))
    #print('res = '+str(res1)+'; Adj Rand Score = '+str(adjusted_rand_score(adata.obs['leiden'],adata.obs['ccAF'])))
    with open('BT145_C/BT145_1W_confusion_matrix_'+str(res1)+'.csv', 'w') as outFile:
        pd.crosstab(adata.obs['leiden'],adata.obs['ccAF']).to_csv(outFile)
    
    
    # Plot UMAP and cluster
    sc.pl.umap(adata, color=['leiden'], save='_BT145_1W_clusters_resolution' + str(res1) + '.pdf')
    
    # Identifying Marker Genes
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon', corr_method='benjamini-hochberg', n_genes=10000)
    #sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False, save='_BT145_1W.pdf')
    #sc.tl.filter_rank_genes_groups(adata, min_fold_change=2)
    # visualize results
    #sc.pl.rank_genes_groups(adata, key='rank_genes_groups_filtered')
    # visualize results using dotplot
    #sc.tl.dendrogram(adata, 'leiden')
    #sc.pl.rank_genes_groups_dotplot(adata, key='rank_genes_groups_filtered', save='_BT145_1W.pdf')
    
    # Writing Marker Genes to CSV
    with open('BT145_1W/BT145_1W_markergenes_dataframe_'+str(res1)+'.csv', 'w') as outFile:
        for x in range(max([int(i) for i in list(adata.obs['leiden'])])+1):
            #cluster = [x]*100
            #df = sc.get.rank_genes_groups_df(adata, str(x))
            df = sc.get.rank_genes_groups_df(adata, str(x), pval_cutoff=0.05, log2fc_min=1).sort_values(by='logfoldchanges', ascending=False)
            df['cluster'] = [x]*df.shape[0]
            if x == 0:
                df.to_csv(outFile, header = True)
            else:
                df.to_csv(outFile, header = False)


fig = plt.figure()
ax = plt.axes()
ax.plot(resolutions, adjRandScores)
ax.plot(resolutions, nmiScore, color='red')
ax.plot(resolutions, amiScore, color='green')
plt.xlabel('Leiden resolution')
plt.ylabel('Adjusted Rand Score')
plt.savefig('BT145_1W/BT145_1W_adjRandScores.pdf')

## Chose 0.475 based on the adjusted rand index
sc.tl.leiden(adata, resolution = 0.475)


from matplotlib.colors import LinearSegmentedColormap
sc.pl.umap(adata, color=['TIMP3', 'GFAP', 'FAM198B', 'PLAT', 'LONRF2', 'PTPRE', 'CRYAB', 'TIMP4'], color_map= LinearSegmentedColormap.from_list('abc',['gainsboro','mistyrose','firebrick'],N=100), save='_1W_GFAP.pdf', vmin=2.5,vmax=6)