scDeepCluster_pytorch代码解读与文章理解

OFF JUMPOL

已于 2024-04-14 12:03:26 修改

阅读量347

点赞数

分类专栏： scRNA-seq 文章标签： pytorch 人工智能 python

于 2023-07-04 21:05:04 首次发布

本文链接：https://blog.csdn.net/qq_34179307/article/details/131269530

版权

scRNA-seq 专栏收录该内容

5 篇文章 0 订阅

订阅专栏

论文：Clustering single-cell RNA-seq data with a
model-based deep learning approach
代码：https://github.com/ttgump/scDeepCluster

环境要求：

pytorch 1.8.0
python 3.9
本电脑conda 名称为python39
使用时需要激活环境： conda activate python39

参数

——n_clusters: 聚类数量，如果设置为0，则通过Louvain算法对预训练后的潜在特征进行估计。如果设置为整数> 0，则模型将使用用户定义的值作为簇数。
——knn: Louvain算法中使用的最近邻数，默认值= 20。当设置n_clusters > 0时不使用
——resolution: Louvain算法的分辨率，默认值= 0.8。较大的值将导致更多的集群数。当设置n_clusters > 0时不使用。
——select_genes: 用于分析的选定基因的数量，default = 0。它将使用均值-方差关系来选择信息丰富的基因。建议选择前2000个基因，但这取决于不同的数据集。
——batch_size: 批处理大小，默认= 256。
——data_file: 数据的文件名。
——maxiter: 集群阶段的最大迭代数，默认= 2000。
——pretrain_epochs: 预训练迭代，默认= 300。
——gamma: 聚类损失系数，default = 1。
——sigma: 随机高斯噪声的系数，默认值为2.5。
——update_interval: 更新集群目标的迭代次数，默认= 1。
——tol: 终止聚类阶段的容忍度，即两个连续迭代之间预测标签的增量，default = 0.001。
——final_latent_file: 输出自动编码器的最终潜在表示的文件名，default = final_latent_file.txt。
——predict_label_file: 用于输出集群标签的文件名，default = pred_labels.txt。Louvain算法中使用的最近邻数，默认值= 20。当设置n_clusters > 0时不使用

输出

final_latent: scRNA-seq数据的低维表示，默认形状(n_cells, 32)，可以通过t-SNE或UMAP可视化。
Predict_label: 预测的聚类标签，形状(n_cells)。

运行

对于单细胞计数数据:
python run_scDeepCluster.py --data_file data.h5 --n_clusters 0
将data_file设置为目标数据(以h5格式存储)，具有两个组件X和Y，其中X是细胞按基因计数矩阵，Y是真实标签。Y是可选的)，n_clusters为聚类的数量(0表示通过Louvain算法对预训练的潜在特征进行自动估计)。

代码

run_scDeepCluster.py

获取数据，将数据转化为python中存储单细胞数据的格式–AnnData

    data_mat = h5py.File(args.data_file, 'r')
    x = np.array(data_mat['X'])
    # y is the ground truth labels for evaluating clustering performance
    # If not existing, we skip calculating the clustering performance metrics (e.g. NMI ARI)
    if 'Y' in data_mat:
        y = np.array(data_mat['Y'])
    else:
        y = None
    data_mat.close()

    if args.select_genes > 0:
        importantGenes = geneSelection(x, n=args.select_genes, plot=False)
        x = x[:, importantGenes]

    # preprocessing scRNA-seq read counts matrix
    adata = sc.AnnData(x)

划分数据集

	adata = adata.copy()
    adata.obs['DCA_split'] = 'train'
    adata.obs['DCA_split'] = adata.obs['DCA_split'].astype('category')

adata.X，行对应观测（即，细胞），列对应特征（即，基因）；
**加粗样式**
对数据归一化

