首页 > 其他分享 >Signac包-2.联合10x多组学分析:scATAC-seq和scRNA-seq

Signac包-2.联合10x多组学分析:scATAC-seq和scRNA-seq

时间:2024-08-05 21:27:06浏览次数:14  
标签:Signac seq scRNA ATAC 基因 pbmc sorted RNA


–https://stuartlab.org/signac/articles/pbmc_multiomic

看文章看累了来看看代码,换换口味。本章主要涉及peaks to genes的联动。

留意更多内容,欢迎关注微信公众号:组学之心

数据下载:

wget https://cf.10xgenomics.com/samples/cell-arc/1.0.0/pbmc_granulocyte_sorted_10k/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-arc/1.0.0/pbmc_granulocyte_sorted_10k/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-arc/1.0.0/pbmc_granulocyte_sorted_10k/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz.tbi

1.导入scATCA-seq和scRNA-seq数据

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)

counts <- Read10X_h5("../vignette_data/multiomic/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")
fragpath <- "../vignette_data/multiomic/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"

从hg38中获取细胞注释

下载获取基因注释文件,包含有关基因组位置的信息,例如染色体、起始和结束位置,以及相关的meta数据。

annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
#给染色体加上前缀“chr”
seqlevels(annotation) <- paste0('chr', seqlevels(annotation))

先创建scRNA-seq seurat项目

pbmc <- CreateSeuratObject(
  counts = counts$`Gene Expression`,
  assay = "RNA"
)

再创建scATAC-seq assay,然后把它加入到seurat项目中

pbmc[["ATAC"]] <- CreateChromatinAssay(
  counts = counts$Peaks,
  sep = c(":", "-"),
  fragments = fragpath,
  annotation = annotation
)
pbmc

## An object of class Seurat 
## 144978 features across 11909 samples within 2 assays 
## Active assay: RNA (36601 features, 0 variable features)
##  1 layer present: counts
##  1 other assay present: ATAC

2.计算质控指标

我们可以使用 DNA 可及性数据计算每个细胞的质量控制指标,并去除这些指标的异常细胞,以及 RNA 或 ATAC 测序中count低或异常高的细胞。

DefaultAssay(pbmc) <- "ATAC"

# 计算每个细胞的nucleosome signal值
pbmc <- NucleosomeSignal(pbmc)

# 计算每个细胞的TSS enrichment值
pbmc <- TSSEnrichment(pbmc)

可以使用 DensityScatter() 函数可视化存储在项目中metadata中的变量之间的关系。通过设置 quantiles=TRUE,还可以快速找到不同 QC 指标的合适截止值:

DensityScatter(pbmc, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

VlnPlot(
  object = pbmc,
  features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal"),
  ncol = 4,
  pt.size = 0
)

小提琴图指控指标的统计分布。用以下阈值来筛选(没有最优解,只有适合与否):

pbmc <- subset(
  x = pbmc,
  subset = nCount_ATAC < 100000 &
    nCount_RNA < 25000 &
    nCount_ATAC > 1800 &
    nCount_RNA > 1000 &
    nucleosome_signal < 2 &
    TSS.enrichment > 1
)
pbmc

## An object of class Seurat 
## 144978 features across 11070 samples within 2 assays 
## Active assay: ATAC (108377 features, 0 variable features)
##  2 layers present: counts, data
##  1 other assay present: RNA

3.基因表达数据预处理

可以使用 SCTransform 对基因表达数据进行归一化,并使用 PCA 降低维数

DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc)
pbmc <- RunPCA(pbmc)

4.染色质可及性数据处理

TF-IDF 和 SVD 的组合步骤称为潜在语义索引 (LSI)。通过LSI,以与处理 scATAC-seq 数据集相同的方式处理 DNA 可及性。

DefaultAssay(pbmc) <- "ATAC"
pbmc <- FindTopFeatures(pbmc, min.cutoff = 5)
pbmc <- RunTFIDF(pbmc)
pbmc <- RunSVD(pbmc)

5.注释细胞类型

将使用 Hao 等人 (2020) 的带注释的 PBMC 参考数据集,可以从这里下载:https://atlas.fredhutch.org/data/nygc/multimodal/pbmc_multimodal.h5seurat

