生信小白学单细胞转录组(sc-RNA)测序数据分析——python语言

最新推荐文章于 2024-06-15 09:40:43 发布

鹓雏历天

最新推荐文章于 2024-06-15 09:40:43 发布

阅读量1.1k

点赞数 9

文章标签：数据分析 python 开发语言

本文链接：https://blog.csdn.net/qq_46397705/article/details/136743069

版权

详细流程已在生信小白学单细胞转录组(sc-RNA)测序数据分析——R语言一文中解读过

import package

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

import matplotlib as mpl
import matplotlib.pyplot as plt

mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams["font.sans-serif"] = "Arial"

%config InlineBackend.figure_format = 'retina'

sc.settings.set_figure_params(dpi=50, dpi_save=300, figsize=(4, 4))

loading data

输入文件类型：mtx tsv tsv

#第一种方法sc.read_10x_mtx（10x的三个标准文件 mtx, tsv, tsv）
adata = sc.read_10x_mtx('data/filtered_gene_bc_matrices/', var_names='gene_symbols', cache=True)
# adata.var_names_make_unique()
adata

#第二种方法，自己构造
obs = pd.read_csv('./data/filtered_gene_bc_matrices/barcodes.tsv', header=None, index_col=0, sep='\t')
var = pd.read_csv('./data/filtered_gene_bc_matrices/genes.tsv', header=None, index_col=1, sep='\t')
obs.index.name = ''
var.index.name = ''
var.columns = ['gene_ids']


from scipy.io import mmread
from scipy.sparse import csr_matrix
mtx = mmread('./data/filtered_gene_bc_matrices/matrix.mtx')
mtx = mtx.T
mtx = csr_matrix(mtx)
adata1 = sc.AnnData(mtx, obs=obs, var=var)
#去除adata1的重复值
adata1.var_names_make_unique()

输入文件类型：h5ad

adata_PBS = sc.read_10x_h5("./GSM5067314_PBS_filtered_feature_bc_matrix.h5", genome=None, gex_only=True)
adata_TNF = sc.read_10x_h5("./GSM5067313_TNF_filtered_feature_bc_matrix.h5", genome=None, gex_only=True)

#参数说明
#'./GSM5067314_PBS_filtered_feature_bc_matrix.h5'：文件路径，指定了要读取的10x Genomics数据集文件的位置
#'genome=None'：指定要使用的基因组信息,设置为None表示不使用任何特定的基因组信息
#'gex_only=True'：指示是否仅读取基因表达数据，设置为True表示只读取基因表达数据，而不读取其他类型的数据，如蛋白质表达数据

QC

sc.pp.filter_cells(adata, min_genes=200)
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)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)

在这里插入图片描述

sc.settings.set_figure_params(dpi_save=300)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', ax=ax[0], show=False)
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', ax=ax[1], show=False)
plt.subplots_adjust(wspace=.4)

在这里插入图片描述

adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]

Normalization

#使用移位对数的方法进行标准化
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

log标准化前后的adata.X数据

feature selection

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

在这里插入图片描述

adata.raw = adata
#备份adata到adata.raw
adata = adata[:, adata.var.highly_variable]
#slice选取adata的前两千个高变基因

#scale标准化，消除测序深度的影响，为降维聚类做准备
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
sc.pp.scale(adata, max_value=10)

reduction 降维

sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True)

在这里插入图片描述

sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50

在这里插入图片描述

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

Clustering

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

在这里插入图片描述

markers = ["MS4A1", "TYROBP", "CD14",'FCGR3A', "FCER1A", "CCR7", "IL7R", "PPBP", "CD8A"]
sc.pl.umap(adata, color=markers, ncols=3)

在这里插入图片描述

new_cluster_names = ['Naive CD4 T', 'Memory CD4', 'CD14 Monocytes','B', 'CD8 T', 'FCGR3A Monocytes','NK', 'DC', 'Platelet']
adata.rename_categories('leiden', new_cluster_names)

sc.settings.set_figure_params(dpi=50, dpi_save=300, figsize=(7, 7))
sc.pl.umap(adata, color='leiden', legend_loc='on data')

在这里插入图片描述

sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

在这里插入图片描述

Understand the data structure of Anndata

In [45]:
adata
Out[45]:
AnnData object with n_obs × n_vars = 2638 × 2000
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'#细胞信息
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'#基因信息
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'leiden', 'leiden_colors', 'rank_genes_groups'#非结构化数据
    obsm: 'X_pca', 'X_umap'#细胞降维信息
    varm: 'PCs'#基因降维信息
    obsp: 'distances', 'connectivities'#距离矩阵

In [46]:
adata_norm = adata.raw.to_adata()
In [47]:
adata_norm
Out[47]:
AnnData object with n_obs × n_vars = 2638 × 13714
    obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'
    var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'leiden', 'leiden_colors', 'rank_genes_groups'
    obsm: 'X_pca', 'X_umap'
    obsp: 'distances', 'connectivities'