专栏首页单细胞天地Seurat教程上新||Mixscape : 用多模态单细胞数据筛选免疫检查点

Seurat教程上新||Mixscape : 用多模态单细胞数据筛选免疫检查点

作者 | 周运来

男,

一个长大了才会遇到的帅哥,

稳健,潇洒,大方,靠谱。

一段生信缘,一棵技能树,

一枚大型测序工厂的螺丝钉,

一个随机森林中提灯觅食的津门旅客。

PD-L1等抑制性免疫检查点分子的表达在人类癌症中较为常见,可导致T细胞介导的免疫应答的抑制。在这里,我们应用ECCITE-seq技术来探索调控PD-L1表达的分子网络。ECCITE-seq技术将混合的CRISPR筛查与单细胞mRNA和表面蛋白测量相结合。我们还开发了一个计算框架,mixscape,它通过识别和去除混杂的变异源,大大提高了单细胞扰动屏幕的信噪比。利用这些工具,我们识别和验证PD-L1的调控因子,并利用我们的多模态数据识别转录和转录后的调控模式。特别是,我们发现kelch样蛋白keap1和转录激活因子NRF2介导了IFN刺激后PD-L1的上调。我们的结果为免疫检查点的调节确定了一个新的机制,并为分析多模态单细胞perturbation screens提供了一个强大的分析框架 。

免疫检查点(IC)分子调节免疫反应中激活和抑制之间的关键平衡。在正常生理条件下,抑制IC分子对于维持自身耐受和防止自身免疫是必不可少的,但在人类癌症中,抑制IC分子的表达常常被错误调控,以逃避免疫监控。例如,抑制性IC PD-L1与T细胞上的PD -1受体相互作用,抑制T细胞活化,在许多癌症中过表达,并作为患者生存和免疫治疗反应的预后因素。因此,人们不仅对识别阻断这些相互作用的治疗途径感兴趣,而且对理解癌细胞上调PD-L1等ICs的分子网络也很感兴趣。

之前的研究已经初步建立了一套能够影响PD-L1 mRNA和表面蛋白水平的分子调控因子。大量研究发现,干扰素(IFN)的暴露可迅速诱导PD-L1在体外和肿瘤微环境中的表达。因此,IFN反应的核心成分是PD-L1表达的上游调控因子,包括转录因子IRF1(直接与PD-L1启动子结合)、JAK-STATsignal转导通路和IFN受体本身。此外,还鉴定了其他IFN的阻断信号、PD-L1启动子染色质状态或对uv介导的应激的调节因子。

此外,最近人们特别关注PD-L1稳定性和降解的转录后调节因子的特性。例如,Cullin 3-SPOPE3-ligase复合物可以细胞周期依赖的方式直接泛素化PD-L1,导致其降解。此外,全基因组CRISPR筛选发现了两个之前未鉴定的调节因子cmtm6和CMTM4,它们通过防止溶酶体介导的降解来稳定PD-L1表面表达。在这些案例中,PD-L1调节剂的干扰被证明可以调节抗肿瘤T细胞的活性,这就突出了理解抑制IC分子调控的治疗兴趣。

我们最近引入了扩展的crispr兼容的CITE-seq (ECCITE-seq),它同时测量转录组、表面蛋白水平和在单细胞分辨率上的扰动,作为一种识别和表征分子调节剂的新方法。ECCITE-seq建立在混合CRISPR屏幕的实验设计上,在单一实验中,多个扰动被复用在一起,但是有明显的优势。首先,单细胞测序读数(即 Perturb-seq, CROP-seq, CRISP-seq)能够测量详细的分子表型,而不是单个表型(单个蛋白的表达或细胞活力)。其次,通过同时耦合测量mRNA、表面蛋白和直接检测同一细胞内的导联,ECCITE-seq允许对每个扰动进行多模态表征。因此,我们认为ECCITE-seq将使我们能够同时测试和识别IC分子的新调控因子,特别是区分转录模式和转录后模式。此外,丰富和高维的读出容易促进网络和基于路径的分析,这可以超越鉴定单个基因和产生对其调控机制的洞察。

在这里,我们应用ECCITE-seq来同时扰乱和表征假定的调节抑制分子对IFN不刺激的反应。当分析我们的单细胞数据时,我们确定了混杂的异质性来源,包括那些接受了靶向引导RNA但没有表现出干扰效应的细胞,这些细胞在下游分析中引入了大量的噪声。我们开发并验证了计算方法来控制这些因素,并大幅增加了我们的统计能力来表征多模态扰动。

