Py学习  »  Python

python单细胞学习(四):lambrechts_2018_luad 数据读取与质控

生信技能树 • 6 天前 • 31 次点击  

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流派的数据分析文章,集合了单细胞数据预处理,数据合并,数据整合,数据注释等工作。跨数据集和大数据量!应该有很多细节值得学习,更特别的是注释工作,做这样的一个大数据量的图谱类细胞类型注释,可以看看里面到底是如何处理各种分析细节的!

note:最近这篇文献我们会开一个直播进行实战训练,进群方式:python单细胞数据分析实战直播交流群

如果下面的代码对你来说有困难,我们生信技能树每月都有一期针对0基础的生信入门培训班,12月份开始招生啦,快来瞧一瞧,详细介绍点击这里生信入门&数据挖掘线上直播课12月班

给你最好的答疑!

前提回顾

作者的核心内容为通过整合19项研究、21个数据集共计298名患者的505份样本(图1A),首次构建了非小细胞肺癌核心图谱。核心图谱共整合 898,422 个单细胞,将其注释为12个粗粒度细胞类型及44个主要细胞亚型/状态(如分裂期细胞),包括169,223个上皮细胞、670,409个免疫细胞及58,790个基质与内皮细胞(图1B)。

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

前面已有内容

step0-2:见稿子

step0:熟悉数据和代码。这一步看文献和下载作者的代码。文献背景介绍:(IF=48.8)顶刊杂志的带有超详细python单细胞代码和数据学习的文献分享

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

step1:下载作者提供的数据并探索。这一步下载好作者的数据。

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

step2:看作者的分析方法细节。这里了解作者的分析思路。步骤见:

  • 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),见稿子

这里得到文件 batch1_3patients.h5ad。

作者已经有提供好的,在下载的数据目录中翻一下就有:data/11_own_datasets/batch1_3patients/

step4:batch2数据处理(11_own_datasets),见稿子

这里会得到文件 ukim_v_batch2.h5ad。

step5:处理 Lambrechts_2018_LUAD_6653 数据

数据列表在文献的 table s1,前去下载即可。

前面作者自测的两个数据都整理成了h5ad格式,这一步整理 Lambrechts_2018_LUAD_6653 数据,来自 2018年8月发表在 Nat Med 杂志上的文献,标题为《Phenotype molding of stromal cells in the lung tumor microenvironment》。DOI:10.1038/s41591-018-0096-5。

作者这里提供了初步的h5ad数据,需要读取进来做一点处理。

看一下这个数据的处理方式

这个数据对应三个h5ad文件:这三个文件我们可以在 作者提供的  input_data.tar.xz 文件里面可以找到。

每个sample对应的各个指标:

每个h5ad对应的过滤参数如下:

加载下面的模块:

环境准备:见这里

# %%
# %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=(55))

lambrechts_2018_luad_6653.h5ad质控

这个数据包含 36589个细胞 patient 678。

1.读取数据:

# %%
#"../../data/20_integrate_scrnaseq_data/02_qc_and_filtering/Lambrechts_2018_LUAD_6653/Lambrechts_2018_LUAD_6653.qc.h5ad",
adata1 = sc.read_h5ad( nxfvars.get( "adata_qc","./data/12_input_adatas/lambrechts_2018_luad_6653.h5ad") )
adata1

adata1的内容如下:

AnnData object with n_obs × n_vars = 36589 × 36601
    obs: 'dataset''sample''patient''origin''condition''platform''age''sex''tissue'

已有的细胞信息:

adata1.obs

探索一下:

# 样本数
adata1.obs['sample'].nunique()
adata1.obs['sample'].value_counts()
# patients 数
adata1.obs['patient'].value_counts()

其他:

adata1.obs['origin'].value_counts()
adata1.obs['tissue'].value_counts()
adata1.obs['condition'].value_counts()
adata1.obs['platform'].value_counts()

2.计算三个qc指标

这里大家应该非常熟悉,不再赘述:

#质控
# annotate the group of mitochondrial genes as "mt"
adata1.var["mt"] = adata1.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata1, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
adata1.obs.head()

未过滤前的qc小提琴图:

