【论文复现】基于多模态深度学习方法的单细胞多组学数据聚类（【生物信息学】实验二：多组学数据融合：scMDC）

一、实验介绍

1. 论文：基于多模态深度学习方法的单细胞多组学数据聚类

Abstract

2. Github链接

二、实验环境

0. 作者要求

1. 环境复现

实验一

实验二（本实验）

2. 库版本介绍

实验一

实验二

3. IDE

三、实验内容

1. 用法

2. 输出

3. 参数

4. run_scMDC

设置超参数

对scRNA-seq数据进行预处理

构建scMultiCluster模型

预训练

使用KMeans确定聚类数k

微调模型,训练聚类层

保存预测结果和embedding到文件

输出预测结果y_pred,计算指标AMI、NMI和ARI

5. 聚类结果可视化

可视化结果

附录：run_scMDC训练过程

一、实验介绍

1. 论文：基于多模态深度学习方法的单细胞多组学数据聚类

Clustering of single-cell multi-omics data with a multimodal deep learning method | Nature Communicationshttps://www.nature.com/articles/s41467-022-35031-9

Abstract

Single-cell multimodal sequencing technologies are developed to simultaneously profile different modalities of data in the same cell. It provides a unique opportunity to jointly analyze multimodal data at the single-cell level for the identification of distinct cell types. A correct clustering result is essential for the downstream complex biological functional studies. However, combining different data sources for clustering analysis of single-cell multimodal data remains a statistical and computational challenge. Here, we develop a novel multimodal deep learning method, scMDC, for single-cell multi-omics data clustering analysis. scMDC is an end-to-end deep model that explicitly characterizes different data sources and jointly learns latent features of deep embedding for clustering analysis. Extensive simulation and real-data experiments reveal that scMDC outperforms existing single-cell single-modal and multimodal clustering methods on different single-cell multimodal datasets. The linear scalability of running time makes scMDC a promising method for analyzing large multimodal datasets.

开发单细胞多模态测序技术以同时分析同一细胞中的不同数据模式。它提供了一个独特的机会，可以在单细胞水平上联合分析多模态数据，以鉴定不同的细胞类型。正确的聚类结果对于下游复杂的生物学功能研究至关重要。然而，组合不同的数据源对单细胞多模态数据进行聚类分析仍然是一个统计和计算挑战。在这里，我们开发了一种新颖的多模态深度学习方法scMDC，用于单细胞多组学数据聚类分析。scMDC 是一个端到端的深度模型，它明确表征不同的数据源，并共同学习深度嵌入的潜在特征以进行聚类分析。大量的仿真和真实数据实验表明，scMDC在不同的单细胞多模态数据集上优于现有的单细胞单模态和多模态聚类方法。运行时的线性可扩展性使scMDC成为分析大型多模态数据集的有前途的方法。

2. Github链接

GitHub - xianglin226/scMDC：单细胞多组学深度聚类https://github.com/xianglin226/scMDC

二、实验环境

0. 作者要求

Python 3.8.1

Pytorch 1.6.0

Scanpy 1.6.0

SKlearn 0.22.1

Numpy 1.18.1

h5py 2.9.0

本研究中scMDC的所有实验都是在Nvidia Tesla P100（16G）GPU上进行的。
我们建议在 conda 环境中安装 conda 环境（conda create -n scMDC）。
scMDC 对包含 5000 cells的数据集进行聚类大约需要 3分钟。

1. 环境复现

未重新配置环境，继续使用前文深度学习系列文章的环境，实践表明可行：

conda create -n DL python=3.7

conda activate DL

pip install torch==1.8.1+cu102 torchvision==0.9.1+cu102 torchaudio==0.8.1 -f https://download.pytorch.org/whl/torch_stable.html

conda install matplotlib

conda install scikit-learn

实验一

conda install pandas

conda install seaborn

conda install networkx

conda install statsmodels

pip install pyHSICLasso

注：本人的实验环境按照上述顺序安装各种库，若想尝试一起安装（天知道会不会出问题）

实验二（本实验）

pip install scanpy

2. 库版本介绍

自己配置的py3.7环境各版本介绍

软件包	本实验版本	目前最新版
matplotlib	3.5.3	3.8.0
numpy	1.21.6	1.26.0
python	3.7.16
scikit-learn	0.22.1	1.3.0
torch	1.8.1+cu102	2.0.1
torchaudio	0.8.1	2.0.2
torchvision	0.9.1+cu102	0.15.2

实验一

networkx	2.6.3	3.1
pandas	1.2.3	2.1.1
pyHSICLasso	1.4.2	1.4.2
seaborn	0.12.2	0.13.0
statsmodels	0.13.5	0.14.0

实验二

Installing collected packages: stdlib_list, natsort, llvmlite, h5py, session-info, numba, pynndescent, anndata, umap-learn, scanpy
Successfully installed anndata-0.8.0 h5py-3.8.0 llvmlite-0.39.1 natsort-8.4.0 numba-0.56.4 pynndescent-0.5.10 scanpy-1.9.3 session-info-1.0.0 stdlib_list-0.9.0 umap-learn-0.5.4

