「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析

本文主要是介绍「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

这一篇文章是回答知识星球中一位星友的提问,她的电脑内存有限,无法直接使用所有数据,只能分析部分数据。

数据来源: https://content.cruk.cam.ac.uk/jmlab/atlas_data.tar.gz 解压缩之后,得到下面数据

数据清单

其中raw_counts.mtx是以稀疏矩阵格式存放的表达量数据,文件为6.5G, 用普通的文本编辑器无法打开,我们可以用Linux命令行的less查看数据存放形式

稀疏矩阵

显然这种格式并不是给人类阅读的,它存放的是非零数据的位置及其具体数值。当然,我们也不需要读懂,只需要R语言或者其他编程语言能够加载即可。

R语言的Matrix包的readMM函数就能够读取该文件

mt <- Matrix::readMM("raw_counts.mtx")
dim(mt)
# [1]  29452 139331
# 行为基因,列为细胞

这一步时间非常的久,我差不多花了10分钟时间。同时占用内存也非常可观,直接占用了8G左右的内存,不到16G内存的电脑可能根本无法读取。

format(object.size(mt), units = "Mb")
# "7377.8 Mb"

稀疏矩阵其实和普通矩阵看起来差不多,除了在显示的时候用.来表示0.

查看部分信息

还有一点就是,对于这种量级的数据,我们无法使用R自带的as.data.frame或者as.matrix将其转成普通的数据库或者矩阵,它会直接报错。因此我也不建议对其进行数据转换。

我们发现这里的矩阵并没有行名和列名,这部分信息需要额外从其他文件中读取

bc <- read.table("barcodes.tsv")
genes <- read.table("genes.tsv", sep = "\t")
dim(bc)
#[1] 139331      1
dim(genes)
#[1] 29452     2

不难发现barcode的行数等于矩阵的列数, gene的行数等于矩阵的行数, 也就是说矩阵的列是细胞,行是基因。

row.names(mt) <- genes$V1
colnames(mt) <- bc$V1
查看部分信息

建议:将此处得到matrix保存为Rds格式,方便后续加载

saveRDS(mt, "raw_matrix.Rds")

接下来就是根据元信息来提取对应的细胞,我们以提取"Mesenchyme"细胞为例进行讲解

meta.info <- read.table("meta.tab",sep = "\t", header = TRUE)cell.info <- meta.info[meta.info$celltype == "Mesenchyme", "cell"]
cell.info <- cell.info[!is.na(cell.info)]mt.sml <- mt[, cell.info]
format(object.size(mt.sml), units = "Mb")
# "280.9 Mb"

代码的核心逻辑为提取出对应行的细胞名,然后根据细胞名提取矩阵中的对应列。

过滤后的细胞就可以用作后续分析。不过在开始分析之前,让我们先把原始的矩阵给删掉,因为它实在是太占用内存了。

rm(mt); gc()

除了用元信息进行过滤外,你还可以通过随机抽样,从原始数据中抽出部分细胞,这样子也能够在内存吃紧的情况进行后续分析。

这篇关于「单细胞转录组系列」如何从稀疏矩阵中提取部分数据进行分析的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1088132

相关文章

Spring Security 从入门到进阶系列教程

Spring Security 入门系列 《保护 Web 应用的安全》 《Spring-Security-入门(一):登录与退出》 《Spring-Security-入门(二):基于数据库验证》 《Spring-Security-入门(三):密码加密》 《Spring-Security-入门(四):自定义-Filter》 《Spring-Security-入门(五):在 Sprin

大模型研发全揭秘:客服工单数据标注的完整攻略

在人工智能(AI)领域,数据标注是模型训练过程中至关重要的一步。无论你是新手还是有经验的从业者,掌握数据标注的技术细节和常见问题的解决方案都能为你的AI项目增添不少价值。在电信运营商的客服系统中,工单数据是客户问题和解决方案的重要记录。通过对这些工单数据进行有效标注,不仅能够帮助提升客服自动化系统的智能化水平,还能优化客户服务流程,提高客户满意度。本文将详细介绍如何在电信运营商客服工单的背景下进行

基于MySQL Binlog的Elasticsearch数据同步实践

一、为什么要做 随着马蜂窝的逐渐发展,我们的业务数据越来越多,单纯使用 MySQL 已经不能满足我们的数据查询需求,例如对于商品、订单等数据的多维度检索。 使用 Elasticsearch 存储业务数据可以很好的解决我们业务中的搜索需求。而数据进行异构存储后,随之而来的就是数据同步的问题。 二、现有方法及问题 对于数据同步,我们目前的解决方案是建立数据中间表。把需要检索的业务数据,统一放到一张M

关于数据埋点,你需要了解这些基本知识

产品汪每天都在和数据打交道,你知道数据来自哪里吗? 移动app端内的用户行为数据大多来自埋点,了解一些埋点知识,能和数据分析师、技术侃大山,参与到前期的数据采集,更重要是让最终的埋点数据能为我所用,否则可怜巴巴等上几个月是常有的事。   埋点类型 根据埋点方式,可以区分为: 手动埋点半自动埋点全自动埋点 秉承“任何事物都有两面性”的道理:自动程度高的,能解决通用统计,便于统一化管理,但个性化定

使用SecondaryNameNode恢复NameNode的数据

1)需求: NameNode进程挂了并且存储的数据也丢失了,如何恢复NameNode 此种方式恢复的数据可能存在小部分数据的丢失。 2)故障模拟 (1)kill -9 NameNode进程 [lytfly@hadoop102 current]$ kill -9 19886 (2)删除NameNode存储的数据(/opt/module/hadoop-3.1.4/data/tmp/dfs/na

异构存储(冷热数据分离)

异构存储主要解决不同的数据,存储在不同类型的硬盘中,达到最佳性能的问题。 异构存储Shell操作 (1)查看当前有哪些存储策略可以用 [lytfly@hadoop102 hadoop-3.1.4]$ hdfs storagepolicies -listPolicies (2)为指定路径(数据存储目录)设置指定的存储策略 hdfs storagepolicies -setStoragePo

Hadoop集群数据均衡之磁盘间数据均衡

生产环境,由于硬盘空间不足,往往需要增加一块硬盘。刚加载的硬盘没有数据时,可以执行磁盘数据均衡命令。(Hadoop3.x新特性) plan后面带的节点的名字必须是已经存在的,并且是需要均衡的节点。 如果节点不存在,会报如下错误: 如果节点只有一个硬盘的话,不会创建均衡计划: (1)生成均衡计划 hdfs diskbalancer -plan hadoop102 (2)执行均衡计划 hd

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题,本文将一起从代码层去分析为什么没有建立索引? 开源ERP项目地址:https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi

科研绘图系列:R语言扩展物种堆积图(Extended Stacked Barplot)

介绍 R语言的扩展物种堆积图是一种数据可视化工具,它不仅展示了物种的堆积结果,还整合了不同样本分组之间的差异性分析结果。这种图形表示方法能够直观地比较不同物种在各个分组中的显著性差异,为研究者提供了一种有效的数据解读方式。 加载R包 knitr::opts_chunk$set(warning = F, message = F)library(tidyverse)library(phyl