scanpy 去批次pipeline

import scanpy as sc
import anndata as an
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import scrublet as scr
import scipy.io
import numpy as np
import os
import pandas as pd
from math import sqrt
import bbknn  
import random
from matplotlib.pyplot import rc_context
import datetime
import DIY
from DIY import get_tsne
sc.logging.print_versions()
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3

random.seed(234)
figure_out = './figures'
results_file = 'write/sce.h5ad'

table_out = './write'
datadir='./CellrangerOut/'
# df=pd.read_csv('sample_info.txt',header=0,index_col=0,sep='\t')
sample_type={}
samplelist = []
def InputData(SampleInfo):
    metadata = pd.read_csv(SampleInfo,header=0,index_col=0,sep='\t')
    filenames = metadata.index 
    adatas = [sc.read_10x_mtx(datadir + filename+ '/filtered_feature_bc_matrix',cache=True) for filename in filenames] # use list to store separate adata
    for i in range(len(adatas)):
        adatas[i].obs['sample'] = metadata.index[i]
        for col in metadata.columns:
            adatas[i].obs[col] = metadata[col][i]
    adata = adatas[0].concatenate(adatas[1:],batch_categories = metadata.index)
    adata.var_names_make_unique()
    sc.pl.highest_expr_genes(adata, n_top=20, save = True )
    sc.pp.filter_cells(adata, min_genes=200)
    sc.pp.filter_genes(adata, min_cells=3)
    
    # ribo_gene_df = pd.read_csv(r'KEGG_RIBOSOME.v2023.1.Hs.csv', sep=',')
    # ribo_gene_df = ribo_gene_df[ribo_gene_df.STANDARD_NAME=='GENE_SYMBOLS']
    # ribo_gene_names = ribo_gene_df.loc[:, 'KEGG_RIBOSOME'].values[0].split(',')
    # adata.var['ribo'] = adata.var.index.isin(ribo_gene_names) 
    
    adata.var['mito'] = adata.var.index.str.startswith(('MT-', 'mt-', 'Mt-'))
    sc.pp.calculate_qc_metrics(adata,
                               expr_type='counts', # default value
                               var_type='genes', # default value
                               qc_vars=['mito'],
                               percent_top=None, 
                               log1p=False, 
                               inplace=True)

    sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mito'],
                 jitter=0.4, multi_panel=True,save=True)
    
    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
    plt.subplots_adjust(wspace=.5)

    sc.pl.scatter(adata, x='total_counts', y='pct_counts_mito',show=False,ax=ax1)
    sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts',show=False,ax=ax2)
    fig.savefig(os.path.join(figure_out, 'scatter.pdf'), 
                dpi=300, bbox_inches='tight')
    return adata

# 计算doublet
sim_doublet_ratio = 2
def Compute_Doublet(adata,rate):
    counts_matrix = adata.to_df()
    n_cells = adata.shape[0]
    scrub = scr.Scrublet(counts_matrix, expected_doublet_rate=rate,
                        n_neighbors = round(0.5 * sqrt(n_cells)),
                        sim_doublet_ratio = sim_doublet_ratio)
    ### annoy-1.15.1
    doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=2, 
                                                          min_cells=3, 
                                                          min_gene_variability_pctl=85, 
                                                          n_prin_comps=30)
    scrub.plot_histogram()
    plt.savefig(os.path.join(figure_out,"histogram.pdf"),dpi=300, bbox_inches='tight')
    plt.show()
    scrub.set_embedding('TSNE', scr.get_tsne(scrub.manifold_obs_, 0.5,10))
    scrub.plot_embedding('TSNE', order_points=True)
    plt.savefig(os.path.join(figure_out,"predicted_doublets.pdf"))
    plt.show()
    
    out_df = pd.DataFrame()
    out_df['barcodes'] = counts_matrix.index
    out_df['doublet_scores'] = doublet_scores
    out_df['predicted_doublets'] = predicted_doublets
    #out_df.to_csv('doublet.txt', index=False,header=True)
    out_df.to_csv(os.path.join(table_out,'doublet.txt'), index=False,header=True)
    out_df.head()
    return out_df,doublet_scores,predicted_doublets 

def Filter_cells(adata,doublet_scores,predicted_doublets,mito):
    adata.obs["doublet_scores"] = doublet_scores
    adata.obs["predicted_doublets"] = predicted_doublets
    #~  可以作为取反的功能
    adata = adata[~adata.obs.predicted_doublets, :]
    
    upper_limit = np.quantile(adata.obs.n_genes_by_counts, 0.98)
    adata = adata[adata.obs.n_genes_by_counts < upper_limit, :]
    adata = adata[adata.obs.pct_counts_mito < mito, :]
    return adata,upper_limit

def CellCycleScoring(adata,cycle,species):
    data = pd.read_csv(cycle,sep=',',header='infer',usecols=[1,2,3,4])
    if species == 'hsa':
        s_genes = list(data.loc[ : ,'hsa_s.genes'])
        g2m_genes  = list(data.loc[ : ,'hsa_g2m.genes'])

    else:
        s_genes = list(data.loc[ : ,'mmu_s.genes'])
        g2m_genes  = list(data.loc[ : ,'mmu_g2m.genes'])
        
    cell_cycle_genes = s_genes + g2m_genes
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    sc.pp.scale(adata,zero_center=False)
    
    sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes, g2m_genes=g2m_genes)
    cell_cycle_genes = [x for x in cell_cycle_genes if x in adata.var_names]
    adata_cc_genes = adata[:, cell_cycle_genes]
    sc.tl.pca(adata_cc_genes)
    sc.pl.pca_scatter(adata_cc_genes, color='phase',save='.Cycle.pdf')
    
    adata.obs['Difference'] = adata.obs['S_score'] - adata.obs['G2M_score']
    adata.write('./write/sce_raw.h5ad')
    return adata



