Signac R|如何合并多个 Seurat 对象 (1)

引言

在本文中演示了如何合并包含单细胞染色质数据的多个 Seurat 对象。为了进行演示,将使用 10x Genomics 提供的四个 scATAC-seq PBMC 数据集:

  1. 500-cell PBMC

  2. 1k-cell PBMC

  3. 5k-cell PBMC

  4. 10k-cell PBMC

实战

在整合多个单细胞染色质数据集的过程中,应当意识到,如果对每个数据集单独进行了峰值检测,那么检测到的峰值可能并不完全一致。为此,需要构建一个适用于所有合并数据集的共通峰值集合。

可以通过GenomicRanges包提供的功能来实现这一点。该包中的reduce函数能够将所有相互重叠的峰值进行合并。此外,disjoin函数也是一个不错的选择,它能够生成互不重叠的独立峰值集合。以下是一个图解示例,用以展示reducedisjoin两种方法的差异:

gr <- GRanges(seqnames = "chr1", ranges = IRanges(start = c(2070300), end = c(120200400)))
gr

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1    20-120      *
##   [2]     chr1    70-200      *
##   [3]     chr1   300-400      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
alt
alt

构建统一的峰值集合

如果在各个实验中分别鉴定了峰值,那么这些峰值可能并不完全一致。可以通过整合所有实验数据中的峰值,来形成一个统一的峰值集合,并在合并这些数据集之前,对每个实验中的峰值集合进行量化分析。

首先,需要导入每个实验的峰值位置信息,并将其转换成基因组范围内的格式。接着,利用 GenomicRanges 包中的 reduce 函数,来创建一个适用于所有数据集的峰值集合,以便在每个实验中进行量化分析。

library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)

