首页 > 文章列表 > 2025年最后一个python学习笔记！明年见~

2025年最后一个python学习笔记！明年见~

📅 发布于 2026年03月10日 21:28

python基本部分学完了之后，见专辑《python生信笔记》，就需要开始进行大量的数据分析实战练习进行掌握了。今天来学习一篇顶刊杂志的单细胞数据处理，使用python！

数据背景

Goveia_Carmeliet_2020_NSCLC 数据处理

今天来分析直播的那篇文献里面里面提到的 Goveia, Carmeliet et al. 数据，来自文献 10.1016/j.ccell.2019.12.001 。文献标题为《An Integrated Gene Expression Landscape Profiling Approach to Identify Lung Tumor Endothelial Cell Heterogeneity and Angiogenic Candidates》。杂志：Cancer Cell. 2020 年1月1号。

数据集为：data/12_input_adatas/goveia_carmeliet_2020_nsclc.h5ad，在文献里面的下载地址为 ：https:///records/7227571，在里面的 input_data.tar.xz

主要结果：

本研究首次在单细胞层面系统绘制了肺肿瘤内皮细胞的跨物种、跨模型异质性图谱，揭示了新型免疫调节性TEC亚群及尖端TEC的功能分化（迁移型 vs 基底膜重塑型）。

  • 对来自人类、小鼠及培养肺肿瘤模型的56,771个内皮细胞进行了单细胞RNA测序
  • Tip ECs 尖端内皮细胞可区分为迁移型与基底膜重塑型两种表型
  • Capillary and venous ECs 毛细血管与静脉内皮细胞表达免疫调节的 gene signatures
  • 整合分析鉴定出胶原修饰作为血管生成的关键通路

人非小细胞肺癌内皮细胞的单细胞图谱

六大内皮：arterial, capillary, venous, lymphatic, tip, and immature stalk phenotypes; listed in Table S5

数据读取

先导入模块：

# %%
# %load_ext autoreload
# %autoreload 2
import scanpy as sc
import numpy as np
import pandas as pd
import scipy.sparse as sp
from scanpy_helpers.annotation import AnnotationHelper
from nxfvars import nxfvars

sc.set_figure_params(figsize=(5, 5))

读取数据：

对应的文件：goveia_carmeliet_2020_nsclc.h5ad

# %%
adata = sc.read_h5ad(filename="data/12_input_adatas/goveia_carmeliet_2020_nsclc.h5ad" )
adata

10w+的细胞：

AnnData object with n_obs × n_vars = 100512 × 18450
    obs: 'patient', 'condition', 'origin', 'Endothelial cell', 'Cluster', 'tissue', 'sample'

其他的属性

看看其他的属性了解数据情况：

adata.obs
adata.obs['patient'].value_counts()
adata.obs['condition'].value_counts()
adata.obs['origin'].value_counts()
adata.obs['Endothelial cell'].value_counts()
adata.obs['Cluster'].value_counts()
adata.obs['sample'].value_counts()

可以看到数据里面有7w多的非内皮细胞：

数据QC

看看qc前的指标：

#表达矩阵里的数值范围
np.min(adata.X),np.max(adata.X)
#过滤前
adata.X.shape
adata[0:6, ['CD3D','TCL1A','MS4A1']].to_df()
adata.var.head()
adata.var_names[:10]
adata.obs.head()

# 计算线粒体基因比例
adata.var['mt']=adata.var_names.str.startswith('MT-')
# 计算核糖体基因比例
adata.var['ribo']=adata.var_names.str.match('^RP[SL]')
# 计算红血细胞基因比例
adata.var['hb']=adata.var_names.str.match('^HB[^(P)]')
# 计算
sc.pp.calculate_qc_metrics(adata,qc_vars=['mt','ribo','hb'],percent_top=None,log1p=False,inplace=True)

# 质控
adata.obs.head()

# 可视化细胞的上述比例情况
# jitter：可以设置散点抖动的宽度
# rotation：可以调整x轴标签的旋转角度，样本多可以设置90度
# size：设置散点的大小
from matplotlib.pyplot import rc_context
with rc_context({"figure.figsize": (11, 10)}):
    sc.pl.violin(
        adata,
        ['n_genes_by_counts','total_counts','pct_counts_mt'],
        groupby="sample",
        jitter=0.4,rotation=90,multi_panel=True,
        #stripplot=False,  # remove the internal dots
        #inner="box",  # adds a boxplot inside violins
    )

直播文献的过滤标准：

这个标准感觉有点宽松。

过滤一下：

## 过滤前的数据
adata.shape

## 过滤
# 方法2：一次性过滤所有条件（更高效）
adata = adata[
    (adata.obs['n_genes_by_counts'] > 200) & 
    (adata.obs['n_genes_by_counts'] < 6000) &
    (adata.obs['total_counts'] > 500) & 
    (adata.obs['total_counts'] < 30000) &
    (adata.obs['pct_counts_mt'] < 20), :
].copy()

## 过滤后的数据
adata.shape
# (61913, 18450)

过滤了4w走了，这。

过滤后的小提琴看看：

from matplotlib.pyplot import rc_context
with rc_context({"figure.figsize": (9, 10)}):
    sc.pl.violin(
        adata,
        ['n_genes_by_counts','total_counts','pct_counts_mt'],
        groupby="sample",
        jitter=0.4,rotation=90,multi_panel=True,
        #stripplot=False,  # remove the internal dots
        #inner="box",  # adds a boxplot inside violins
    )
    
## 保存数据
adata.write("01_qc/goveia_carmeliet_2020_nsclc_qc.h5ad")

降维聚类

# 首先将数据矩阵标准化（校正文库大小）：
sc.pp.normalize_total(adata,target_sum=1e4)

# 对数据进行log
sc.pp.log1p(adata)
adata.raw = adata.copy()

# 高变基因
sc.pp.highly_variable_genes(adata)

# 归一化
sc.pp.scale(adata)

# pca分析
sc.tl.pca(adata,use_highly_variable=True)

# harmony整合
sc.external.pp.harmony_integrate(adata,'sample')

# 聚类
sc.pp.neighbors(adata,use_rep='X_pca_harmony',n_pcs=30)
sc.tl.umap(adata)
# Using the igraph implementation and a fixed number of iterations can be significantly faster, especially for larger datasets
sc.tl.leiden(adata, flavor="igraph", n_iterations=2,resolution=0.5)

可视化看看：

# umap图
from matplotlib.pyplot import rc_context
with rc_context({"figure.figsize": (6.5, 6)}):
    sc.pl.umap(adata, color=["leiden",'sample'], legend_loc="on data", size=7)

整合结果还可以：

画一下marker看看：

sc.pl.umap(adata, color=["EPCAM", "PTPRC","NKG7", "CD3D","VWF","ACTA2","DCN","CD79A","CD68"])

各个主要的marker也出来了，数据质量还行。

最后保存一下数据：

## 保存数据
adata.write("01_qc/goveia_carmeliet_2020_nsclc_umap.h5ad")

今天分享到这里，祝大家元旦快乐！

明年见~