Hard 随机选择subset @CareerCup

本博客介绍了一个函数,用于随机从大小为n的数组中选取m个整数,确保每个元素被选中的概率相等。通过递归和交换实现随机选择,并在最后输出选择的整数集合。

算法同上题


package Hard;

import CtCILibrary.AssortedMethods;


/**
 * Write a method to randomly generate a set of m integers from an array of size n. Each element must have equal probability of being chosen.

译文:

写一个函数,随机地从大小为n的数组中选取m个整数。要求每个元素被选中的概率相等。
 *
 */
public class S18_3 {

	/* Random number between lower and higher, inclusive */
    public static int rand(int lower, int higher) { 
            return lower + (int)(Math.random() * (higher - lower + 1));
    }

    public static void swap(int[] a, int n, int m){
    	int tmp = a[n];
    	a[n] = a[m];
    	a[m] = tmp;
    }
    
    
    /* pick M elements from original array, using only elements 0 through i (inclusive).*/
    public static int[] pickMRecursively(int[] original, int m, int i) {
            if (i + 1 < m) { // Not enough elements
                return null; 
            } else if (i + 1 == m) { // Base case -- copy first m elements into array
                int[] set = new int[m];
                for (int k = 0; k < m; k++) {
                    set[k] = original[k];
                }
                return set;
            } else {
                int[] set = pickMRecursively(original, m, i - 1);
                int k = rand(0, i);
                if (k < m) {
                    set[k] = original[i];
                }
                return set;
            }
    }
    
    /* pick M elements from original array.  Clone original array so that
     * we don’t destroy the input. */
    public static int[] pickMRandomly(int[] original, int m) {
        for (int i = 0; i < m; i++) { 
        	int k = rand(i, original.length-1);		// 产生i到n-1间的随机数
            swap(original, i, k);
        }
        int[] subset = new int[m];
        for(int i=0; i<m; i++){
        	subset[i] = original[i];
        }
        return subset;
    }
    
    public static void main(String[] args) {
            int[] cards = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
            System.out.println(AssortedMethods.arrayToString(cards));
            int[] set = pickMRandomly(cards, 4);
            System.out.println(AssortedMethods.arrayToString(set));
    }
}