加载PBMC的参考数据

reference <- LoadH5Seurat("00practice/PBMC10K_scATAC_scRNA/pbmc_multimodal.h5seurat", assays = list("SCT" = "counts"), reductions = 'spca')
reference <- UpdateSeuratObject(reference)
DefaultAssay(pbmc) <- "SCT"

细胞类型注释转移

transfer_anchors <- FindTransferAnchors(
  reference = reference,
  query = pbmc,
  normalization.method = "SCT",
  reference.reduction = "spca",
  recompute.residuals = FALSE,
  dims = 1:50
)

predictions <- TransferData(
  anchorset = transfer_anchors, 
  refdata = reference$celltype.l2,
  weight.reduction = pbmc[['pca']],
  dims = 1:50
)

pbmc <- AddMetaData(
  object = pbmc,
  metadata = predictions
)

设置细胞类型为预测值,并去除低质量的预测

Idents(pbmc) <- "predicted.id"
pbmc <- pbmc[, pbmc$prediction.score.max > 0.5]

6.联合UMAP图展示

pbmc <- FindMultiModalNeighbors(
  object = pbmc,
  reduction.list = list("pca", "lsi"), 
  dims.list = list(1:50, 2:40),
  modality.weight.name = "RNA.weight",
  verbose = TRUE
)

注意的是ATAC数据dims.list从第2个开始,第一个 LSI component 通常捕获测序深度(技术变化)而不是生物变化。

#构建联合UMAP图
pbmc <- RunUMAP(
  object = pbmc,
  nn.name = "weighted.nn",
  assay = "RNA",
  verbose = TRUE
)

DimPlot(pbmc, label = TRUE, repel = TRUE, reduction = "umap") + NoLegend()

7.将峰值与基因联系起来

对于每个基因,signac可以通过计算基因表达与附近峰值的可及性之间的相关性,并校正由于 GC 含量、总体可及性和峰值大小而导致的偏差,来找到可能调节该基因的峰值。

在整个基因组上运行此步骤会非常耗时,在这里用部分基因的峰值基因关联为例进行演示。通过省略 genes.use 参数,使用相同的函数来查找所有基因的关联:

# 首先计算每个峰的 GC 含量
pbmc <- RegionStats(pbmc, genome = BSgenome.Hsapiens.UCSC.hg38)

# 将峰值与基因联系起来
pbmc <- LinkPeaks(
  object = pbmc,
  peak.assay = "ATAC",
  expression.assay = "SCT",
  genes.use = c("LYZ", "MS4A1")
)

idents.plot <- c("B naive", "B intermediate", "B memory",
                 "CD14 Mono", "CD16 Mono", "CD8 TEM", "CD8 Naive")

作图展示

p1 <- CoveragePlot(
  object = pbmc,
  region = "MS4A1",
  features = "MS4A1",
  expression.assay = "SCT",
  idents = idents.plot,
  extend.upstream = 500,
  extend.downstream = 10000
)

p2 <- CoveragePlot(
  object = pbmc,
  region = "LYZ",
  features = "LYZ",
  expression.assay = "SCT",
  idents = idents.plot,
  extend.upstream = 8000,
  extend.downstream = 5000
)

patchwork::wrap_plots(p1, p2, ncol = 1)

图展示了7种细胞类型种MS4A1和LYZ两个基因的ATAC-seq和RNA-seq情况,中间的小山峰是ATAC-seq数据的展示,代表了染色质开放的位置区域,每个基因名字下的深蓝色条带显示该基因的结构和转录方向。深蓝色粗矩形代表外显子;深蓝色细线代表内含子;箭头表示转录方向(RNA聚合酶的运动方向)。

连线部分是小山峰之间的相关性预测,代表它们之间相互作用的强度,score越高,互作越强。后续可以进一步做FIMO分析,来给小山峰注释上调控元件,进一步探究调控元件互作的内容/方式。

右侧的小提琴图展示该基因在不同细胞亚群中的表达情况(RNA)。

在这里插入图片描述

标签:Signac,seq,scRNA,ATAC,基因,pbmc,sorted,RNA
From: https://blog.csdn.net/weixin_56751316/article/details/140889417