利用这些工具,我们确定了一组基因,其扰动会影响PD-L1转录水平、表面蛋白水平,或两者都影响,并确定了每个调控器所使用的潜在分子通路。特别是,我们发现在人类癌症中经常突变的kelch样蛋白keap1和转录激活因子NRF2可以改变PD-L1水平。我们将这些发现与CUL3的一个新的调节机制联系起来,并表明该基因通过稳定NRF2通路作为PD-L1mRNA的间接转录激活剂。综上所述,我们的发现确定了免疫检查点调节的一个重要途径,并提出了一个强大的和广泛适用的分析框架来分析ECCITE-seq数据。

在ECCITE-seq实验中,我们运行了10x Genomics 5’(Chromium Single Cell Immune Profiling Solution v1.0, #1000014, #1000020, #1000151)的8个lanes ,目标是每个lanes 捕获10000个细胞。在运行之前,细胞存活率被确定和细胞数量估计之前描述。为了跟踪每个生物学重复身份,样本按照细胞hashing protocol 进行。

  • mRNA
  • hashtags (Hashtag-derived oligos, HTOs)
  • protein (Antibody-derived oligos, ADTs)
  • gRNA (Guide-derived oligos, GDOs)

文库均由10x genomics and ECCITE-seq protocols方法构建。在NovaSeq运行时,所有库都在两个lanes上进行测序。利用Cellranger (V2.1.1)将来自mRNA文库的测序片段映射到hg19参考基因组。为了生成HTO、ADT和GDO库的计数矩阵,使用了CITE-seq-count package (https://github.com/Hoohm/CITE-seq-Count)。然后使用计数矩阵作为Seurat R包的输入来执行所有的下游分析。

下面我们跟着官网教程来看看是如何达到目的的。

首先,我们需要安装新的Seurat,下载示例数据:

remotes::install_github(“satijalab/seurat”, ref = “mixscape”)

# Load packages.
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(scales)
library(dplyr)

# Download dataset using SeuratData.
InstallData(ds = "thp1.eccite")

# Setup custom theme for plotting.
custom_theme <- theme(plot.title = element_text(size = 16, hjust = 0.5), legend.key.size = unit(0.7,
    "cm"), legend.text = element_text(size = 14))

老规矩,我们来看一下数据格式。

# Load object.
eccite <- LoadData(ds = "thp1.eccite")
eccite
An object of class Seurat
18776 features across 20729 samples within 4 assays
Active assay: RNA (18649 features, 0 variable features)
 3 other assays present: ADT, HTO, GDO

我们看到有4个assay: RNA , ADT, HTO, GDO,一个Seurat玩转4套数据,在才叫多模态啊。我们分别看看这四套数据的内容:

eccite@assays$RNA@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
              l1_AAACCTGAGCCAGAAC l1_AAACCTGAGTGGACGT l1_AAACCTGCATGAGCGA l1_AAACCTGTCTTGTCAT
AL627309.1                      .                   .                   .                   .
AP006222.2                      .                   .                   .                   .
RP4-669L17.10                   .                   .                   .                   .
RP11-206L10.3                   .                   .                   .                   .
> eccite@assays$ADT@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
      l1_AAACCTGAGCCAGAAC l1_AAACCTGAGTGGACGT l1_AAACCTGCATGAGCGA l1_AAACCTGTCTTGTCAT
CD86                  118                 243                  99                 110
PDL1                  522                 139                 109                 334
PDL2                   94                 104                  56                  48
CD366                  67                  59                  80                  47
> eccite@assays$HTO@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
          l1_AAACCTGAGCCAGAAC l1_AAACCTGAGTGGACGT l1_AAACCTGCATGAGCGA l1_AAACCTGTCTTGTCAT
rep1-tx                    82                  18                  55                  14
rep1-ctrl                   3                   1                   4                   .
rep2-tx                    13                  11                   6                   6
rep2-ctrl                   1                   4                   .                   .
> eccite@assays$GDO@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
       l1_AAACCTGAGCCAGAAC l1_AAACCTGAGTGGACGT l1_AAACCTGCATGAGCGA l1_AAACCTGTCTTGTCAT
eGFPg1                   1                   1                   1                   1
CUL3g1                   1                   1                   1                   1
CUL3g2                   1                   1                   1                   1
CUL3g3                   1                   1                   1                   1

metadata的内容:

 head(eccite@meta.data)
                    orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO nCount_GDO nFeature_GDO nCount_ADT nFeature_ADT percent.mito MULTI_ID MULTI_classification
l1_AAACCTGAGCCAGAAC      Lane1      17207         3942         99            4        576          111        801            4     2.295577  rep1-tx              rep1-tx
l1_AAACCTGAGTGGACGT      Lane1       9506         2948         35            5        190          111        545            4     4.512939  rep1-tx              rep1-tx
l1_AAACCTGCATGAGCGA      Lane1      15256         4258         66            4        212          111        344            4     4.116413  rep1-tx              rep1-tx
l1_AAACCTGTCTTGTCAT      Lane1       5135         1780         22            3        243          111        539            4     5.491723  rep1-tx              rep1-tx
l1_AAACGGGAGAACAACT      Lane1       9673         2671         99            5        198          111       1053            4     3.359868  rep1-tx              rep1-tx
l1_AAACGGGAGACAGAGA      Lane1      14941         3918         97            5        124          111        487            4     3.379961  rep1-tx              rep1-tx
                    MULTI_classification.global HTO_classification guide_ID guide_ID.global  gene con      NT    crispr replicate     S.Score   G2M.Score Phase
l1_AAACCTGAGCCAGAAC                     rep1-tx            rep1-tx  STAT2g2         Singlet STAT2  tx STAT2g2 Perturbed      rep1 -0.25271565 -0.77130934    G1
l1_AAACCTGAGTGGACGT                     rep1-tx            rep1-tx   CAV1g4         Singlet  CAV1  tx  CAV1g4 Perturbed      rep1 -0.12380192 -0.33260303    G1
l1_AAACCTGCATGAGCGA                     rep1-tx            rep1-tx  STAT1g2         Singlet STAT1  tx STAT1g2 Perturbed      rep1 -0.15463259 -0.69441836    G1
l1_AAACCTGTCTTGTCAT                     rep1-tx            rep1-tx   CD86g1         Singlet  CD86  tx  CD86g1 Perturbed      rep1 -0.06126198 -0.03781951    G1
l1_AAACGGGAGAACAACT                     rep1-tx            rep1-tx   IRF7g2         Singlet  IRF7  tx  IRF7g2 Perturbed      rep1 -0.13218850 -0.35315597    G1
l1_AAACGGGAGACAGAGA                     rep1-tx            rep1-tx     NTg1         Singlet    NT  tx      NT        NT      rep1 -0.20998403 -0.58247261    G1

禁不住用我们的桑吉图看看分属关系:

library(ggforce)
eccite@meta.data %>%
  gather_set_data(17:21) %>%
  ggplot(aes(x, id = id, split = y, value = 1))  +
  geom_parallel_sets(aes(fill = NT), show.legend = FALSE, alpha = 0.3) +
  geom_parallel_sets_axes(axis.width = 0.1, color = "lightgrey", fill = "white") +
  geom_parallel_sets_labels(angle = 0) +
  theme_no_axes()

要理解这里的每一列是什么就要学习CRISPR的基本知识了。重要的一点是理解perturbation 的概念。

# Normalize protein.
eccite <- NormalizeData(object = eccite, assay = "ADT", normalization.method = "CLR", margin = 2)

为了获得全局的观点,我们先对RNA的数据执行Seurat的一般流程:基于rna的聚类是由混杂的变异源驱动的。

# Prepare RNA assay for dimensionality reduction: Normalize data, find variable features and
# scale data.
DefaultAssay(object = eccite) <- "RNA"
eccite <- NormalizeData(object = eccite) %>% FindVariableFeatures() %>% ScaleData()

# Run Principle Component Analysis (PCA) to reduce the dimensionality of the data.
eccite <- RunPCA(object = eccite)

# Run Uniform Manifold Approximation and Projection (UMAP) to visualize clustering in 2-D.
eccite <- RunUMAP(object = eccite, dims = 1:40)

# Generate plots to check if clustering is driven by biological replicate ID, cell cycle phase
# or target gene class.
p1 <- DimPlot(eccite, group.by = "replicate", label = F, pt.size = 0.2, reduction = "umap") + scale_color_brewer(palette = "Dark2") +
    ggtitle("Biological Replicate") + xlab("UMAP 1") + ylab("UMAP 2") + custom_theme

p2 <- DimPlot(eccite, group.by = "Phase", label = F, pt.size = 0.2, reduction = "umap") + ggtitle("Cell Cycle Phase") +
    ylab("UMAP 2") + xlab("UMAP 1") + custom_theme

p3 <- DimPlot(eccite, group.by = "crispr", pt.size = 0.2, reduction = "umap", split.by = "crispr",
    ncol = 1, cols = c("grey39", "goldenrod3")) + ggtitle("Perturbation Status") + ylab("UMAP 2") +
    xlab("UMAP 1") + custom_theme

# Visualize plots.
((p1/p2 + plot_layout(guides = "auto")) | p3)

接下来,计算局部扰动特征减轻混杂效应。

 Calculate local perturbation signature.
eccite <- CalcPerturbSig(object = eccite, assay = "RNA", slot = "data", gd.class = "gene", nt.cell.class = "NT",
    reduction = "pca", ndims = 40, num.neighbors = 20, split.by = "replicate", new.assay.name = "PRTB")

# Prepare PRTB assay for dimensionality reduction: Normalize data, find variable features and
# center data.
DefaultAssay(object = eccite) <- "PRTB"

# Use variable features from RNA assay.
VariableFeatures(object = eccite) <- VariableFeatures(object = eccite[["RNA"]])
eccite <- ScaleData(eccite, do.scale = F, do.center = T)

# Run PCA to reduce the dimensionality of the data.
eccite <- RunPCA(object = eccite, reduction.key = "prtbpca", reduction.name = "prtbpca")

# Run UMAP to visualize clustering in 2-D.
eccite <- RunUMAP(object = eccite, dims = 1:40, reduction = "prtbpca", reduction.key = "prtbumap",
    reduction.name = "prtbumap")

# Generate plots to check if clustering is driven by biological replicate ID, cell cycle phase
# or target gene class.
q1 <- DimPlot(eccite, group.by = "replicate", reduction = "prtbumap", pt.size = 0.2) + scale_color_brewer(palette = "Dark2") +
    ggtitle("Biological Replicate") + ylab("UMAP 2") + xlab("UMAP 1") + custom_theme

q2 <- DimPlot(eccite, group.by = "Phase", reduction = "prtbumap", pt.size = 0.2) + ggtitle("Cell Cycle Phase") +
    ylab("UMAP 2") + xlab("UMAP 1") + custom_theme

q3 <- DimPlot(eccite, group.by = "crispr", reduction = "prtbumap", split.by = "crispr", ncol = 1,
    pt.size = 0.2, cols = c("grey39", "goldenrod3")) + ggtitle("Perturbation Status") + ylab("UMAP 2") +
    xlab("UMAP 1") + custom_theme

# Visualize plots.
(q1/q2 + plot_layout(guides = "auto") | q3)

Mixscape识别没有检测扰动的细胞。

#install.packages('mixtools')
# Run mixscape to classify cells based on their perturbation status.
eccite <- RunMixscape(object = eccite, assay = "PRTB", slot = "scale.data", labels = "gene", nt.class.name = "NT",
    min.de.genes = 5, iter.num = 10, de.assay = "RNA", verbose = F)

# Show that only IFNG pathway KO cells have a reduction in PD-L1 protein expression.
Idents(object = eccite) <- "gene"

VlnPlot(object = eccite, features = "adt_PDL1", idents = c("NT", "JAK2", "STAT1", "IFNGR1", "IFNGR2",
    "IRF1"), group.by = "gene", pt.size = 0.2, sort = T, split.by = "mixscape_class.global", cols = c("coral3",
    "grey79", "grey39")) + ggtitle("PD-L1 protein") + theme(axis.text.x = element_text(angle = 0,
    hjust = 0.5))

KO,NP,NT 是啥意思?

  • If the cell received a NT gRNA, it retains its assignment as non-targeting (NT)
  • If the cell received a targeting gRNA, and is classified in step 8 as NT, it is assigned a non-perturbed (NP) label
  • If the cell received a targeting gRNA, and is classified in step 8 as KO, it receives a perturbed/knock-out (KO) label

用线性判别分析(LDA)可视化扰动响应。

# Remove non-perturbed cells and run LDA to reduce the dimensionality of the data.
Idents(eccite) <- "mixscape_class.global"
sub <- subset(eccite, idents = c("KO", "NT"))

sub <- MixscapeLDA(object = sub, assay = "RNA", pc.assay = "PRTB", labels = "gene", nt.label = "NT",
    npcs = 10, logfc.threshold = 0.25, verbose = F)

# Use LDA results to run UMAP and visualize cells on 2-D.
sub <- RunUMAP(sub, dims = 1:11, reduction = "lda", reduction.key = "ldaumap", reduction.name = "ldaumap")

# Visualize UMAP clustering results.
Idents(sub) <- "mixscape_class"
sub$mixscape_class <- as.factor(sub$mixscape_class)
p <- DimPlot(sub, reduction = "ldaumap", label = T, repel = T, label.size = 5)

col = setNames(object = hue_pal()(12), nm = levels(sub$mixscape_class))
names(col) <- c(names(col)[1:7], "NT", names(col)[9:12])
col[8] <- "grey39"

p + scale_color_manual(values = col, drop = FALSE) + ylab("UMAP 2") + xlab("UMAP 1") + custom_theme

代码和降维图均来自mixscape官网,为什么没有自己画,因为画到一般发现自己的PC计算资源不够用了。对原理感兴趣看一看文章:Characterizing the molecular regulation of inhibitory immune checkpoints with multi-modal single-cell screens.

在分析的过程中注意Seurat数据的assay之间的切换,这是四套数据了。大部分函数是我们熟悉的Seurat的函数,几个关键的函数需要我们亲自查看help文档,如RunMixscape ``CalcPerturbSig。当然,最为关键的还是我们对ECCITE-seq和CRISPR技术的原理和细节。

Mixscape 的出现再次验证了:单细胞只是概念,我们可以基于此开发相应的技术。


References

[1] mixscape: https://github.com/satijalab/seurat/tree/mixscape [2] 基因筛查进入单细胞时代: https://www.lifeomics.com/?p=36110 [3] CRISPR-Cas9基因编辑技术简介: https://zhuanlan.zhihu.com/p/137760447

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:周运来

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-07-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 综述:高维单细胞RNA测序数据分析工具(中)

    当你的才华还撑不起你的野心时,请潜下心来,脚踏实地,跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了,通过文献速递这个栏目很幸运聚集了一些小伙伴...

    生信技能树jimmy
  • 聚类算法之PCA与tSNE

    没有体现任何的基因差异或者样本聚类(热图中的聚类是自然层次聚类),可以看到样本名都是无规律的交叉显示

    生信技能树jimmy
  • dyno使用教程--1个R包实现59种单细胞轨迹推断分析

    A comparison of single-cell trajectory inference methods[1](Saelens et al., 2019...

    生信技能树jimmy
  • 修改CSS ,解决帝国分页错位的问题

    用户1272546
  • CentOS yum安装或者升级GCC到4.8

    之前一直没怎么用C++11,最近的项目都开始使用,但是Centos使用yum工具安装GCC的版本是4.4.7,因此需要升级才行。

    帘卷西风
  • node.js中的nodemon详解及安装方法

    nodemon是一种工具,可在检测到目录中的文件更改时通过自动重新启动节点应用程序来帮助开发基于node.js的应用程序。

    逍遥灰灰
  • IT文档中总出现的"hooks" 是什么? "钩子" ? 回调 ?

    就像一些外来的钩子,在源代码之间钩取一些信息,当捕捉到感兴趣的事时,就拦截下来,让自己的代码执行一下,处理一下这个信息,然后再放出去继续之前的进程。 这样就可...

    JavaEdge
  • OpenGL与CUDA互操作方式总结

    一、介绍 CUDA是Nvidia推出的一个通用GPU计算平台,对于提升并行任务的效率非常有帮助。本人主管的项目中采用了OpenGL做图像渲染,但是在数据处理方面...

    24K纯开源
  • UOS系统安装-配置Qt环境

    1、UOS系统相对来说,较为纯净,像GCC、G++等编译器是没有安装的,需要自己去安装

    Sky_Mao
  • 被黑的Drupal网站被用来挖矿,传播远控,发送诈骗邮件

    两个漏洞编号为CVE-2018-7600和CVE-2018-7602的远程代码执行漏洞之前已经被Drupal开发人员修复。

    FB客服

扫码关注云+社区

领取腾讯云代金券