<think>这个错误信息表明:在数据框(data.frame)子集选择时,指定的列名(score_cols)在数据框中不存在。具体来说,是在`gp3_subset@meta.data`中试图选择`score_cols`所指定的列,但其中有些列名并不存在。可能的原因:1.`AddModuleScore`函数没有成功添加模块分数,因此元数据中没有生成以"CT_"为前缀的列。2.`score_cols`的生成方式有误,导致列名与实际添加的列名不匹配。解决步骤:1.首先检查`AddModuleScore`函数执行后,元数据中是否添加了新的列。可以通过以下命令查看:```rcolnames(gp3_subset@meta.data)```2.确认`AddModuleScore`函数生成的列名。根据Seurat的`AddModuleScore`函数,它会生成一组列,列名由`name`参数和从1开始的数字组成(注意:在R中,索引通常从1开始,但这里生成的是从0开始?实际上,在代码中我们使用了`seq_along(markers_list)-1`,所以生成的是"CT_0","CT_1",...)。然而,根据Seurat文档,`AddModuleScore`生成的列名是`name`参数加上数字(从1开始,而不是0)。但是,我们注意到在代码中指定了`name="CT_"`,并且我们在生成`score_cols`时使用了`seq_along(markers_list)-1`(即从0开始)。这可能是问题的根源。让我们查看`AddModuleScore`的文档:实际上,它生成的列名是`name`后面接一个数字,这个数字从1开始,到标记基因集的数量。例如,如果有3个标记基因集,那么列名就是`CT_1`,`CT_2`,`CT_3`。因此,我们应该生成列名为`paste0("CT_",1:length(markers_list))`。修正代码:将score_cols<-paste0("CT_",seq_along(markers_list)-1)改为score_cols<-paste0("CT_",1:length(markers_list))但是,注意:`seq_along(markers_list)`返回的是1,2,...,n,所以`seq_along(markers_list)-1`返回0,1,...,n-1。而`AddModuleScore`生成的是从1开始的数字。因此,我们只需要使用1到n。另外,还有一种可能是`AddModuleScore`在最新版本中改变了列名的生成方式。因此,最可靠的方式是查看函数执行后元数据中的列名。修正后的代码片段:```rif(length(markers_list)>0){gp3_subset<-AddModuleScore(object=gp3_subset,features=markers_list,name="CT_",#更清晰的列名前缀assay="RNA",nbin=12#减少bin数量提高计算效率)#6.分配细胞类型#注意:AddModuleScore生成的列名是"CT_1","CT_2",...直到"CT_n"score_cols<-paste0("CT_",1:length(markers_list))#确保这些列存在if(!all(score_cols%in%colnames(gp3_subset@meta.data))){stop("Somemodulescorecolumnsaremissing.CheckAddModuleScoreoutput.")}gp3_subset$cell_type<-names(markers_list)[apply(gp3_subset@meta.data[,score_cols,drop=FALSE],1,which.max)]#7.分配大类gp3_subset$broad_category<-"Unassigned"for(categoryinnames(cell_categories)){gp3_subset$broad_category[gp3_subset$cell_type%in%cell_categories[[category]]]<-category}#检查未分配的细胞unassigned<-sum(gp3_subset$broad_category=="Unassigned")cat("Unassignedcells:",unassigned,"(",round(unassigned/ncol(gp3_subset)*100,2),"%)\n")}else{stop("Novalidmarkersetsfound.Checkmarkergenefile.")}```此外,为了增强代码的健壮性,我们在使用`score_cols`之前检查了这些列是否存在。如果问题仍然存在,请检查`AddModuleScore`是否成功运行,以及是否因为某种原因没有添加列(例如,标记基因集为空?)。另外,注意`AddModuleScore`函数在计算时可能会因为基因不存在而报错,但这里没有捕获错误。所以,最好在运行前确保每个标记基因集中至少有一些基因存在于数据集中。如果标记基因集中没有基因存在于数据集中,那么该基因集将不会产生有效的分数列。因此,在创建`markers_list`时,应该先过滤掉数据集中不存在的基因。我们可以添加一个步骤:在调用`AddModuleScore`之前,过滤每个标记基因集,只保留在`gp3_subset`中存在的基因。但注意:`AddModuleScore`内部已经处理了这种情况(它只会使用存在的基因),所以即使标记基因集中有基因不存在,函数也会运行,只是使用存在的基因计算。因此,不会导致列不生成。所以,主要问题还是列名的索引从1开始,而不是0。让我们按照上述修正代码,并添加必要的检查。如果问题解决,就可以继续可视化部分。修正后的完整代码(包括之前的修正)如下:</think>这个错误表明 `score_cols` 中的列名在 `gp3_subset@meta.data` 中不存在。这通常是由 `AddModuleScore` 函数未能成功创建分数列导致的。以下是完整的解决方案: ```r # 确保 Seurat 对象有效 library(Seurat) # 1. 检查 AddModuleScore 是否成功添加了分数列 if (length(markers_list) > 0) { # 执行 AddModuleScore gp3_subset <- AddModuleScore( object = gp3_subset, features = markers_list, name = "CT_", # 列名前缀 assay = "RNA", nbin = 12 ) # 2. 验证分数列是否存在 # 获取所有以"CT_"开头的列名 actual_score_cols <- grep("^CT_", colnames(gp3_subset@meta.data), value = TRUE) if (length(actual_score_cols) == length(markers_list)) { # 3. 分配细胞类型 - 使用实际存在的列名 gp3_subset$cell_type <- names(markers_list)[apply( gp3_subset@meta.data[, actual_score_cols, drop = FALSE], 1, which.max )] # 4. 分配大类 gp3_subset$broad_category <- "Unassigned" for (category in names(cell_categories)) { gp3_subset$broad_category[gp3_subset$cell_type %in% cell_categories[[category]]] <- category } # 5. 检查未分配的细胞 unassigned <- sum(gp3_subset$broad_category == "Unassigned") cat("Unassigned cells:", unassigned, "(", round(unassigned/ncol(gp3_subset)*100, 2), "%)\n") } else { stop(paste("AddModuleScore failed to create expected columns.", "Expected:", length(markers_list), "columns,", "Found:", length(actual_score_cols))) } } else { stop("No valid marker sets found. Check marker gene file.") } # 6. 可视化 col_palette <- c( "Stromal cells" = "#1F77B4", "Immune cells" = "#FF7F0E", "Tumor cells" = "#D62728", "Unassigned" = "grey80" ) # 使用实际存在的分组列进行绘图 if ("broad_category" %in% colnames(gp3_subset@meta.data)) { p1 <- DimPlot(gp3_subset, reduction = "umap_harmony", group.by = "broad_category", cols = col_palette, pt.size = 1.2, order = c("Tumor cells", "Immune cells", "Stromal cells", "Unassigned")) + ggtitle("GP3 Broad Categories") + theme_classic() print(p1) } else { warning("broad_category column not found. Skipping broad category plot.") } if ("cell_type" %in% colnames(gp3_subset@meta.data)) { p2 <- DimPlot(gp3_subset, reduction = "umap_harmony", group.by = "cell_type", label = TRUE, repel = TRUE, pt.size = 1.0) + ggtitle("Detailed Cell Types") + theme_classic() print(p2) } else { warning("cell_type column not found. Skipping detailed cell type plot.") } ``` ### 错误原因和解决方案: 1. **`AddModuleScore` 未创建预期列**: - 原因:标记基因在数据集中不存在或格式不正确 - 解决:使用 `grep("^CT_", colnames(gp3_subset@meta.data), value = TRUE)` 获取实际创建的列名 - 添加验证步骤检查列数是否匹配 2. **标记基因集问题**: ```r # 正确格式示例 - 命名列表的列表 markers_list <- list( "T cells" = c("CD3D", "CD3E", "CD8A"), "B cells" = c("CD19", "MS4A1", "CD79A"), "Macrophages" = c("CD14", "CD68", "AIF1") ) ``` - 确保每个子列表都有名称 - 确认基因名与数据集中的一致(大小写敏感) 3. **数据预处理问题**: - 确保数据已正确归一化:`gp3_subset <- NormalizeData(gp3_subset)` - 确认所有标记基因在数据集中存在: ```r all_genes <- rownames(gp3_subset) for (ct in names(markers_list)) { missing <- setdiff(markers_list[[ct]], all_genes) if (length(missing) > 0) { warning(paste("Missing genes for", ct, ":", paste(missing, collapse = ", "))) } } ``` 4. **备用细胞类型分配方法**: 如果 `AddModuleScore` 仍然失败,可以使用替代方法: ```r # 方法1:使用 Seurat 的 AddModuleScore_UCell if (!require(UCell)) install.packages("UCell") library(UCell) scores <- ScoreSignatures_UCell(gp3_subset@assays$RNA@data, features = markers_list) gp3_subset@meta.data <- cbind(gp3_subset@meta.data, scores) # 方法2:手动计算平均表达 avg_exp <- sapply(markers_list, function(genes) { valid_genes <- intersect(genes, rownames(gp3_subset)) if (length(valid_genes) == 0) return(rep(NA, ncol(gp3_subset))) rowMeans(GetAssayData(gp3_subset, slot = "data")[valid_genes, ]) }) gp3_subset$cell_type <- names(markers_list)[apply(avg_exp, 1, which.max)] ``` ### 完整工作流程建议: 1. 数据预处理:归一化、缩放、PCA 2. 降维和聚类:UMAP/tSNE、FindClusters 3. 标记基因验证:检查基因是否存在于数据中 4. 运行 `AddModuleScore` 并验证输出列 5. 分配细胞类型和大类 6. 可视化并验证结果
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值