plan("multicore", workers = 4)
options(future.globals.maxSize = 50000 * 1024^2# for 50 Gb RAM

# read in peak sets
peaks.500 <- read.table(
  file = "pbmc500/atac_pbmc_500_nextgem_peaks.bed",
  col.names = c("chr""start""end")
)
peaks.1k <- read.table(
  file = "pbmc1k/atac_pbmc_1k_nextgem_peaks.bed",
  col.names = c("chr""start""end")
)
peaks.5k <- read.table(
  file = "pbmc5k/atac_pbmc_5k_nextgem_peaks.bed",
  col.names = c("chr""start""end")
)
peaks.10k <- read.table(
  file = "pbmc10k/atac_pbmc_10k_nextgem_peaks.bed",
  col.names = c("chr""start""end")
)

# convert to genomic ranges
gr.500 <- makeGRangesFromDataFrame(peaks.500)
gr.1k <- makeGRangesFromDataFrame(peaks.1k)
gr.5k <- makeGRangesFromDataFrame(peaks.5k)
gr.10k <- makeGRangesFromDataFrame(peaks.10k)

# Create a unified set of peaks to quantify in each dataset
combined.peaks <- reduce(x = c(gr.500, gr.1k, gr.5k, gr.10k))

# Filter out bad peaks based on length
peakwidths <- width(combined.peaks)
combined.peaks <- combined.peaks[peakwidths  < 10000 & peakwidths > 20]
combined.peaks

## GRanges object with 89951 ranges and 0 metadata columns:
##           seqnames            ranges strand
##              <Rle>         <IRanges>  <Rle>
##       [1]     chr1     565153-565499      *
##       [2]     chr1     569185-569620      *
##       [3]     chr1     713551-714783      *
##       [4]     chr1     752418-753020      *
##       [5]     chr1     762249-763345      *
##       ...      ...               ...    ...
##   [89947]     chrY 23422151-23422632      *
##   [89948]     chrY 23583994-23584463      *
##   [89949]     chrY 23602466-23602779      *
##   [89950]     chrY 28816593-28817710      *
##   [89951]     chrY 58855911-58856251      *
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

构建片段对象

为了对合并的峰值集合进行量化分析,需要针对每个实验创建一个片段对象。片段对象是一个在 Signac 中特别定义的类,它负责存储与单个片段文件相关的所有数据。

首先需要导入每个实验的细胞元数据,这样就能了解每个文件包含哪些细胞的条形码信息。之后,利用 CreateFragmentObject 函数来生成片段对象。这个函数会进行一系列验证,确保文件不仅存在于硬盘上,而且已经过压缩和索引处理,同时计算文件及其 tabix 索引的 MD5 校验值,以便于能够检测到文件在任何时间点的变更情况,并确认文件中确实包含了预期的细胞类型。

# load metadata
md.500 <- read.table(
  file = "pbmc500/atac_pbmc_500_nextgem_singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ] # remove the first row

md.1k <- read.table(
  file = "pbmc1k/atac_pbmc_1k_nextgem_singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ]

md.5k <- read.table(
  file = "pbmc5k/atac_pbmc_5k_nextgem_singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ]

md.10k <- read.table(
  file = "pbmc10k/atac_pbmc_10k_nextgem_singlecell.csv",
  stringsAsFactors = FALSE,
  sep = ",",
  header = TRUE,
  row.names = 1
)[-1, ]

# perform an initial filtering of low count cells
md.500 <- md.500[md.500$passed_filters > 500, ]
md.1k <- md.1k[md.1k$passed_filters > 500, ]
md.5k <- md.5k[md.5k$passed_filters > 500, ]
md.10k <- md.10k[md.10k$passed_filters > 1000, ] # sequenced deeper so set higher cutoff

# create fragment objects
frags.500 <- CreateFragmentObject(
  path = "pbmc500/atac_pbmc_500_nextgem_fragments.tsv.gz",
  cells = rownames(md.500)
)


frags.1k <- CreateFragmentObject(
  path = "pbmc1k/atac_pbmc_1k_nextgem_fragments.tsv.gz",
  cells = rownames(md.1k)
)

frags.5k <- CreateFragmentObject(
  path = "pbmc5k/atac_pbmc_5k_nextgem_fragments.tsv.gz",
  cells = rownames(md.5k)
)

frags.10k <- CreateFragmentObject(
  path = "pbmc10k/atac_pbmc_10k_nextgem_fragments.tsv.gz",
  cells = rownames(md.10k)
)

在各数据集中对峰值进行量化

利用 FeatureMatrix 函数,现在能够为每个样本生成一个以峰值和细胞为维度的矩阵。此函数通过 future 包支持并行计算。

pbmc500.counts <- FeatureMatrix(
  fragments = frags.500,
  features = combined.peaks,
  cells = rownames(md.500)
)

pbmc1k.counts <- FeatureMatrix(
  fragments = frags.1k,
  features = combined.peaks,
  cells = rownames(md.1k)
)

pbmc5k.counts <- FeatureMatrix(
  fragments = frags.5k,
  features = combined.peaks,
  cells = rownames(md.5k)
)

pbmc10k.counts <- FeatureMatrix(
  fragments = frags.10k,
  features = combined.peaks,
  cells = rownames(md.10k)
)

总结

本文[1]提供了一个详细的流程来合并单细胞染色质数据集,包括数据下载、预处理、合并以及后续的分析和可视化步骤。强调了在合并过程中创建共有峰值集合的重要性,并提供了在没有片段文件时的替代方法。

Reference
[1]

Source: https://stuartlab.org/signac/articles/merging

本文由 mdnice 多平台发布

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/52540.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

SAP与生产制造MPM系统集成案例

一、需求介绍 某公司为保证企业内部生产管理系统的多项基础数据的同步更新&#xff0c;确保各模块间信息的一致性和准确性&#xff0c;对后续的生产计划和物料管理打下基础&#xff0c;该公司将MPM系统和SAP系统经过SAP PO中间件集成平台进行了集成。MPM全称为Manufacturing…

超实用的8个无版权、免费、高清图片素材网站整理

不管是设计、文章配图&#xff0c;还是视频制作&#xff0c;图片都至关重要。但是图片版权一直都是困扰很多设计、自媒体以及企业的大问题。现在&#xff0c;因为图片侵权被告的案例已经是司空见惯了&#xff0c;有的公众号甚至因为图片版权问题遭受致命打击。 1. Pexels Pexe…

Spring Boot 全局异常@ControllerAdvice和@RestControllerAdvice的区别

在Spring Boot中&#xff0c;ControllerAdvice和RestControllerAdvice都可以用于实现全局异常处理&#xff0c;但它们在处理方式和返回值类型上略有不同。至于为什么在某些情况下可能更偏向于使用RestControllerAdvice而不是ControllerAdvice&#xff0c;这主要取决于你的应用类…

前后端传参@RequestParam使用上的一个小坑

RequestParam(value "art") String art 默认情况下&#xff0c;value "art"表示前端传入参数的名字叫art&#xff0c;如果前端不传的话会报错 RequestParam(value "art" required false) String art 当equired false 时表示前端可以不传这个…

Spring框架:开发者的得力伙伴,魅力何在?

目录 一. Spring介绍 二. Spring搭建 三. Spring Bean管理 ▐ 管理方式 ▐ 依赖注入 四. Spring数据访问层管理 五. Spring集成MyBatis 海漫浩浩,我亦苦作舟!大家一起学习,一起进步! 一. Spring介绍 Spring是什么? Spring 是一个轻量级的, IOC 和 AOP 的一站式 J…

2024年最新Java面试宝典系列-Collections集合篇1

Java中的集合类有哪些&#xff1f;它们的特点是什么 List&#xff1a;有序集合&#xff0c;允许重复元素&#xff0c;实现类如ArrayList、LinkedList。Set&#xff1a;无序集合&#xff0c;不允许重复元素&#xff0c;实现类如HashSet、TreeSet。Map&#xff1a;键值对集合&am…

阿里云服务器 篇七:服务器热备份/定时备份

文章目录 系列文章bypy其他备选项目安装bypy使用bypy命令帮助绑定百度账号查看网盘空间大小和账号信息显示在百度网盘中的根目录在百度网盘中创建备份目录百度网盘中的其他文件操作命令从本地目录同步到百度网盘中在云服务器设置定时同步创建云端备份目录为sudo命令做准备创建同…

如何使用ssm实现基于java web的计算机office课程平台设计与实现+vue

TOC ssm277基于java web的计算机office课程平台设计与实现vue 绪论 1.1 研究背景 现在大家正处于互联网加的时代&#xff0c;这个时代它就是一个信息内容无比丰富&#xff0c;信息处理与管理变得越加高效的网络化的时代&#xff0c;这个时代让大家的生活不仅变得更加地便利…

博弈论(Nim游戏的扩展)

公平组合游戏ICG 若一个游戏满足: 1.由两名玩家交替行动; 2.在游戏进程的任意时刻&#xff0c;可以执行的合法行动与轮到哪名玩家无关; 3.不能行动的玩家判负; 则称该游戏为一个公平组合游戏。 NIM博弈属于公平组合游戏&#xff0c;但城建的棋类游戏&#xff0c;比如围棋&…

大刀阔斧改革之后,阅文距离“东方迪士尼”更近了吗?

当前&#xff0c;网文IP的确是“富矿”。中国社会科学院文学研究所发布的《2023中国网络文学发展研究报告》显示&#xff0c;截至2023年底&#xff0c;网络文学IP市场规模2605亿元&#xff0c;同比增长近百亿元。 近日&#xff0c;网文产业中的头部企业阅文集团也披露数据称&a…

虚拟内存和linux(操作系统part1)

一个操作系统的虚拟内存和linux部分知识点的笔记整理&#xff0c;资料大多参考于&#xff1a;小林coding和Javaguide。 虚拟内存的作用 第一&#xff0c;虚拟内存可以使得进程运行内存超过物理内存大小&#xff0c;因为程序运行符合局部性原理&#xff0c;CPU 访问内存会有很…

【iOS安全】iPhone8 iOS14.4.2 越狱教程

环境配置 iPhone 8&#xff1a; 固件版本 iOS 14.4.2 (18D70) 产品类型 iPhone10,1 (A1906) 销售型号 MQ862J/A MacBook Pro&#xff1a; macOS 10.15.7 装有CheckRa1n beta 0.12.4 概述 尝试了几个版本的unc0ver和Taurine&#xff0c;发现都不好使 unc0ver显示unsupported…

Spring Cloud Eureka与Kubernetes的集成:服务发现的混合方案

Spring Cloud Eureka与Kubernetes的集成&#xff1a;服务发现的混合方案 引言 随着微服务架构的流行&#xff0c;服务发现&#xff08;Service Discovery&#xff09;已经成为构建分布式系统的关键组件之一。在分布式系统中&#xff0c;服务实例的数量和位置是动态变化的&…

95.SAP MII功能详解(08)Workbench-Transaction介绍

目录 1.Transaction 2.Properties of transaction 1.Transaction You use transactions to access data from multiple sources and execute processes, which are triggered synchronously or asynchronously.您可以使用事务从多个源访问数据并执行同步或异步触发的流程。…

React Hooks 的使用场景有哪些?

React Hooks是React 16.8引入的一项特性&#xff0c;它允许你在不编写类组件的情况下使用state和其他React特性。以下是React Hooks的一些主要使用场景&#xff1a; 状态管理&#xff1a;使用useState Hook在函数组件中添加本地状态。 副作用处理&#xff1a;使用useEffect Ho…

代码随想录——两个字符串的删除操作(Leetcode 583)

题目链接 动态规划 思路&#xff1a; 确定dp数组&#xff08;dp table&#xff09;以及下标的含义 dp[i][j]&#xff1a;以i-1为结尾的字符串word1&#xff0c;和以j-1位结尾的字符串word2&#xff0c;想要达到相等&#xff0c;所需要删除元素的最少次数。 确定递推公式 当…

第 4 章 ECMASript 8 新特性

4.1 async 和 await async 和 await 两种语法结合可以让异步代码像同步代码一样 4.1.1 async 函数 async 函数的返回值为 promise 对象&#xff0c;promise 对象的结果由 async 函数执行的返回值决定 如果返回的是一个成功的Promise&#xff0c;result的结果就是一个成功的…

PD取电快充协议方案

PD快充协议是通过调整电压和电流来提供不同的充电功率。它采用了一种基于USB-C端口的通信协议&#xff0c;实现了充电器于设备之间的信息交换。在充电过程中设备会向充电器发出请求&#xff0c;要求提供不同的电压和电流&#xff0c;充电器接收到请求后&#xff0c;会根据设备的…

算法训练营|图论第一天 98. 所有可达路径

题目&#xff1a;所有可到达路径 题目链接&#xff1a; 98. 所有可达路径 (kamacoder.com) 解题思路&#xff1a; 邻接矩阵&#xff0c;注意没有result -1时候特判 #include<bits/stdc.h> using namespace std; vector<vector<int>>result; vector<…

cordova手动更新

1&#xff1a;依赖 cordova-plugin-file cordova-plugin-file-transfer cordova-plugin-file-opener2 第二个参数&#xff1a;application/vnd.android.package-archive来源 cordova plugin add cordova-plugin-app-version//获取cordova版本号 cordova plugin add cordova-p…