3. IDE

建议使用Pycharm

win11 安装 Anaconda（2022.10）+pycharm（2022.3/2023.1.4）+配置虚拟环境_QomolangmaH的博客-CSDN博客https://blog.csdn.net/m0_63834988/article/details/128693741https://blog.csdn.net/m0_63834988/article/details/128693741https://blog.csdn.net/m0_63834988/article/details/128693741https://blog.csdn.net/m0_63834988/article/details/128693741

三、实验内容

1. 用法

准备 h5 格式的输入数据。（参阅“数据集”文件夹中的自述文件）
根据“脚本”文件夹中的运行脚本运行scMDC（如果您处理mRNA+ATAC数据并使用run_scMDC_batch.py进行多批次数据聚类，请注意参数设置）
基于训练良好的 scMDC 模型通过run_LRP.py运行 DE 分析（请参阅“脚本”文件夹中的 LRP 运行脚本）

2. 输出

scMDC输出数据的潜在表示，可用于进一步的下游分析，并通过t-SNE或Umap可视化;
多批次 scMDC 输出集成数据集的潜在表示，在其上校正批处理效应。
LRP输出一个基因等级，指示基因对给定簇的重要性，可用于通路分析。

3. 参数

--n_clusters：簇数（K）;如果此参数设置为 -1，scMDC 将估计 K。
--cutoff：一个周期的比率，在此之前，模型只训练低级自动编码器。
--batch_size：批量大小。
--data_file：数据输入的路径。
- 数据格式：H5。
- 结构：X1（RNA），X2（ADT或ATAC），Y（标签，如果退出），批次（多批次数据聚类的批次指示符）。
--maxiter：训练的最大时期。默认值：10000。
--pretrain_epochs：预训练的周期数。默认值：400。
--gamma：聚类损失系数。默认值：0.1。
--phi1 和 phi2：预训练和聚类阶段的 KL 损失系数。默认值：CITE-Seq 为 0.001;0.005 表示 SMAGE-Seq*。
--update_interval：检查性能的间隔。默认值：1。 --tol：停止模型的条件，即更改标签的百分比。
--tol：0.001。
--ae_weights：权重文件的路径。
--save_dir：存储输出的目录。
--ae_weight_file：存储权重的目录。
--resolution：用于估计 k 的分辨率参数，默认值：0.2。
--n_neighbors：估计 K 的n_neighbors参数，默认值：30。
--embedding_file：如果保存嵌入文件。默认值：否
--prediction_file：如果保存预测文件。默认值：否

--encodeLayer：RNA 的低电平编码器层：默认值：CITE-Seq 的 [256,64,32,16];[256,128,64] 对于 SMAGE-seq。

--decodeLayer1：ADT 的低级编码器层：默认：CITE-Seq 的 [16,64,256]。[64,128,256] 为 SMAGE-seq。
--decodeLayer2：高级编码器的层。默认值：[16,20] 表示 CITE-Seq。[64,128,256] 为 SMAGE-seq。
--sigma1：RNA数据上的噪声。默认值：2.5。
--sigma2：ADT数据上的噪声。默认值：CITE-Seq 为 1.5;2.5 对于 SMAGE-Seq
--filter1：如果对基因进行特征选择。默认值：否。
--filter2：如果在ATAC上进行功能选择。默认值：否。
--f1：如果进行羽毛选择，则用于聚类的高可变基因（在X1中）的数量。默认值：2000 -
-f2：如果进行羽毛选择，则来自 ATAC（在 X2 中）的高变量基因数用于聚类。默认值：2000
*为方便起见，我们将 10X 单细胞多组 ATAC + 基因表达技术表示为 SMAGE-seq。

4. run_scMDC

利用scMultiCluster模型联合利用scRNA-seq和其他模态数据进行细胞类型的无监督聚类,并评估了聚类效果。

设置超参数

