SARS3(monocle)

# ## 安装所需要的R包
# chooseCRANmirror(graphics=FALSE)
# install.packages('Seurat')
# install.packages("tidyverse")
# install.packages("ggpubr")
# install.packages("cowplot")
# install.packages("dplyr")
# install.packages("ggplot2")
# install.packages("tidyr")
# install.packages("msigdbr")
# 
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("piano")
# BiocManager::install("HSMMSingleCell")
# BiocManager::install("monocle")


## 载入R包
library(Seurat)
library(tidyverse)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(tidyr)
library(piano)
library(msigdbr)
library(monocle)
#设置镜像
#options(repos="http://mirrors.tuna.tsinghua.edu.cn/CRAN/")
#options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")

##载入巨噬细胞
setwd("~/Rdata/sars_scRNA/")
MP <- readRDS("GSE171828_Seurat_object_monocyte_macrophage_DC.Rds")
## define cell cluster name 
MP$Annotation <- as.character(MP$Annotation)
MP$Annotation_2 <- as.character(MP$Annotation)

MP$Annotation[MP$Annotation_2 == "APOE pos FABP4 high tissue M2"] <- "APOE+ tissue macrophage"
MP$Annotation[MP$Annotation_2 == "SPP1 high fibrogenic M2"] <- "SPP1hi CHIT1int profibrogenic M2"
MP$Annotation[MP$Annotation_2 == "Transitional M1"] <- "Weakly activated M1"
MP$Annotation[MP$Annotation_2 == "Interferon stimulated M1"] <- "Highly activated M1"
MP$Annotation[MP$Annotation_2 == "Infiltrating macrophage"] <- "Monocyte-derived infiltrating macrophage"
#MP$Annotation[MP$Annotation_2 == "APOE pos FABP4 high tissue M2"] <- "APOE+ tissue macrophage"

## Fig5 b

#伪时序分析计算

SO.traj1 <- MP[,MP$Annotation %in% c("Monocyte-derived infiltrating macrophage", "Weakly activated M1", "Highly activated M1")]
pd <- new('AnnotatedDataFrame', data = SO.traj1@meta.data)
fd <- new('AnnotatedDataFrame', data = data.frame(gene_short_name = row.names(SO.traj1), row.names = row.names(SO.traj1)))
data <- as(as.matrix(SO.traj1@assays$SCT@data), 'sparseMatrix')
cds <- newCellDataSet(data, phenoData = pd, featureData = fd, expressionFamily = negbinomial.size())
#cds <- newCellDataSet(as(SO.traj1@assays$SCT@data, "matrix"), phenoData = pd, featureData = fd, expressionFamily = negbinomial.size())
cds$clusters <- SO.traj1$Annotation

#估计规模因素和离散度
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

#过滤低质量的细胞
set.seed(123)
cds <- detectGenes(cds, min_expr = 0.1)
markers <- FindAllMarkers(object = SO.traj1, min.pct = 0.25, thresh.use = 0.25)
markers <- subset(markers, p_val_adj < 0.05)
order.genes <- unique(as.character(markers$gene))

#构建单细胞轨迹
cds <- setOrderingFilter(cds, order.genes)
cds <- reduceDimension(cds = cds, max_components = 3,method = 'DDRTree')
cds <- orderCells(cds)
pdf("Fig5d.MP_M1_pseudotime_trajectory.pdf", 400/72,200/72)
plot_cell_trajectory(cds, color_by = "clusters", show_branch_points = F, cell_size = 0.5) + scale_x_reverse()
dev.off()

#构建单细胞UMAP轨迹
head(MP@reductions$umap@cell.embeddings)
head(cds@phenoData@data)
> head(MP@reductions$umap@cell.embeddings)
                             UMAP_1     UMAP_2
Ctrl-2:AAACCCACAATAGTAGx -4.5643387   4.824413
Ctrl-2:AAACCCACACCAGCTGx -5.9036785   3.614603
Ctrl-2:AAACCCACATGACTGTx  3.2230615 -11.148943
Ctrl-2:AAACCCAGTAGAGATTx  0.8020717   1.683502
Ctrl-2:AAACCCAGTCATAAAGx  3.4056079   4.289331
Ctrl-2:AAACCCATCAAGGACGx -5.1570769   4.094346
> head(cds@phenoData@data)
                            orig.ident nCount_RNA nFeature_RNA Sample_name Experimental_group nCount_SCT
