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=( 5 , 5 )) 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 ( 4062 , 19039 ) ## 过滤后的小提琴分布 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 个细胞。
今天先分享到这里~(未完待续)。