读取数据，数据包括两个模态:scRNA-seq读数矩阵X1和另一种基因组数据X2(如ADT/ATAC),以及标签y。

    import argparse

    parser = argparse.ArgumentParser(description='train',
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('--n_clusters', default=12, type=int)
    parser.add_argument('--cutoff', default=0.5, type=float,
                        help='Start to train combined layer after what ratio of epoch')
    parser.add_argument('--batch_size', default=256, type=int)
    # parser.add_argument('--data_file', default=f'../datasets/SMAGESeq_10X_pbmc_10k_granulocyte_plus.h5')
    parser.add_argument('--data_file', default=f'../datasets/output500.h5')
    parser.add_argument('--maxiter', default=5000, type=int)
    parser.add_argument('--pretrain_epochs', default=200, type=int)
    parser.add_argument('--gamma', default=.1, type=float,
                        help='coefficient of clustering loss')
    parser.add_argument('--tau', default=.1, type=float,
                        help='fuzziness of clustering loss')
    parser.add_argument('--phi1', default=0.005, type=float,
                        help='coefficient of KL loss in pretraining stage')
    parser.add_argument('--phi2', default=0.005, type=float,
                        help='coefficient of KL loss in clustering stage')
    parser.add_argument('--update_interval', default=1, type=int)
    parser.add_argument('--tol', default=0.001, type=float)
    parser.add_argument('--lr', default=1., type=float)
    parser.add_argument('--ae_weights', default=None)
    parser.add_argument('--save_dir', default='atac_pbmc10k')
    parser.add_argument('--ae_weight_file', default='AE_weights_pbmc10k.pth.tar')
    parser.add_argument('--resolution', default=0.2, type=float)
    parser.add_argument('--n_neighbors', default=30, type=int)
    parser.add_argument('--embedding_file', default=True)
    parser.add_argument('--prediction_file', default=True)
    parser.add_argument('-el', '--encodeLayer', nargs='+', default=[256, 128, 64])
    parser.add_argument('-dl1', '--decodeLayer1', nargs='+', default=[64, 128, 256])
    parser.add_argument('-dl2', '--decodeLayer2', nargs='+', default=[64, 128, 256])
    parser.add_argument('--sigma1', default=2.5, type=float)
    parser.add_argument('--sigma2', default=2.5, type=float)
    parser.add_argument('--f1', default=1000, type=float, help='Number of mRNA after feature selection')
    parser.add_argument('--f2', default=1000, type=float, help='Number of ADT/ATAC after feature selection')
    parser.add_argument('--filter1', default=True, help='Do mRNA selection')
    parser.add_argument('--filter2', default=True, help='Do ADT/ATAC selection')
    parser.add_argument('--run', default=1, type=int)
    parser.add_argument('--device', default='cpu')
    args = parser.parse_args()
    print(args)

输出：

Namespace(ae_weight_file='AE_weights_pbmc10k.pth.tar', ae_weights=None, batch_size=256, cutoff=0.5, data_file='../datasets/output500.h5', decodeLayer1=[64, 128, 256], decodeLayer2=[64, 128, 256], device='cpu', embedding_file=True, encodeLayer=[256, 128, 64], f1=1000, f2=1000, filter1=True, filter2=True, gamma=0.1, lr=1.0, maxiter=5000, n_clusters=12, n_neighbors=30, phi1=0.005, phi2=0.005, prediction_file=True, pretrain_epochs=200, resolution=0.2, run=1, save_dir='atac_pbmc10k', sigma1=2.5, sigma2=2.5, tau=0.1, tol=0.001, update_interval=1)

对scRNA-seq数据进行预处理

规范化,批量效应校正,log变换等。

    data_mat = h5py.File(args.data_file)
    x1 = np.array(data_mat['X1'])
    x2 = np.array(data_mat['X2'])
    y = np.array(data_mat['Y'])
    data_mat.close()

    # #Gene filter
    # if args.filter1:
    #     importantGenes = geneSelection(x1, n=args.f1, plot=False)
    #     x1 = x1[:, importantGenes]
    # if args.filter2:
    #     importantGenes = geneSelection(x2, n=args.f2, plot=False)
    #     x2 = x2[:, importantGenes]

    print('------------------------------------')
    print(x1.shape)
    print(x2.shape)

    # preprocessing scRNA-seq read counts matrix
    adata1 = sc.AnnData(x1)
    adata1.obs['Group'] = y

    adata1 = read_dataset(adata1,
                          transpose=False,
                          test_split=False,
                          copy=True)

    adata1 = normalize(adata1,
                       size_factors=True,
                       normalize_input=True,
                       filter_min_counts=False,
                       logtrans_input=True,
                       nor=0)

    adata2 = sc.AnnData(x2)
    adata2.obs['Group'] = y
    adata2 = read_dataset(adata2,
                          transpose=False,
                          test_split=False,
                          copy=True)

    adata2 = normalize(adata2,
                       size_factors=True,
                       normalize_input=True,
                       logtrans_input=True)

    # adata2 = clr_normalize_each_cell(adata2)

    # adata1.write_h5ad('adata1.h5ad')  # 保存adata1为h5ad文件
    # adata2.write_h5ad('adata2.h5ad')  # 保存adata2为h5ad文件
    # print('------------------------------------')
    # print(adata1.X.shape)
    # print(adata2.X.shape)

    input_size1 = adata1.n_vars
    input_size2 = adata2.n_vars

    print(args)

输出：

(11020, 500)
(11020, 500)
Namespace(ae_weight_file='AE_weights_pbmc10k.pth.tar', ae_weights=None, batch_size=256, cutoff=0.5, data_file='../datasets/output500.h5', decodeLayer1=[64, 128, 256], decodeLayer2=[64, 128, 256], device='cpu', embedding_file=True, encodeLayer=[256, 128, 64], f1=1000, f2=1000, filter1=True, filter2=True, gamma=0.1, lr=1.0, maxiter=5000, n_clusters=12, n_neighbors=30, phi1=0.005, phi2=0.005, prediction_file=True, pretrain_epochs=200, resolution=0.2, run=1, save_dir='atac_pbmc10k', sigma1=2.5, sigma2=2.5, tau=0.1, tol=0.001, update_interval=1)