python单细胞学习（二）：使用python批量加速读取count矩阵

python基本部分学完了之后，见专辑《python生信笔记》，就需要开始进行大量的数据分析实战练习进行掌握了。这次给大家找了个数据，文献于2022年12月12日发表在 Cancer Cell 杂志（IF=48.8）上，标题为《High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer》。

这篇文献是一篇极其典范的python流派的数据分析文章，集合了单细胞数据预处理，数据合并，数据整合，数据注释等工作。跨数据集和大数据量！应该有很多细节值得学习，更特别的是注释工作，做这样的一个大数据量的图谱类细胞类型注释，可以看看里面到底是如何处理各种分析细节的！

资源介绍：

python单细胞代码和数据学习的文献分享
python单细胞学习：（IF=48.8）顶刊杂志复现（一）
希望一起学习python单细胞打卡的，加我拉进python学习群，微信：Biotree123，暗号“python学习打卡”，不发暗号不通过好友！！！

前提回顾

书接上回，我们了解了文献的背景，数据，分析方法步骤等信息，今天来学习其中的代码。

作者的核心内容为通过整合19项研究、21个数据集共计298名患者的505份样本（图1A），首次构建了非小细胞肺癌核心图谱。

核心图谱共整合 898,422 个单细胞，将其注释为12个粗粒度细胞类型及44个主要细胞亚型/状态（如分裂期细胞），包括169,223个上皮细胞、670,409个免疫细胞及58,790个基质与内皮细胞（图1B）。

核心图谱对应的数据集合见：文献的Table S1。

梳理好的数据下载链接：https://doi.org/10.5281/zenodo.6411867

可复现学习的代码：https://github.com/icbi-lab/luca

step1-2：见稿子

分析步骤：

QC of the individual datasets based on detected genes, read counts and mitochondrial fractions
Merging of all datasets into a single AnnData object. Harmonization of gene symbols.
Annotation of two "seed" datasets as input for scANVI.
Integration of datasets with scANVI
Doublet removal with SOLO
Annotation of cell-types based on marker genes and unsupervised leiden clustering.
Integration of additional datasets with transfer learning using scArches.

step3：读取数据并构建h5ad对象（11_own_datasets）

作者自己的数据集1里面包括3个patient，提供的为count矩阵：

ls -1  data/11_own_datasets/batch1_3patients/processed/
'NSCLC P1_Lung.csv'
'NSCLC P1_Tumor.csv'
'NSCLC P2_Lung.csv'
'NSCLC P2_Tumor.csv'
'NSCLC P3_Lung.csv'
'NSCLC P3_Tumor.csv'

然后样本表型矩阵在这里：https://github.com/icbi-lab/luca/tree/master/tables

3.1 加载python模块

我这里为conda小环境sc2，使用vscode运行python代码，如果提示你缺少啥就安装啥好了。

import scanpy as sc
import pandas as pd
from glob import glob
from pathlib import Path
import re
import scipy.sparse
from multiprocessing import Pool
import anndata

3.2 读取样本meta信息

这里又学了一个新知识，python读取xlsx文件：

meta = pd.read_excel("tables/patient_table_batch1_3_patients.xlsx", engine="openpyxl", skiprows=1)
meta

修改一下表头：

meta.rename(columns={"Tumornummer": "tumor_id", 
      "Patient": "patient", 
      "Alter": "age", 
      "Geschlecht": "sex", 


    

      "Tumor": "tumor_type"}, inplace=True)
meta["sex"] = [{"W": 'f', "M":'m'}[x] for x in meta["sex"]]
meta

3.3 读取count矩阵

先列出文件：六个样本

filenames = glob("data/11_own_datasets/batch1_3patients/processed/*.csv")
filenames

预定义一个读取count矩阵的函数：

## 定义一个读取count矩阵的函数
def load_counts(filename):
    p = Path(filename)
    ((patient, tissue), ) = re.findall("(P\d+)_(.*)\.csv", p.name)
    expr = pd.read_csv(p, skiprows=5, index_col=0)
    gene_expr = expr.loc[(~expr.index.str.contains("\(Ab\)") & ~expr.index.str.startswith("Lex_")), :]
    ab_expr = expr.loc[expr.index.str.contains("\(Ab\)"), :]
    obs = pd.DataFrame().assign(cell_id=expr.columns)
    obs["patient"] = patient
    obs["tissue"] = tissue
    obs.set_index("cell_id", inplace=True)
    adata = sc.AnnData(X=scipy.sparse.csc_matrix(gene_expr).T, obs=obs)
    adata.var_names = gene_expr.index
    adata.obsm["surface_protein"] = ab_expr.T
    return adata

这里面有几句python代码我不是很明白意思，然后让deepseek帮我解读一下：

这段Python代码使用了多进程并行的方式来批量加载数据文件

# Pool(16)创建包含16个工作进程的进程池
# 数字16表示同时运行的进程数量
# map函数将任务分配给所有进程并行执行
# 存储所有进程返回结果的列表
with Pool(16) as p: 
    adatas = p.map(load_counts, filenames)

厉害哇！又学到一个知识点！

3.4 创建h5ad对象

将上面的adatas合并在一起，生成h5ad对象：

adata = anndata.concat(adatas, index_unique="_", join="outer")
adata
# 修改一下meta信息
meta["patient"] = [x.strip() for x in meta["patient"]]
adata.obs = adata.obs.reset_index().merge(meta, on=["patient"], how="left").set_index("cell_id")
adata
adata.obs

作者还修改了一点细节：

adata.obs.drop_duplicates()
adata.obs["condition"] = "NSCLC"
adata.obs["origin"] = ["tumor_primary" if c == "Tumor" else "normal_adjacent" for c in adata.obs["tissue"]]
adata.obs["sample"] = [f"{patient}_{origin}" for patient, origin in zip(adata.obs["patient"], adata.obs["origin"])]
adata.obs["sex"] = [{"m": "male", "f": "female"}[s] for s in adata.obs["sex"]]
adata.obs["tissue"] = "lung"
adata.obs.drop_duplicates()

最后保存刚刚创建好的h5ad数据：

# 保存
adata.write_h5ad("batch1_3patients.h5ad", compression="lzf")

python起飞的一天！

如果上面的内容对你有帮助，欢迎一键三连！

友情转发：