# 保留至少具有min_counts计数的基因、细胞
sc.pp.filter_genes(adata, min_counts=1)
sc.pp.filter_cells(adata, min_counts=1)
adata.raw = adata.copy()
# 通过所有基因的总计数对每个细胞进行归一化，这样每个细胞在归一化后都有相同的总计数。
sc.pp.normalize_per_cell(adata)
adata.obs['size_factors'] = adata.obs.n_counts / np.median(adata.obs.n_counts)
# Computes X = log(X+1)
sc.pp.log1p(adata)
# 将数据缩放为单位方差和零均值。
sc.pp.scale(adata)

# scanpy.pp.filter_genes 根据细胞数量或计数过滤基因。
scanpy.pp.filter_genes(data, min_counts=None, min_cells=None, max_counts=None, max_cells=None, inplace=True, copy=False)

预处理步骤和DCA一样：
首先，在任何细胞中没有计数的基因都被过滤掉。
然后，计算大小因子（size factor），并规范化读取计数，因此每个细胞的count number的总数是相同的。形式上，如果我们将细胞 $i$ 的总读取计数表示为 $l_i$ ，则细胞 $i$ 的大小因子为 $l_i/median(l)$ 。
最后一步是对读取计数进行对数变换和scale，使数据遵循单位方差和零均值。

模型

model = scDeepCluster(input_dim=adata.n_vars, z_dim=32, 
                encodeLayer=[256, 64], decodeLayer=[64, 256], sigma=args.sigma, gamma=args.gamma, device=args.device)
                
model.pretrain_autoencoder(X=adata.X, X_raw=adata.raw.X, size_factor=adata.obs.size_factors, 
                                batch_size=args.batch_size, epochs=args.pretrain_epochs, ae_weights=args.ae_weight_file)

y_pred, _, _, _, _ = model.fit(X=adata.X, X_raw=adata.raw.X, size_factor=adata.obs.size_factors, n_clusters=args.n_clusters, init_centroid=None, 
                    y_pred_init=None, y=y, batch_size=args.batch_size, num_epochs=args.maxiter, update_interval=args.update_interval, tol=args.tol, save_dir=args.save_dir)

scDeepCluster.py

model.pretrain_autoencoder
在这里插入图片描述

model.fit
输入 X.raw(没有归一化，不带噪声) 到Encoder-ec_mu 得到数据的潜在特征编码（潜在空间的数据特征维度32 远远小于原始数据 16653），用K-means对潜在的特征聚类，得到初始聚类中心

输入 X（归一化，不带噪声)到Encoder-ec_mu 得到数据的潜在特征编码 latent z,计算软分类q来度量z和聚类中心之间的相似度， $q_{ij}$ 的计算其实就是 softmax，所以也就表示在当前的 embedding 情况下，第i个数据属于第j个簇的概率。
在这里插入图片描述
根据 q计算目标分布ｐ
Ｐ是为了满足加强预测结果，更多强调有高置信度的点，正常化损失函数的要求被作者设计出来的。
$f_j=\sum q_{ij}^2$ 就表示第 j 类的大小。因此设计这个的意思就是 $q_{ij}$ 先平方，再 normalize （对类的规模惩罚），再 softmax 的概率。

那么为什么要平方呢？我觉得是为了满足作者提出分布需要满足的前两个性质：平方之后，让大的更大、小的更小，满足了 1；同时，如果样本 i 属于 j 类的概率比其他样本属于 j 类的概率高得多，那么算出来的也会大得多，这样就满足了 2。因此在最终的优化过程中，”DEC通过从高置信度预测中学习来改进每次迭代中的初始估计，这反过来又有助于改进低置信度预测”。

计算两次迭代之间ｙ的预测不同占比为多少，如果变化不大（＜0.001）则停止训练。
ｑ是软分类，q argmax后为ｙ的预测结果。
第一个Epoch用Kmeans得到的分类结果和当前ｑ预测的结果判断停止条件。

按照Batch的大小更新 z 和对应的q, 以及 $\mu, \theta, \pi$ ，再次运用AE模型，
计算当前epoch的p和此batch的KL散度作为cluster loss,计算当前batch的 x.raw和 $\mu, \theta, \pi$ 生成的分布模型的zinb loss作为重建loss,