Ctrl-2:AAACCCAGTCATAAAGx SeuratProject   7550.153         2237      Ctrl-2                 Ct       8394
Ctrl-2:AAAGTGAAGAAGTGTTx SeuratProject   7832.050         2136      Ctrl-2                 Ct       8429
Ctrl-2:AACAAGAAGTTGAAGTx SeuratProject   2246.858          991      Ctrl-2                 Ct       7883
Ctrl-2:AACAGGGCATATCTCTx SeuratProject  12854.146         3141      Ctrl-2                 Ct       9512
Ctrl-2:AACGAAACAAGACCTTx SeuratProject   4119.166         1723      Ctrl-2                 Ct       7721
Ctrl-2:AAGAACAGTGAGTGACx SeuratProject   6715.242         1945      Ctrl-2                 Ct       8183
                         nFeature_SCT SCT_snn_res.0.3 seurat_clusters SCT_snn_res.0.5 SCT_snn_res.0.4
Ctrl-2:AAACCCAGTCATAAAGx         2225               2               9               8               4
Ctrl-2:AAAGTGAAGAAGTGTTx         2135               2               9               8               4
Ctrl-2:AACAAGAAGTTGAAGTx         1552               5               8               6               4
Ctrl-2:AACAGGGCATATCTCTx         3106               2               9               8               0
Ctrl-2:AACGAAACAAGACCTTx         1813               2               8               6               4
Ctrl-2:AAGAACAGTGAGTGACx         1945               5               8               6               4
                                                       Annotation SCT_snn_res.1 SCT_snn_res.0.6
Ctrl-2:AAACCCAGTCATAAAGx Monocyte-derived infiltrating macrophage            10               9
Ctrl-2:AAAGTGAAGAAGTGTTx Monocyte-derived infiltrating macrophage             5               9
Ctrl-2:AACAAGAAGTTGAAGTx                      Weakly activated M1             9               8
Ctrl-2:AACAGGGCATATCTCTx Monocyte-derived infiltrating macrophage            16               9
Ctrl-2:AACGAAACAAGACCTTx                      Weakly activated M1             9               8
Ctrl-2:AAGAACAGTGAGTGACx                      Weakly activated M1             9               8
                                    Annotation_2 Size_Factor                                 clusters
Ctrl-2:AAACCCAGTCATAAAGx Infiltrating macrophage   1.0953478 Monocyte-derived infiltrating macrophage
Ctrl-2:AAAGTGAAGAAGTGTTx Infiltrating macrophage   1.0700055 Monocyte-derived infiltrating macrophage
Ctrl-2:AACAAGAAGTTGAAGTx         Transitional M1   0.8270419                      Weakly activated M1
Ctrl-2:AACAGGGCATATCTCTx Infiltrating macrophage   1.4714588 Monocyte-derived infiltrating macrophage
Ctrl-2:AACGAAACAAGACCTTx         Transitional M1   0.9654185                      Weakly activated M1
Ctrl-2:AAGAACAGTGAGTGACx         Transitional M1   1.0523062                      Weakly activated M1
                         num_genes_expressed Pseudotime State
Ctrl-2:AAACCCAGTCATAAAGx                2225   28.72614     5
Ctrl-2:AAAGTGAAGAAGTGTTx                2135   27.32175     5
Ctrl-2:AACAAGAAGTTGAAGTx                1552   16.39644     5
Ctrl-2:AACAGGGCATATCTCTx                3106   31.99751     5
Ctrl-2:AACGAAACAAGACCTTx                1813   26.84087     5
Ctrl-2:AAGAACAGTGAGTGACx                1945   27.36084     5


MP@meta.data<-MP@meta.data[rownames(cds@phenoData@data),]
MP@meta.data$Pseudotime <- cds@phenoData@data$Pseudotime
head(MP@meta.data)
 head(MP@meta.data)
                            orig.ident nCount_RNA nFeature_RNA Sample_name Experimental_group nCount_SCT
