给你 UMAP 坐标重复文章一模一样的图?

Original JunJunLab 老俊俊的生信笔记 2023-06-15

收录于合集 #单细胞笔记 21个

想回到从前

1引言

不懂的人说一两句话觉得很轻松的事情,敲敲代码而已?有可能你得花费大量时间和精力去搞明白,到最后人家还不想要,相当于白浪费一两天时间,还得自己买账,生活不易,猫猫叹气。

我还不如写篇推文分享给大家学习,多好。一篇 Nature 文章:

文章标题:

解码人肾纤维化中的肌成纤维细胞起源 (Decoding myofibroblast origins in human kidney fibrosis)

文章代码:

虽然文章给了代码,但是去看每一步做了什么还是需要大量软件知识和 R 语言知识的,此外这篇文章也不是使用我们常见的 seurat 标准流程走的,文章里的很多软件我都没没用过。还需要学习啊!

索性直接去拿文章给的处理好的数据画个图算了。

2了解

文章给了处理好的矩阵和 metadata 信息,但是是没有 barcode 的,所以需要自己起个名字,此外作者还给了 umap 的坐标信息,只有两列,并没有匹配的细胞信息,所以 默认给的细胞顺序和矩阵的是一样的。

矩阵及 metadta 信息:

https://zenodo.org/record/4059315#.YrauWEZBy3A

我们只拿 CD10negetive 数据测试。

3构建 seurat 对象

我们拿矩阵及 metadta 信息构建 seurat 对象:

library(Seurat)
library(ggplot2)
library(tidyverse)
library(Matrix)
library(ggsci)

# load gene
ginfo <- read.csv('CD10negative/kidneyMap_UMI_counts_rowData.txt')

# load mtx
mtx <- readMM('CD10negative/kidneyMap_UMI_counts.mtx')

# add gene and cell
colnames(mtx) <- paste('cell',1:ncol(mtx),sep = '')
rownames(mtx) <- ginfo$Gene.Symbol

dim(mtx)
# [1] 23446 51849

head(mtx[1:3,1:3])
# 3 x 3 sparse Matrix of class "dgTMatrix"
#             cell1 cell2 cell3
# WASH7P          .     .     .
# MIR1302-2HG     .     .     .
# AL627309.1      .     .     .

# load mata
mata <- read.csv('CD10negative/kidneyMap_UMI_counts_colData.txt')

# create seurat object
CD10neg <- CreateSeuratObject(counts = mtx,
                              meta.data = mata,
                              project = 'CD10negtive')
CD10neg
# An object of class Seurat
# 23446 features across 51849 samples within 1 assay
# Active assay: RNA (23446 features, 0 variable features)

4UMAP 亚群

我们拿 metadata 和作者提供的 umap 坐标信息可以绘制原文中一样的图:

# load pc
umap_pc <- read.csv('human_CD10negative_umapCoords.csv',header = F)
colnames(umap_pc) <- c('pc1','pc2')

# combine
meta_mer <- cbind(mata,umap_pc)

rownames(meta_mer) <- colnames(mtx)

# plot
ggplot(meta_mer,aes(x = pc1,y = pc2,color = Annotation.Level.1)) +
  geom_point(size = 0.2) +
  theme_bw(base_size = 16) +
  theme(panel.grid = element_blank(),
        aspect.ratio = 0.9) +
  scale_color_aaas(name = '') +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  xlab('UMAP 1') + ylab('UMAP 2')

5计算 normalize 数据

# standard processing
CD10neg <- NormalizeData(CD10neg, normalization.method = "LogNormalize", scale.factor = 10000)
CD10neg <- FindVariableFeatures(CD10neg, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(CD10neg)
CD10neg <- ScaleData(CD10neg, features = all.genes)

你也可以拿这个对象 重新降维聚类画图 等等:

# CD10neg <- RunPCA(CD10neg, features = VariableFeatures(object = CD10neg))
# ElbowPlot(CD10neg)
# CD10neg <- FindNeighbors(CD10neg, dims = 1:30)
# CD10neg <- FindClusters(CD10neg, resolution = 0.5)
#
# CD10neg <- RunUMAP(CD10neg, dims = 1:30)
# CD10neg <- RunTSNE(CD10neg, dims = 1:30)

# DimPlot(CD10neg, reduction = "umap")
# DimPlot(CD10neg, reduction = "tsne")
#
# DimPlot(CD10neg, reduction = "umap",group.by = "Patient.ID")
#
# DimPlot(CD10neg, reduction = "umap",group.by = "Annotation.Level.1")
# DimPlot(CD10neg, reduction = "tsne",group.by = "Annotation.Level.1")

# gene
# FeaturePlot(CD10neg,
#             features = c('ACTA2','DCN','PDGFRB','EPCAM','PTPRC','EGFL7','PLP1'),
#             ncol = 3)

6基因表达

library(reshape2)

gene <- c('ACTA2','DCN','PDGFRB','EPCAM','PTPRC','EGFL7','PLP1')

# featch data
exp <- FetchData(CD10neg,vars = gene)

df <- cbind(umap_pc,exp) %>%
  melt(.,id.vars = c('pc1','pc2'),value.name = 'expr',variable.name = 'gene')

# plot
ggplot(df,aes(x = pc1,y = pc2,color = expr)) +
  geom_point(size = 0.5) +
  # scale_color_gradient2(low = 'grey90',mid = 'green',high = 'red',
  #                       midpoint = 2.5) +
  scale_color_gradient(low = 'grey90',high = 'red') +
  theme_bw(base_size = 14) +
  theme(panel.grid = element_blank(),
        strip.background = element_rect(color = NA,fill = NA),
        panel.spacing = unit(0,'mm'),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        aspect.ratio = 1) +
  facet_wrap(~gene,nrow = 2) +
  xlab('UMAP 1') + ylab('UMAP 2')