## 未过滤前的小提琴分布
sc.pl.violin(
    adata1,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

3.过滤

根据作者给的指标进行过滤:

# 基本过滤
sc.pp.filter_cells(adata1, min_genes=250)
sc.pp.filter_genes(adata1, min_cells=3)

## 过滤
# 方法2:一次性过滤所有条件(更高效)
adata1 = adata1[
    (adata1.obs['n_genes_by_counts'] > 250) & 
    (adata1.obs['n_genes_by_counts'] 10000) &
    (adata1.obs['total_counts'] > 1200) & 
    (adata1.obs['total_counts'] 40000) &
    (adata1.obs['pct_counts_mt'] 20), :
].copy()

过滤后的小提琴:




    
## 过滤后的数据
adata1.shape
# (26511, 25514)

## 过滤后的小提琴分布
sc.pl.violin(
    adata1,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

4.保存数据

过滤后还有 26511 个细胞,过滤掉了  36589 -  26511 = 10078。这里确实不太知道作者那个最小total_counts为什么设置为 1200。过滤掉的细胞有点多!

## 保存数据
adata1.write("lambrechts_2018_luad_6653_qc.h5ad")

lambrechts_2018_luad_6149v1.h5ad 质控

这个数据包含,lambrechts_2018_luad_6149v1.h5ad:4994个细胞 patient 1 2

这里跟上面一样的流程,全部代码放在一起:

## 1.读取数据
adata_v1 = sc.read_h5ad( nxfvars.get( "adata_qc","./data/12_input_adatas/lambrechts_2018_luad_6149v1.h5ad") )
adata_v1
#  4994 × 36601
adata_v1.obs
adata_v1.obs['sample'].unique()
adata_v1.obs['sample'].value_counts()
adata_v1.obs['patient'].value_counts()
adata_v1.obs['platform'].value_counts()

## 2.计算三个qc指标
# 质控
# annotate the group of mitochondrial genes as "mt"
adata_v1.var["mt"] = adata_v1.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata_v1, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
adata_v1.obs.head()

## 未过滤前的小提琴分布
sc.pl.violin(
    adata_v1,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

## 3.过滤
# 基本过滤
sc.pp.filter_cells(adata_v1, min_genes=200)
sc.pp.filter_genes(adata_v1, min_cells=3)

## 过滤
# 方法2:一次性过滤所有条件(更高效)
adata_v1 = adata_v1[
    (adata_v1.obs[ 'n_genes_by_counts'] > 200) & 
    (adata_v1.obs['n_genes_by_counts'] 10000) &
    (adata_v1.obs['total_counts'] > 600) & 
    (adata_v1.obs['total_counts'] 30000) &
    (adata_v1.obs['pct_counts_mt'] 15), :
].copy()

## 过滤后的数据
adata_v1.shape
(406219039)

## 过滤后的小提琴分布
sc.pl.violin(
    adata1,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

## 保存数据
adata_v1.write("lambrechts_2018_luad_6149v1_qc.h5ad")

未过滤前的小提琴分布:

过滤后的小提琴分布:

过滤后还有 4062 个细胞。

总共过滤细胞数:4994 - 4062 = 932个细胞。

lambrechts_2018_luad_6149v2.h5ad 质控

这个数据包含,lambrechts_2018_luad_6149v2.h5ad:59562个细胞 3 4 5

这里跟上面一样的流程,全部代码放在一起:

## 1.读取数据
adata_v2 = sc.read_h5ad( nxfvars.get( "adata_qc","./data/12_input_adatas/lambrechts_2018_luad_6149v1.h5ad") )
adata_v2
#  AnnData object with n_obs × n_vars = 59562 × 36601
adata_v2.obs
adata_v2.obs['sample'].unique()
adata_v2.obs['sample'].value_counts()
adata_v2.obs['patient'].value_counts()
adata_v2.obs['platform'].value_counts()

## 2.计算三个qc指标
# 质控
# annotate the group of mitochondrial genes as "mt"
adata_v2.var["mt"] = adata_v2.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata_v2, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)
adata_v2.obs.head()

## 未过滤前的小提琴分布
sc.pl.violin(
    adata_v2,
    ["n_genes_by_counts""total_counts" "pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

## 3.过滤
# 基本过滤
sc.pp.filter_cells(adata_v2, min_genes=250)
sc.pp.filter_genes(adata_v2, min_cells=3)

## 过滤
# 方法2:一次性过滤所有条件(更高效)
adata_v2 = adata_v2[
    (adata_v2.obs['n_genes_by_counts'] > 250) & 
    (adata_v2.obs['n_genes_by_counts'] < 10000) &
    (adata_v2.obs['total_counts'] > 600) & 
    (adata_v2.obs['total_counts'] < 30000) &
    (adata_v2.obs['pct_counts_mt'] < 20), :
].copy()

## 过滤后的数据
adata_v2.shape
# (48397, 26927)

## 过滤后的小提琴分布
sc.pl.violin(
    adata1,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True
)

## 过滤的小提琴分布
sc.pl.violin(
    adata_v2,
    ["n_genes_by_counts""total_counts""pct_counts_mt"],
    jitter=0.4,
 groupby="sample",
   # **{"color": "#f8766d"},  # 使用额外的参数来设置颜色
    multi_panel=True,
 rotation=90  # 旋转x轴标签90度
)

## 保存数据
adata_v2.write("lambrechts_2018_luad_6149v2_qc.h5ad")

未过滤前的小提琴分布:

过滤后:

过滤后的细胞数还有  48397。

过滤掉的细胞数:59562 - 48397 = 11165个细胞。

总共获得 26511 + 4062 + 48397 =  78970 个细胞。

今天先分享到这里~(未完待续)。

友情转发:


Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/190086