单细胞转录组实战02: 数据整理与之质控

整理数据

从cellranger的输出目录中读取filtered_feature_bc_matrix.h5表达量矩阵，并把多个样本合并为1个anndata对象。

导库建立工作目录

from pathlib import Path
import scanpy as sc

OUTPUT_DIR='output/01.preprocess'
Path(OUTPUT_DIR).mkdir(parents=True,exist_ok=True)

adata_dict = {}
for i in Path('quantify').glob("**/filtered_feature_bc_matrix.h5"):
    adata = sc.read_10x_h5(i)
    adata.var_names_make_unique()
    adata.obs_names_make_unique()
    adata_dict[str(i.parent).split('/')[1]] = adata
    
adata = sc.concat(adata_dict,label='Sample',axis=0)
adata.obs_names_make_unique()

Quality control

数据质控包括过滤细胞和过滤基因。

• 过滤细胞

评估一个细胞的质量：总的UMIs数量, 基因数量, 线粒体比例

• 过滤基因

一个基因至少在几个细胞中表达

简单的过滤一下细胞和基因，一个细胞至少表达300个基因，一个基因至少在10个细胞中表达。

sc.pp.filter_cells(adata,min_genes=300)
sc.pp.filter_genes(adata,min_cells=10)

计算线粒体基因的比例，线粒体基因以MT开头

adata.var['Mito'] = adata.var_names.str.startswith(r'MT-',r"mt-")
sc.pp.calculate_qc_metrics(adata, qc_vars=['Mito'], percent_top=None, log1p=False, inplace=True)

可视化count数量、基因数量、线粒体基因百分比

ks = ('total_counts','n_genes_by_counts','pct_counts_Mito')
_, axes = plt.subplots(1, 3, figsize=(6, 3))
axes = axes.flatten()
for a,k in enumerate(ks):
    sc.pl.violin(adata, keys=k,jitter=False,show=False,ax=axes[a],ylabel='')
plt.subplots_adjust(wspace = 0.5);

每个样本的count数量、基因数量、线粒体基因百分比

_, axes = plt.subplots(1, 2, figsize=(6, 3))
for a,k in enumerate(ks[1:]):
    sc.pl.scatter(adata, x='total_counts', y=k,color="pct_counts_Mito",ax=axes[a],legend_loc="right margin",show=False)
plt.subplots_adjust(wspace = 0.3);

_, axes = plt.subplots(1, 3, figsize=(6, 3))
axes = axes.flatten()
for a,k in enumerate(ks):
    sc.pl.violin(adata,keys=k,groupby="Sample",rotation=45,jitter=False,show=False,stripplot=False,ax=axes[a])
plt.subplots_adjust(wspace = 0.5);

细胞中有较高的线粒体比例有两种可能，1是该细胞受到刺激或者是正在死亡的细胞，2是细胞活性比较强。

• 简单的过滤

线粒体基因比例小于20%，n_genes_by_counts小于整体的98%，total_counts_upper小于整体的98%

n_genes_by_counts_upper_lim = np.quantile(adata.obs.n_genes_by_counts.values, .98)
total_counts_upper_lim = np.quantile(adata.obs.total_counts.values, .98)
afilter = {
  "pct_counts_Mito":"x<20",
  "n_genes_by_counts":f"x<{n_genes_by_counts_upper_lim}",
  "total_counts":f"x<{total_counts_upper_lim}"
}
for k in afilter:
    v = afilter.get(k)
    _lg = adata.obs[k].apply(lambda x: eval(v))
    adata = adata[_lg, :]

数据标准化

可能有多种原因使得每个细胞在测序时候总的count不同，为了细胞间可以比较因此将他们都统一为10000。

adata.layers["counts"] = adata.X.copy() # preserve counts

sc.pp.normalize_total(adata, target_sum=10000)
sc.pp.log1p(adata,base=None)

筛选高变基因

一般筛选1000-5000个高变基因（HVG），可以降低数据维度。

HVG用于降维、聚类及其可视化等。

同时也在adata存一份全部基因的rawdata用于marker鉴定、差异分析、细胞周期分数计算、基因表达量可视化等。

sc.pp.highly_variable_genes(adata, n_top_genes = 3000,
    flavor="seurat_v3", subset=False, batch_key="Sample",layer='counts')
sc.pl.highly_variable_genes(adata,log=True)
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

regress_out用普通的线性回归消除pct_counts_Mito变量的影响。

scale可选，把数据中心化。

sc.pp.regress_out(adata, keys='pct_counts_Mito', n_jobs=12)
sc.pp.scale(adata, max_value=10)

批次校正

首先经过PCA降维然后harmony改变维度信息，从而达到批次校正目的，因此adata.X中的表达量是不会变的。

先展示一下批次效应PCA图

sc.tl.pca(adata, use_highly_variable=True, svd_solver='arpack',n_comps=50)
sc.pl.pca(adata, color="Sample",legend_loc="on data",legend_fontsize="small")
sc.pl.pca_variance_ratio(adata,n_pcs=50)

再展示一下批次效应UMAP图

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.pl.umap(adata, color="Sample",legend_loc="on data",legend_fontsize="small")

最后批次校正，看一下整合后的UMAP图

sc.external.pp.harmony_integrate(adata,key='Sample',adjusted_basis='X_pca')
sc.pp.neighbors(adata,n_neighbors=15)
sc.tl.umap(adata)
sc.pl.umap(adata, color='Sample',legend_loc="on data",legend_fontsize="small")

保存数据

adata.write_h5ad(OUTPUT_DIR+"/adata.h5",compression='lzf')