从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()
数据质控包括过滤细胞和过滤基因。
评估一个细胞的质量：总的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')联系客服
