ADAR1基因敲除前后肿瘤免疫微环境单细胞水平变化

在单细胞大行其道的近两年，我也安排了学徒们做了几百个有表达量矩阵可以下载的单细胞转录组文献图表复现，挑选其中100个成功的案例，提供代码给大家，希望对大家有帮助！

100个单细胞转录组图表复现

今天要介绍的文章是：Loss of ADAR1 in tumours overcomes resistance to immune checkpoint blockade. Nature 2019 ，它的数据链接是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE110746

可以看到是4个样品：

GSM3015874 Adar_replicate_1
GSM3015875 Adar__replicate_2
GSM3015876 Control_replicate_1
GSM3015877 Control_replicate_2


然后作者提供了这4个10X样品的合并后的3个文件，如下所示

Supplementary file Size Download File type/resource
GSE110746_barcodes.tsv.gz 41.2 Kb (ftp)(http) TSV
GSE110746_genes.tsv.gz 212.7 Kb (ftp)(http) TSV
GSE110746_matrix.mtx.gz 73.2 Mb (ftp)(http) MTX


对于 10x的数据，很简单的 Read10X(data.dir = "GSE110746/") 就可以读取啦！

这个数据的研究目标是：Knockout of Adar1 in B16 melanoma results in global reprogramming of the tumor immune microenvironment

本次研究总共是8000多个单细胞，但是研究者绘制tSNE图的时候，去除了肿瘤细胞，仅仅是展示了7000多个微环境细胞，如下所示：

tSNE图

作者在附件列出来了详细的细胞亚群注释依据：

细胞亚群注释

可以看到每个亚群的特异性基因的高表达情况：

高表达情况热图

开始复现

首先保证 运行环境下面有一个 GSE110746 文件夹，里面有3个文件，如下所示：

$ ls -lh GSE110746
total 295M
-rw-r--r-- 1 win10 197121 230K Feb 17  2018 barcodes.tsv
-rw-r--r-- 1 win10 197121 724K Feb 16  2018 genes.tsv
-rw-r--r-- 1 win10 197121 293M Feb 16  2018 matrix.mtx


数据链接是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE110746

step1构建对象

rm(list=ls())
options(stringsAsFactors = F) 
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(patchwork) # The Composer of Plots
library(Seurat) # Tools for Single Cell Genomics
## 读取数据
scrna_data <- Read10X(data.dir = "GSE110746/")
#构建Seurat对象
seob <- CreateSeuratObject(
  counts = scrna_data, # 表达矩阵可以是稀疏矩阵也可以是普通的矩阵
  min.cells = 3, # 筛选去除小于3个细胞中表达的基因
  min.features = 200) # 去除只有200个一下基因的表达细胞
meat.data <- seob[[]]
dim(seob)
seob[["sample"]] <- gsub(".*-","",rownames(meat.data))
seob = NormalizeData(seob) 
dim(seob)
save(seob,file = "201.RData")


step2进行多样品整合

因为是 4个样品：

GSM3015874 Adar_replicate_1
GSM3015875 Adar__replicate_2
GSM3015876 Control_replicate_1
GSM3015877 Control_replicate_2


这里采用 SCT 方法 ：


# SCT法整合数据
rm(list=ls())
options(stringsAsFactors = F) 
load(file = "201.RData")
# 1. SCTransform
seob_list <- SplitObject(seob, split.by = "sample") # split.by 与52行对应
seob_list
# 查找指定基因
# grep("Mlana",rownames(seob))
# grep("Apoe",rownames(seob))
for(i in 1:length(seob_list)){
  seob_list[[i]] <- SCTransform(
    seob_list[[i]], 
    variable.features.n = 3000,
    # vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"),
    verbose = FALSE)
}

# 2. 数据整合（Integration）
## 选择要用于整合的基因
features <- SelectIntegrationFeatures(object.list = seob_list,
                                      nfeatures = 3000)
## 准备整合
seob_list <- PrepSCTIntegration(object.list = seob_list, 
                                anchor.features = features)
## 找 anchors，15到30分钟
anchors <- FindIntegrationAnchors(object.list = seob_list, 
                                  # reference = 3 # 当有多个样本时，制定一个作为参考可加快速度
                                  normalization.method = "SCT", # 选择方法“SCT”
                                  anchor.features = features
)
## 整合数据
seob <- IntegrateData(anchorset = anchors, 
                      normalization.method = "SCT")
DefaultAssay(seob) <- "integrated"
save(seob,file = "meta_PCA.RData")
# 我们测试了seurat的V3和V4版本的包，跑这个代码，结果居然是不一样的！


上面的代码以 seurat的V3版本为准，如果是V4，后面的聚类分群数量不一致。

多个样品的整合

这里虽然是选择了SCT，其实CCA也可以，或者说harmony也可以，并没有太大的差异哈！仅仅是一个示例哦

step3 降维聚类分群

## 降维分析----
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(patchwork) # The Composer of Plots
library(Seurat) # Tools for Single Cell Genomics
load(file = "meta_PCA.RData")
seob <- RunPCA(seob)  
ElbowPlot(seob, ndims = 50)
## t-SNE 不是基于表达矩阵做的，而是基于PCA结果做的
seob <- RunTSNE(seob, 
                # 使用多少个PC，如果使用的是 SCTransform 建议多一些
                dims = 1:30)

## UMAP
seob <- RunUMAP(seob, dims = 1:30)