def Run_Normalization(adata,n_neighbors,n_pcs,resolution):
    adata.layers['count'] = adata.X.copy()
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(adata, 
                            n_top_genes=3000,
                            flavor='seurat',
                            subset=False, 
                            batch_key='sample')
    sc.pl.highly_variable_genes(adata,save=True)
    adata = adata[:, adata.var.highly_variable]
    
    sc.pp.regress_out(adata, ['Difference'])
    #scale数据
    sc.pp.scale(adata, max_value=10)
    
    sc.tl.pca(adata)
    sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)
    sc.tl.umap(adata)
    sc.tl.leiden(adata, resolution=resolution, key_added='leiden')                                                                                                                                                 
    
    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
    plt.subplots_adjust(wspace=.5)
    sc.pl.umap(adata, color=['sample'], frameon=False, ax=ax1, show=False)
    sc.pl.umap(adata, color=['leiden'], frameon=False, legend_loc='on data', ax=ax2, show=False)

    plt.savefig(os.path.join(figure_out,"umap.pdf"))
    plt.show()
    return adata

def Run_BBKNN(adata,n_neighbors,n_pcs,resolution):
    sc.tl.pca(adata)
    sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)
    sc.external.pp.bbknn(adata, batch_key='sample')
    sc.tl.umap(adata)
    sc.tl.leiden(adata, resolution=resolution, key_added='leiden')

    fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
    plt.subplots_adjust(wspace=.5)
    sc.pl.umap(adata, color=['sample'], frameon=False, ax=ax1, show=False)
    sc.pl.umap(adata, color=['leiden'], frameon=False, legend_loc='on data', ax=ax2, show=False)

    plt.savefig(os.path.join(figure_out,"umap-BBKNN.pdf"))    
    plt.show()
    
    with rc_context({'figure.figsize': (5, 5)}):
        sc.pl.umap(adata, color='leiden', add_outline=True, legend_loc='on data',
                   legend_fontsize=12, legend_fontoutline=2,frameon=False,
                   title='clustering of cells', palette='Set1',save ='-BBKNN-re.pdf')
    
    adata.write(results_file)
    return adata

def Show_Markers(adata):
    ax = sc.pl.correlation_matrix(adata, 'leiden', figsize=(5,3.5),save= True)
    adata_raw = sc.read_h5ad('./write/sce_raw.h5ad')
    adata_raw.obs = adata.obs
    adata_raw.obsm['X_umap'] = adata.obsm['X_umap']
    adata_raw.obsm['seurat_clusters'] = adata.obsm['leiden']
    adata_raw.obsm['nCount_RNA'] = adata.obsm['total_counts']
    adata_raw.obsm['nFeature_RNA'] = adata.obsm['n_genes']
    adata_raw.obsm['percent.mt'] = adata.obsm['pct_counts_mito']
    
    adata_raw.write('./write/sce_raw.h5ad')
    adata = adata_raw
    
    adata.layers['count'] = adata.X.copy()
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    
    sc.settings.verbosity=2    
    sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')
    markers = sc.get.rank_genes_groups_df(adata, group=None, pval_cutoff=.05, log2fc_min=.25)
    markers.to_csv(os.path.join(table_out,'all_markers.csv'), index=False,header=True)
    
    top5 = pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)
    top5.to_csv(os.path.join(table_out,'top5_markers.csv'),index=False,header =True)
    fig=plt.figure(figsize=(6,24),dpi=100)
    for i in  top5.columns: 
        plt.subplot(2, 7, int(i)+1) #做一个3*3的图 range（9）从0开始，因需要从1开始，所以i+1
        sc.pl.rank_genes_groups_violin(adata, groups=str(i), n_genes=5,show=False)
        plt.tight_layout()
        plt.axis = 'off' #关闭坐标 让图更美观
    plt.savefig(os.path.join(figure_out,"top5-markers.png"))      
    plt.show()
    
    adata.layers['scaled'] = sc.pp.scale(adata, copy=True).X
    sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')
    
    sc.pl.rank_genes_groups_matrixplot(adata, n_genes=3, use_raw=False, vmin=-3, vmax=3, cmap='bwr', layer='scaled',save =True)
    sc.pl.rank_genes_groups_stacked_violin(adata, n_genes=3, cmap='bwr',save = True)
    
    sc.pl.rank_genes_groups_dotplot(adata, n_genes=3, values_to_plot='logfoldchanges', min_logfoldchange=3, vmax=7, vmin=-7, cmap='bwr',save = True)
    sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, use_raw=False, swap_axes=True, show_gene_labels=False,
                                vmin=-3, vmax=3, cmap='bwr',save =True)
    os.system('/PERSONALBIO/work/singlecell/s01/.conda/envs/py4/bin/Rscript ./sceasy.R')
    return adata

    
if __name__ == '__main__':
    start = datetime.datetime.now()
    SampleInfo = './sample_info.txt'
    cycle = './cycle.gene.csv'
    species = 'hsa' ## or mmu
    adata = InputData(SampleInfo)
    out_df,doublet_scores,predicted_doublets = Compute_Doublet(adata,0.06) # 双细胞比率 默认0.06
    adata,upper_limit = Filter_cells(adata,doublet_scores,predicted_doublets,10) # 分辨率
    adata = CellCycleScoring(adata,cycle,species)
    adata = Run_Normalization(adata,50,50,0.99) # neibor pc resolution
    adata = Run_BBKNN(adata,50,50,0.99)
    Show_Markers(adata)
    end = datetime.datetime.now()
    print("程序运行时间："+str((end-start).seconds/3600)+"h")