相关文章

  • Pytorch笔记|小土堆|P16-22|神经网络基本骨架、卷积层、池化层、非线性激活层、归一化
    torch.nnContainers是神经网络骨架,含6个类,最常用的是Module——BaseclassforallNNmodulesModule所有神经网络模型(子类)都必须继承Module(父类),Module相当于给所有的神经网络提供了模板,但可进行修改官方示例:importtorch.nnasnnimporttorch.nn.functionalasFclass......
  • 易基因:MeRIP-seq+RNA-seq揭示家禽(鸡)脂肪沉积中的m6A RNA甲基化调控机制|项目文章
    大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。脂肪组织是主要的脂肪沉积和代谢场所,在人类和动物的健康、代谢平衡和免疫稳态中起着重要作用。肥胖引起的多种代谢疾病已成为全球性健康问题。家禽腹部脂肪的过量沉积也会导致代谢疾病和增加饲料浪费,从而增加家禽生产......
  • 洛谷题单指南-前缀和差分与离散化-P4552 [Poetize6] IncDec Sequence
    原题链接:https://www.luogu.com.cn/problem/P4552题意解读:对一组数字序列,进行若干次区间+1或者-1操作,最终使得所有数字一样,计算最少的操作次数,以及能得到多少种不同序列。解题思路:要使得序列每一个数字都相同,则其差分除了第一项之外其余项都是0。因此,问题转化为:给定一个差分数......
  • nodejs 使用 sequelize 实现 mysql数据库的批量插入
    直接上代码:/***设置mysql连接,返回连接实例。连接格式:账户:密码@数据库地址/具体数据库名称***/constsetConnect=()=>{constsequelize=newSequelize(`mysql://${你的mysql地址}`,{logging:(...msg)=>Logger.INSTANCE.inf......
  • Invertible Bracket Sequences
    看到这种类似的括号匹配的题目,一定要想到卡特兰数的证明过程呀(将(看成\(1\),)看成\(-1\),于是不难得出充分条件,虽然这道题目并不是直接这么给的,但是我看没人证明)剩下的看官方题解就好了,之所以可以删掉官方题解所说的\(x\),是因为接下来如果\(x\)是\(r_1\)的答案候选项的话,由于\(r_1>......
  • macOS Sequoia 15.1 beta (24B5009l) Boot ISO 原版可引导镜像下载
    macOSSequoia15.1beta(24B5009l)BootISO原版可引导镜像下载iPhone镜像、Safari浏览器重大更新、备受瞩目的游戏和AppleIntelligence等众多全新功能令Mac使用体验再升级请访问原文链接:https://sysin.org/blog/macOS-Sequoia-boot-iso/,查看最新版。原创作品,转载请保......
  • macOS Sequoia 15.1 beta (24B5009l) ISO、IPSW、PKG 下载
    macOSSequoia15.1beta(24B5009l)ISO、IPSW、PKG下载iPhone镜像、Safari浏览器重大更新、备受瞩目的游戏和AppleIntelligence等众多全新功能令Mac使用体验再升级请访问原文链接:https://sysin.org/blog/macOS-Sequoia/,查看最新版。原创作品,转载请保留出处。作者主页......
  • Bracket Sequences II
    原题链接题解一个合法的括号序列,满足长度为偶数前缀和处处不小于0左括号等于右括号数量code#include<bits/stdc++.h>#definelllonglong#definelowbit(x)((x)&(-x))usingnamespacestd;constllinf=1e18;constllmod=1e9+7;llqpow(lla,lln){......
  • DiffSeq
    目录概符号说明流程代码 GongS.,LiM.,FengJ.,WuZ.andKongL.DiffuSeq:Sequencetosequencetextgenerationwithdiffusionmodels.InInternationalConferenceonLearningRepresentations(ICLR),2023概本文提出了一种用于Seq2Seq的不需要clas......
  • 论文阅读:Sequence to sequence learning for joint extraction of entities and relat
    用以解决重叠关系问题GGNNs模型GGNNs(门控图神经网络,GatedGraphNeuralNetworks)是一种处理图结构数据的神经网络模型。它是图神经网络(GNN)的一个变体,使用了类似于长短时记忆网络(LSTM)中的门控机制来更有效地处理图中的信息流。GGNNs的核心机制GGNNs的核心思想是通过在图结构中......