# # 检查连续型协变量：线粒体比例
FeaturePlot(seob,
            reduction = "umap", # pca, umap, tsne
            features = "percent.mt")

## 聚类分析 -----
seob <- FindNeighbors(seob,
                      k.param = 5,# 最小值为5
                      dims = 1:30)
seob <- FindClusters(seob,
                     resolution = 0.4, # 分辨率值越大，cluster 越多0~1之间
                     random.seed = 1) #设置随机种子
DimPlot(seob, 
              reduction = "umap", # pca, umap, tsne
              group.by  = "seurat_clusters",
              label = T)

如下所示：

如果是seurat V4 ，同样的代码得到细胞亚群并不一致哦！如果是harmony或者CCA整合，也不一样，感兴趣的可以自行去测试哈！

step4 细胞亚群注释

首先看文章里面的标记基因在各个亚群的表达情况：


library(SummarizedExperiment)
## 构建CellMarker文件
# 这个文件来源于文章，强烈 依赖于生物学背景
library(readr)
cell_marker <- read.table(file = "marker.txt", 
                          sep = "\t",
                          header = T) 
cell_marker
## marker 基因表达

## 可视化
library(tidyr)
cell_marker <- separate_rows(cell_marker, Cell_Marker, sep = ",") %>% # 转为长数据
  distinct() %>% # 去除重复
  arrange(Cell_Type) # 排序

p1 <- DimPlot(seob, 
              reduction = "umap", # pca, umap, tsne
              group.by  = "seurat_clusters",
              label = T)

p2 <- DotPlot(seob, features = unique(cell_marker$Cell_Marker)) +
  theme(axis.text = element_text(size = 8,# 修改基因名的位置
                                 angle = 45,
                                 hjust = 1))
p1 + p2
ggsave(plot = p1,filename = 'umap.pdf')

ggsave(plot = p2,filename = 'markers.pdf')

DotPlot(seob, 
        features = c("Pmel","Mlana")) +
  theme(axis.text = element_text(size = 8,# 修改基因名的位置
                                 angle = 45,
                                 hjust = 1))

如下所示：

合适的标记基因在各个亚群的表达情况

高表达不同基因的亚群可以根据生物学背景进行命名：

## 将细胞类型信息加到meat.data表中
## 构建细胞类型向量

# 这个是肉眼看，依据生物学背景
# 所以每次都需要修改

# 一定要自己肉眼看各个标记基因哦！
# 然后决定下面的分群
# 首先我们看 seurat V4版本的是 0-13的群
grep("ll1b",rownames(seob))
cluster2type <- c("0"="MDSC",
                  "1"="M1 Macrophage",
                  "2"="tumour cells",
                  "3"="Mki67-CD8+T cell",
                  "4"="M2 Macrophage",
                  "5"="Neutrophil",
                  "6"="Mki67+CD8+T cell",
                  "7"="Monocyte",
                  "8"="tumour cells",
                  "9"="CD103+cDC",
                  "10"="Migratory cDC", 
                  "11"="MoDC",
                  "12"="Stromal cell",
                  "13"="Plasmacytoid")

# 然后看seurat V3版本的是 0-18的群


# Seurat_v3
cluster2type <- c("0"="M1 Macrophage",
                  "1"="MDSC",
                  "2"="Neutrophil",
                  "3"="M2 Macrophage",
                  "4"="tumour cells",
                  "5"="Mki67-CD8+T cell",
                  "6"="Mki67+CD8+T cell",
                  "7"="tumour cells",
                  "8"="CD11b+ cDC",
                  "9"="Treg",
                  "10"="tumour cells", 
                  "11"="CD103+cDC",
                  "12"="NK cell",
                  "13"="Migratory cDC",
                  "14"="MoDC",
                  "15"="Stromal cell",
                  "16"="Plasmacytoid DC",
                  "17"="", #未找到高变基因
                  "18"="Stromal cell")


## 可视化
## 二维图
FeaturePlot(seob, 
            reduction = "umap", 
            features = c("Pmel","Mlana"), 
            #split.by = "sample",
            label = TRUE)

## 小提琴图
VlnPlot(object = seob, 
        features = c("Pmel","Mlana")) + plot_layout(guides = "collect") & theme(legend.position = "top")

## 添加细胞类型

seob[['cell_type']] = unname(cluster2type[[e-mail]$seurat_clusters])

DimPlot(seob, 
        reduction = "umap", 
        group.by = "cell_type",
        label = TRUE, pt.size = 1,
        label.size = 6, 
        repel = TRUE)


最后拿到了有生物学意义的亚群：

有生物学意义的亚群

这就是完成了最基础的聚类分群和注释，也可以参考其它例子，比如：人人都能学会的单细胞聚类分群注释

如果你没有单细胞转录组数据分析的基础知识，可以先一下基础10讲：

• 01. 上游分析流程
• 02.课题多少个样品，测序数据量如何
• 03. 过滤不合格细胞和基因（数据质控很重要）
• 04. 过滤线粒体核糖体基因
• 05. 去除细胞效应和基因效应
• 06.单细胞转录组数据的降维聚类分群
• 07.单细胞转录组数据处理之细胞亚群注释
• 08.把拿到的亚群进行更细致的分群
• 09.单细胞转录组数据处理之细胞亚群比例比较

栏目感想

写100个单细胞的图表复现，工作量肯定是不小，幸好我有那么多学徒。确实是非常感谢大家对我们的《生信技能树》的支持，希望这个栏目能在3个月内更新完毕（2021-5到8月），敬请期待哦！