Ctrl-2:AAACCCAGTCATAAAGx SeuratProject   7550.153         2237      Ctrl-2                 Ct       8394
Ctrl-2:AAAGTGAAGAAGTGTTx SeuratProject   7832.050         2136      Ctrl-2                 Ct       8429
Ctrl-2:AACAAGAAGTTGAAGTx SeuratProject   2246.858          991      Ctrl-2                 Ct       7883
Ctrl-2:AACAGGGCATATCTCTx SeuratProject  12854.146         3141      Ctrl-2                 Ct       9512
Ctrl-2:AACGAAACAAGACCTTx SeuratProject   4119.166         1723      Ctrl-2                 Ct       7721
Ctrl-2:AAGAACAGTGAGTGACx SeuratProject   6715.242         1945      Ctrl-2                 Ct       8183
                         nFeature_SCT SCT_snn_res.0.3 seurat_clusters SCT_snn_res.0.5 SCT_snn_res.0.4
Ctrl-2:AAACCCAGTCATAAAGx         2225               2               9               8               4
Ctrl-2:AAAGTGAAGAAGTGTTx         2135               2               9               8               4
Ctrl-2:AACAAGAAGTTGAAGTx         1552               5               8               6               4
Ctrl-2:AACAGGGCATATCTCTx         3106               2               9               8               0
Ctrl-2:AACGAAACAAGACCTTx         1813               2               8               6               4
Ctrl-2:AAGAACAGTGAGTGACx         1945               5               8               6               4
                                                       Annotation SCT_snn_res.1 SCT_snn_res.0.6
Ctrl-2:AAACCCAGTCATAAAGx Monocyte-derived infiltrating macrophage            10               9
Ctrl-2:AAAGTGAAGAAGTGTTx Monocyte-derived infiltrating macrophage             5               9
Ctrl-2:AACAAGAAGTTGAAGTx                      Weakly activated M1             9               8
Ctrl-2:AACAGGGCATATCTCTx Monocyte-derived infiltrating macrophage            16               9
Ctrl-2:AACGAAACAAGACCTTx                      Weakly activated M1             9               8
Ctrl-2:AAGAACAGTGAGTGACx                      Weakly activated M1             9               8
                                    Annotation_2 Pseudotime
Ctrl-2:AAACCCAGTCATAAAGx Infiltrating macrophage   28.72614
Ctrl-2:AAAGTGAAGAAGTGTTx Infiltrating macrophage   27.32175
Ctrl-2:AACAAGAAGTTGAAGTx         Transitional M1   16.39644
Ctrl-2:AACAGGGCATATCTCTx Infiltrating macrophage   31.99751
Ctrl-2:AACGAAACAAGACCTTx         Transitional M1   26.84087
Ctrl-2:AAGAACAGTGAGTGACx         Transitional M1   27.36084


mydata<- FetchData(MP,vars = c("UMAP_1","UMAP_2","Pseudotime"))

p <- ggplot(mydata,aes(x = UMAP_1,y =UMAP_2,colour = Pseudotime))+
  geom_point(size = 1)+
  scale_color_gradientn(values = seq(0,1,0.2),
                        colours = c('blue','green','yellow'))

p1 <- p + theme_bw() + theme(panel.border = element_blank(),
                             panel.grid.major = element_blank(),
                             panel.grid.minor = element_blank(),
                             axis.line = element_line(colour = "black"))
p1

## Fig 5c伪时序分析热图
order.genes <- order.genes[!grepl("ENSMPUG", order.genes)]
pdf("Fig5e.MP_M1_pseudotime_heatmap.pdf", 5/2.54*2,9/2.54*1.5)
plot_pseudotime_heatmap(cds[order.genes,], num_clusters = 4, cores = 1, show_rownames = F, return_heatmap = T)
dev.off()

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 219,589评论 6赞 508
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 93,615评论 3赞 396
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 165,933评论 0赞 356
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,976评论 1赞 295
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,999评论 6赞 393
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,775评论 1赞 307
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,474评论 3赞 420
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 39,359评论 0赞 276
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,854评论 1赞 317
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 38,007评论 3赞 338
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 40,146评论 1赞 351
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,826评论 5赞 346
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,484评论 3赞 331
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 32,029评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 33,153评论 1赞 272
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 48,420评论 3赞 373
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 45,107评论 2赞 356

SARS3(monocle)

推荐阅读更多精彩内容