聚类迭代

[xiangpin]

根据上面两段话的描述，将这个过程用代码实现。

主要过程：

1. 根据bulk数据的肿瘤与对照鉴定差异基因，获取肿瘤（恶性）与对照（非恶性）的高表达基因。
2. 根据上面获得的基因，在单细胞数据集中的上皮细胞（或其它类型）进行打分（恶性得分与非恶性得分），随后用K-mean进行聚类
   得到伪恶性与非恶性两类。
3. 再根据伪恶性与非恶性两类的单细胞上皮细胞类群，然后重复1,2步骤，即进行差异分析，获得恶性与非恶性的高表达基因，然后
   再打分，分类，直到上一轮的上皮细胞类群标签（恶性与非恶性）与本轮的上皮细胞（恶性与非恶性）标签数基本一致（比如>=90%）
   再停止迭代。

[GuangchuangYu]

第一步写一个函数，第二步和第三步写成一个函数。

参考文献为: doi:10.1136/gutjnl-2019-320368

[P10911004-NPUST]

你好，這裡提供一個可能的解決方案。

rm(list = ls())
suppressPackageStartupMessages({
    if ( ! require(tidyverse) ) install.packages("tidyverse")
    if ( ! require(BiocManager) ) install.packages("BiocManager")
    if ( ! require(Seurat) ) BiocManager::install("Seurat", ask = FALSE)
    if ( ! require(edgeR) ) BiocManager::install("edgeR", ask = FALSE)
})

data("pbmc_small") # A small dataset for demonstration


init_cluster <- function(
        obj,  # An `Seurat` object
        top_n = 100,  # How many top DEGs
        ctrl = 5,  # How many genes will be sampled from each bin
        centers = 2  # How many clusters we assume
){
    # 第一个 loop 先生成初始标签
    
    # 取得 Differentially Expressed Genes (DEGs)
    obj <- Seurat::FindVariableFeatures(obj, verbose = FALSE)
    top_genes <- head(Seurat::VariableFeatures(obj), top_n)
    
    # 根据提供的DEGs打分数
    score <- Seurat::AddModuleScore(obj, features = list(top_genes), ctrl = ctrl)
    
    # Seurat::AddModuleScore()生成的分数会自动插入meta.data这个slot的最后一栏，
    # 其为dataframe，因此仍需[[1]]以转换成atomic vector
    score <- score[[ncol(score@meta.data)]][[1]]
    
    # kmeans
    km <- stats::kmeans(score, centers = centers)
    ret <- list(
        "score" = score,
        "cluster" = km$cluster
    )
    
    return(ret)
}


get_DEGs <- function(
        obj,  # An `Seurat` object
        cluster,  # The cluster labels obtain from kmeans clustering
        top_n = 100  # How many top DEGs
){
    # 我不熟悉scRNA-seq，这个函数的DEGs获取方式大概率不恰当。
    # 因此，把这个部分独立出来成一个函数，方便更换。
    y <- edgeR::DGEList(counts = obj[["RNA"]]$data, group = cluster)
    group <- factor(cluster)
    design <- model.matrix( ~ group)
    fit <- edgeR::glmQLFit(y, design)
    qlf <- edgeR::glmQLFTest(fit, coef = 2)
    DEGs <- topTags(qlf, n = top_n)
    return(rownames(DEGs))
}


get_score <- function(
        obj,  # An `Seurat` object
        DEGs,  # Differentially expressed genes (a character vector)
        ctrl = 5  # How many genes will be sampled from each bin
){
    if ( ! inherits(DEGs, "list") ) DEGs <- list(DEGs)
    score <- Seurat::AddModuleScore(obj, features = DEGs, ctrl = ctrl)
    score <- score[[ncol(score@meta.data)]][[1]]
    return(score)
}


#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Main function -- Iterative kmeans
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
kmeans_iter <- function(
        obj,  # An `Seurat` object
        top_n = 100,  # How many top DEGs
        ctrl = 5,  # How many genes will be sampled from each bin
        centers = 2,  # How many clusters we assume
        accuracy = 0.7,  # The iteration stop when accuracy reached this value
        max_iter = 200  # The iteration stop when iteration exceed this value
){
    i <- 0 # 迭代次数
    j <- 0 # accuracy 达标之计数
    DEGs_vct <- character(0) # 記錄多次出現的基因
    
    while (i <= max_iter)
    {
        i <- i + 1
        
        if (i == 1) {
            # 第一个 loop 先生成初始标签
            ret <- init_cluster(
                obj = obj, 
                top_n = top_n, 
                ctrl = ctrl, 
                centers = centers
            )
            ret[["accuracy"]] <- numeric(0)
            next
        }
        
        clust1 <- ret$cluster # 上一轮生成的分群标签
        DEGs <- get_DEGs(obj, cluster = clust1, top_n = top_n)
        score <- get_score(obj, DEGs, ctrl)
        
        km <- kmeans(score, centers)
        clust2 <- km$cluster # 这一轮得到的分群标签
        
        ret$score <- score
        ret$cluster <- clust2
        
        cfm <- table(clust1, clust2) # 上一輪與這一輪的標籤匹配情況
        acc <- sum(diag(cfm)) / sum(cfm) # 有多少標籤和上一輪是一致的
        ret[["accuracy"]] <- append(ret[["accuracy"]], acc)
        
        # 當連續三次準確率大於目標閾值 (accuracy)，則視為模型已穩定，結束迭代
        # 这三次参与的 DEGs 也记录下来，按出现次数倒序排列，出现越多次，可能越重要
        if (acc < accuracy)
        {
            j <- 0
            DEGs_vct <- character(0)
        } else {
            j <- j + 1
            DEGs_vct <- append(DEGs_vct, DEGs)
        }
        # j <- ifelse(acc < accuracy, 0, j + 1)
        if (j == 3) {
            cat("=====================\n")
            cat("Done!\n\n")
            cat(sprintf("Accuracy: %s\nIterations: %s/%s\n", acc, i, max_iter))
            cat("=====================\n")
            break
        }
    }
    
    ret[["DEGs"]] <- sort(table(DEGs_vct), decreasing = TRUE)[1:top_n]
    
    return(ret)
}

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
# Testing
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
out <- kmeans_iter(pbmc_small)
out