空转第三节课多样本整合的补充2（python版本）

原创

追风少年i

发布于 2023-10-03 20:23:44

5142

发布于 2023-10-03 20:23:44

文章被收录于专栏：空间转录组数据分析空间转录组数据分析

作者，Evil Genius

大家国庆过得如何了？？？？？

如果很开心，不妨分享一下

上一篇文章讲过了，对空间转录组进行整合分析，比较CCA和harmony的结果可以发现

harmony整合的过程中把正常区域整合在了一起，而CCA在整合的过程把部分正常区域和部分肿瘤区域合并成了一个cluster。

CCA对肿瘤区域的聚类结果更为混乱一点，这一点符合认知。

整体效果而言倾向于harmony的分析结果。

关于空间整合分析的结果仍然建议大家采用文献的思路，整合的过程分别使用CCA和harmony进行分析，根据形态学划分来判断采用哪种整合结果。

我们在空间转录组第三节课上的分享的多样本整合基本都是基于R版本，但是对于空间转录组而言很多分析都是基于python的，那么这一篇我们来更新python版本的空间多样本整合。

在这之前了解一下SpaceRanger的分析参数

参数介绍

输出文件夹名字，两个样本的输出文件分别命名为A和P

slide

Visium slide serial number. Refer to the Slide Parameters for information on supported slide versions. Required unless --unknown-slide is passed.

area

Visium capture area identifier. Required unless --unknown-slide is passed. Options for Visium are A1, B1, C1, D1.

loupe-alignment

spaceranger count运行时图片对其有两种方式，一种软件自动识别图片进行对齐，另外一种就是先用Loupe软件手动对齐，生成对于json文件提供给后面的软件并用loupe-alignment参数指定。

测序流程示意图

样本定量实例

#>>>A.sh>>>
human_index_dir=~/DataHub/10X/refdata-gex-GRCh38-2020-A
mouse_index_dir=~/DataHub/10X/refdata-gex-mm10-2020-A
fastqs_dir=~/Project/ST/data/V1_Mouse_Brain_Sagittal_Anterior_Section_1_fastqs
image_path=~/Project/ST/data/V1_Mouse_Brain_Sagittal_Anterior_image.tif
output_dir=~/Project/ST/data

cd ${output_dir}

spaceranger count \
    --id A \
    --description Mouse_Brain_Sagittal_Anterior_Section_1 \
    --transcriptome ${mouse_index_dir} \
    --fastqs ${fastqs_dir} \
    --image ${image_path} \
    --slide V19L29-035 \
    --area B1 \
    --localcores 20 \
    --localmem 128
#<<<A.sh<<<
cd ~/Project/ST/data
nohup zsh A.sh &> A.sh.log &

输出文件

├── filtered_feature_bc_matrix.h5
├── spatial
│   ├── aligned_fiducials.jpg
│   ├── detected_tissue_image.jpg
│   ├── scalefactors_json.json
│   ├── spatial_enrichment.csv
│   ├── tissue_hires_image.png
│   ├── tissue_lowres_image.png
│   ├── tissue_positions.csv
└── web_summary.html

outs文件夹下需要查看的文件

web_summary.html：这个是必须要看的，粗略浏览本次10x样本走SpaceRanger count流程的运行质量
filtered_feature_bc_matrix.h5: Python读取表达量矩阵
spatial：图像信息；文件夹内包括Visium-specific outs: QC images to check image processing pipeline, downsampled input images, and files that describe spot barcode locations in the images

其他输出文件描述可以查看官网介绍https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/output/overview

接下来进行python版本的多样本整合

from pathlib import Path
import re
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import seaborn as sns
import squidpy as sq

import bioquest as bq
import sckit as sk

批量读取数据

adata_dct = {}
for i in Path("data").glob("**/outs"):
  _s = str(i).split('/')[1]
  _a = sc.read_visium(i,library_id=_s)
  _a.var_names_make_unique()
  adata_dct[_s] = _a

adata = sc.concat(adata_dct,label="library_id",uns_merge="unique")
adata.var_names_make_unique()
adata.obs_names_make_unique()
adata
# AnnData object with n_obs × n_vars = 6049 × 32285
#     obs: 'in_tissue', 'array_row', 'array_col', 'library_id'
#     uns: 'spatial'
#     obsm: 'spatial'

过滤，可选（不推荐）

sk.qcMetrics(adata,batch_key="library_id",output_dir=OUTPUT_DIR,mitochondrion=True)
sk.qc_hist(adata,batch_key="library_id",n_genes_by_counts=(0,6000),total_counts=(0,60000),output_dir=OUTPUT_DIR)

gene_bottom,gene_up = np.quantile(adata.obs.n_genes_by_counts.values, [.2,.98])
count_bottom,count_up = np.quantile(adata.obs.total_counts.values, [.2,.98])
afilter = {
    "pct_counts_Mito":"x<30",
    "n_genes_by_counts":f"{gene_bottom}<x<{gene_up}",
    "total_counts":f"{count_bottom}<x<{count_up}"
}
adata = sk.subset(adata,afilter,inplace=False)

基础处理，均一化(normalize),没有进行批次处理

adata = sk.normalise(adata,batch_key='library_id',flavor="Seurat_v3",n_top_genes=3000,output_dir=OUTPUT_DIR,pca_use_hvg=True,n_jobs=24, inplace=False)

聚类

sc.tl.leiden(adata, key_added="Cluster")
sc.tl.tsne(adata,n_jobs=24)
sc.tl.umap(adata,n_jobs=24)
sc.pl.tsne(adata,color="Cluster",legend_loc="on data",show=False,legend_fontoutline=3);
sc.pl.umap(adata, color="Cluster",legend_loc="on data",show=False,legend_fontoutline=3);

空间可视化

fig, axes = plt.subplots(1, 2, figsize=(6, 3))
for x,y in zip(axes,sorted(adata.uns['spatial'].keys())):
    sc.pl.spatial(adata[adata.obs.library_id==y],
        frameon=False,
        color="Cluster",
        size=1.5,
        library_id=y,
        title=y,
        ax=x,
        show=False,
        legend_loc="right margin" if y=='P' else  None);
plt.subplots_adjust(wspace=0.1)

借助harmony去批次

sk.harmony(adata,batch_key='library_id',output_dir=OUTPUT_DIR)

重新聚类

sc.tl.leiden(adata, key_added="Cluster")
sc.tl.tsne(adata, use_rep="X_harmony",n_jobs=24)
sc.tl.umap(adata, use_rep="X_harmony",n_jobs=24)
sc.pl.tsne(adata,color="Cluster",legend_loc="on data",show=False,legend_fontoutline=3);
sc.pl.umap(adata, color="Cluster",legend_loc="on data",show=False,legend_fontoutline=3);

空间可视化

fig, axes = plt.subplots(1, 2, figsize=(6, 3))
for x,y in zip(axes,sorted(adata.uns['spatial'].keys())):
    sc.pl.spatial(adata[adata.obs.library_id==y],
        frameon=False, color="Cluster",
        size=1.5,
        library_id=y,
        title=y,
        ax=x,
        show=False,
        legend_loc="right margin" if y=='P' else  None);
plt.subplots_adjust(wspace=0.1)

adata.write_h5ad(f"{OUTPUT_DIR}/adata_after_harmony.h5ad",compression='lzf')

至于聚类信息python和R如果转换，课上已经讲过了，大家要灵活运用分析方法。

生活很好，有你更好

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据分析

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据分析

#空间转录组

#整合

登录后参与评论

0 条评论

热度