单细胞韧皮部研究代码解析3-comparison_brady.R

原创

小胡子刺猬的生信学习123

修改于 2023-05-05 20:33:24

22200

代码可运行

文章被收录于专栏：文献分享及代码学习文献分享及代码学习

运行总次数：0

代码可运行

单细胞韧皮部研究代码解析1-QC_filtering.R:https://cloud.tencent.com/developer/article/2256814?areaSource=&traceId=

单细胞韧皮部研究代码解析2--comparison_denyer2019.R：https://cloud.tencent.com/developer/beta/article/2260381?areaSource=&traceId=

代码链接：https://github.com/tavareshugo/publication_Otero2022_PhloemPoleAtlas/blob/master/analysis/05-comparison_brady.R

在做单细胞分析的时候，很多研究学者会将以前做的显微切割的RNA-seq的数据与自己的研究内容进行比对，则可以得到在单细胞和亚细胞水平上的数据结果。一般是通过将已经人工注释好的亚群与显微切割数据进行相关性分析，去判断自己的分群是否准确，但是这位作者的代码是可以在进行降维获得亚群后，直接可以与普通的RNA-SEQ数据进行整合分析，为后面的人工亚群注释进行相关的参考。

这篇文章的作者也是通过相关得内容进行了代码演示，通过改写作者的部分内容，是可以做到自己研究的内容的图片可视化的。

今天主要是对上面的图片的代码进行解析。

代码解析

# Setup -------------------------------------------------------------------
# R包的加载
library(scater)
library(tidyverse)
#在R中，经常会需要读入excel文件，这个包也是经常用到的
library(readxl)
library(dplyr)

# set ggplot2 theme
theme_set(theme_bw() + theme(text = element_text(size = 14)))
# 这是作者提的R的函数，在前面的教程里面也提到了这个脚本的下载地址
source("analysis/functions/utils.R")


# Tidy Brady data ------------------------------------------------------
# 这一段是作者要进行比较的数据集，主要是根据自己相关的内容进行更改
# read both sheets
# Brady提供的是probe和gene对应的数据集，是需要将不同的表达组织进行对应的
# 文章中作者选用的数据链接来源：https://science.sciencemag.org/highwire/filestream/588629/field_highwire_adjunct_files/1/1146265BradySTables.zip，可以先自己下载作者的数据进行尝试后，对自己的相关的内容进行改写
brady <- list(
  longitudinal = read_xls("data/external/brady_et_al/ST12.xls",
                          sheet = "LONGITUDINAL", skip = 1) %>%
    pivot_longer(c(-Probe, -Gene),
                 names_to = "slice", values_to = "expr"),
  radial = read_xls("data/external/brady_et_al/ST12.xls",
                    sheet = "RADIAL", skip = 1) %>%
    pivot_longer(c(-Probe, -Gene),
                 names_to = "slice", values_to = "expr")
) %>%
  bind_rows(.id = "set") %>%
  rename_with(tolower) %>%
  # some probes don't have gene
  drop_na(gene)

# annotations for radial sections
radial_annot <- tribble(
  ~slice, ~tissue,
  "AGL42_MEAN", "QC",
  "APL_MEAN", "Phloem + CC",
  "COBL9_MEAN", "Hair cell",
  "CORTEX_MEAN", "Cortex",
  "gl2_MEAN", "Non-hair cell",
  "JO121_MEAN", "Xylem pole pericyle",
  "J0571_MEAN", "Ground tissue",
  "J2661_MEAN", "Mature pericycle",
  "LRC_MEAN", "LRC",
  "pet111_MEAN", "Columella",
  "rm1000_MEAN", "Lateral root",
  "S17_MEAN", "Phloem Pole Pericycle",
  "S18_MEAN", "Maturing xylem",
  "S4_MEAN", "Developing xylem",
  "S32_MEAN", "All protophloem",
  "scr5_MEAN", "Endodermis",
  "SUC2_MEAN", "Phloem CC",
  "wol_MEAN", "Stele",
  "xylem_2501_MEAN", "Stele"
)
## 以上的相关的内容主要是对每一列的值进行宽表改成长表，如何进行表格长宽的改变，也在以前的R语言教程里面有，有不会的友友可以取翻一下以前的教程
# Read SCE data -----------------------------------------------------------

# sce object
# 对降维分析的数据进行加载
ring <- readRDS("data/processed/SingleCellExperiment/ring_batches_strictfilt.rds")


# Longitudinal sections ----------------------------------------------------
# assign each cell to a section in Brady using correlation

## 根据gene内容对两个数据集进行合并
# matrix of expression
lon_matrix <- brady %>%
  filter(set == "longitudinal") %>%
  pivot_wider(names_from = slice, values_from = expr) %>%
  # retain only genes with a single probe
  group_by(gene) %>%
  filter(n() == 1) %>%
  ungroup() %>%
  # retain genes that also occur in our dataset
  filter(gene %in% rownames(ring)) %>%
  # convert to matrix
  select(-set, -probe) %>%
  column_to_rownames("gene") %>%
  as.matrix()

# 选择spearman进行相关系数计算
# correlation matrix
lon_cor <- cor(as.matrix(assay(ring, "logvst")[rownames(lon_matrix), ]),
               lon_matrix,
               method = "spearman")

## 得到每个细胞和brady数据的相关性数据
# get the Brady section with maximum correlation for each cell
lon_classes <- tibble(cell = rownames(lon_cor),
                      section = colnames(lon_cor)[max.col(lon_cor)],
                      cor = rowMax(lon_cor))

## 随后是选用作者写的函数getReducedDim进行绘图
# Plot
ring %>%
  getReducedDim("UMAP30_MNN_logvst") %>%
  mutate(cell = paste(Sample, Barcode, sep = "_")) %>%
  full_join(lon_classes, by = "cell") %>%
  mutate(section = str_remove(section, "SB_MEAN")) %>%
  mutate(section = as.numeric(str_remove(section, "L"))) %>%
  ggplot(aes(V1, V2)) +
  geom_point(data = getReducedDim(ring, "UMAP30_MNN_logvst"), colour = "lightgrey") +
  geom_point(aes(colour = factor(section))) +
  scale_colour_viridis_d() +
  coord_equal() + theme_void() +
  labs(colour = "Brady et. al\nSection")

ring %>%
  getReducedDim("UMAP30_MNN_logvst") %>%
  mutate(cell = paste(Sample, Barcode, sep = "_")) %>%
  full_join(lon_classes, by = "cell") %>%
  mutate(section = str_remove(section, "SB_MEAN")) %>%
  mutate(section = as.numeric(str_remove(section, "L"))) %>%
  ggplot(aes(V1, V2)) +
  geom_point(data = getReducedDim(ring, "UMAP30_MNN_logvst"), colour = "lightgrey") +
  ggpointdensity::geom_pointdensity() +
  facet_grid(~ section) +
  scale_colour_viridis_c(option = "B") +
  coord_fixed() +
  theme_void() +
  theme(legend.position = "none")