occdn.limour.top/source/_posts/单细胞水平的细胞周期评分.md

---

title: 单细胞水平的细胞周期评分 tags: [] id: '1573' categories:

• - 未分类
• - 生信 date: 2022-02-26 13:04:36 ---

定义一些辅助函数

library(ggsci)
f_levels2col <- function(df_f, col){
    col[as.numeric(df_f)]
}

f_getdfcol <- function(df, cN){
    res <- df[[cN]]
    names(res) <- rownames(df)
    res
}
f_metaG2G <- function(metaG, matrixN=F){
    res  <- list()
    alltype <- unique(metaG[[1]])
    for(type in alltype){
        res[[type]] <- rownames(metaG)[metaG[[1]] == type]
        if (matrixN){
            res[[type]] <- gsub('-','.',res[[type]])
        }
    }
    res
}
 
f_br_cluster_f <- function(sObject, lc_groupN, is_sce){
    if(is_sce){
        lc_filter <- unlist(unique(sObject[[lc_groupN]]))
    }else{
        lc_filter <- unlist(unique(sObject[lc_groupN]))
    }
    lc_filter <- lc_filter[!is.na(lc_filter)]
    lc_filter
}
 
f_br_cluster <- function(sObject, lc_groupN, lc_labelN, lc_prop = F, is_sce=T){
    if(is_sce){
        lc_g <- f_metaG2G(sObject[[lc_groupN]])
        lc_l <- sObject[[lc_labelN]]
    }else{
        lc_g <- f_metaG2G(sObject[lc_groupN])
        lc_l <- sObject[lc_labelN]
    }
    lc_l[[1]] <- as.character(lc_l[[1]])
    res <- data.frame(row.names = f_br_cluster_f(lc_l, lc_labelN, is_sce))
    if(lc_prop){
        for(Nm in names(lc_g)){
            tmp <- prop.table(table(lc_l[lc_g[[Nm]],]))
            res[[Nm]] <- tmp[rownames(res)]
        }
    }else{
        for(Nm in names(lc_g)){
            tmp <- table(lc_l[lc_g[[Nm]],])
            res[[Nm]] <- tmp[rownames(res)]
        }
    }
    res[is.na(res)] = 0
    res
}

数据准备

library(Seurat)
sce <- CellCycleScoring(sce, s.features = cc.genes$s.genes, g2m.features = cc.genes$g2m.genes, set.ident = TRUE)
df <- sce@meta.data
df <- subset(df, immune_annotation != 'immune')
df$Phase_status <- df$Phase
idx <- with(df, S.Score < (median(S.Score) + 4 * mad(S.Score)) & G2M.Score < (median(G2M.Score) + 4 * mad(G2M.Score)))
df[idx, 'Phase_status'] = 'Low cycling'
df[!idx, 'Phase_status'] = 'High cycling'

火山图绘制

p <- ggplot() + geom_point(aes(S.Score, G2M.Score, col=patient_id, shape=cell_type, size=Phase_status, alpha=Phase_status), data = df, alpha=f_levels2col(df$Phase_status, c(0.1, 0.5)))
p <- p + theme_bw()
p
ggsave(p, filename = 'fig1.E_12inch.pdf', width = 12, height = 12)

细胞比例检验

f_test_percent <- function(dfA, dfB, colN, rowN){
    n1 <- sum(dfA[[colN]])
    n2 <- sum(dfB[[colN]])
    a <- dfA[rowN, colN]
    c <- dfB[rowN, colN]
    b <- n1 - a 
    d <- n2 - c
    mat <- matrix(c(a,c,b,d), ncol = 2, nrow = 2)
    res <- chisq.test(mat)
    if(min(res$expected) < 5){
        res <- fisher.test(mat)
    }
    res
}

mat4_C <- f_br_cluster(subset(df, group=='CRPC'), 'cell_type', 'Phase_status', is_sce=F)
mat4_H <- f_br_cluster(subset(df, group=='HSPC'), 'cell_type', 'Phase_status', is_sce=F)

f_test_percent(mat4_C, mat4_H, 'Luminal', 'High cycling')
f_test_percent(mat4_C, mat4_H, 'Fibroblasts', 'High cycling')
f_test_percent(mat4_C, mat4_H, 'Endothelial', 'High cycling')
f_test_percent(mat4_C, mat4_H, 'Basal cell', 'High cycling')
a <- t(mat4_H)
t(a/rowSums(a))
a <- t(mat4_C)
t(a/rowSums(a))