更多精彩內(nèi)容請至我的公眾號---KS科研分享與服務(wù)
先加載需要的R包赁项，都加載了葛躏，沒毛病。
setwd("/home/shpc_100828/Pyscenic/")
#加載分析包
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
#可視化相關(guān)包悠菜，多加載點沒毛病
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
library(grid)
library(ComplexHeatmap)
library(data.table)
library(ggplot2)
library(pheatmap)
將loom文件讀入R舰攒，提取數(shù)據(jù)。
sce_SCENIC <- open_loom("sce_SCENIC.sce_SCENIC")
# exprMat <- get_dgem(sce_SCENIC)#從sce_SCENIC文件提取表達(dá)矩陣
# exprMat_log <- log2(exprMat+1) # log處理
regulons_incidMat <- get_regulons(sce_SCENIC, column.attr.name="Regulons")
regulons <- regulonsToGeneLists(regulons_incidMat)
class(regulons)
regulonAUC <- get_regulons_AUC(sce_SCENIC, column.attr.name='RegulonsAUC')
regulonAucThresholds <- get_regulon_thresholds(sce_SCENIC)
第一個可視化：
RSS分析悔醋，查看細(xì)胞類型特異性轉(zhuǎn)錄因子摩窃，需要先加載seurat對象，提取metadata信息芬骄，并進(jìn)行分析猾愿！默認(rèn)是點圖！
human_data <- readRDS("~/Pyscenic/human_data.rds")
cellinfo <- human_data@meta.data[,c('celltype','group',"nFeature_RNA","nCount_RNA")]#細(xì)胞meta信息
colnames(cellinfo)=c('celltype', 'group','nGene' ,'nUMI')
######計算細(xì)胞特異性TF
cellTypes <- as.data.frame(subset(cellinfo,select = 'celltype'))
selectedResolution <- "celltype"
sub_regulonAUC <- regulonAUC
rss <- calcRSS(AUC=getAUC(sub_regulonAUC),
cellAnnotation=cellTypes[colnames(sub_regulonAUC),
selectedResolution])
rss=na.omit(rss)
rssPlot <-
plotRSS(
rss,
zThreshold = 3,
cluster_columns = FALSE,
order_rows = TRUE,
thr=0.1,
varName = "cellType",
col.low = '#330066',
col.mid = '#66CC66',
col.high = '#FFCC33')
rssPlot
image.png
我們也可以提取數(shù)據(jù)账阻，用熱圖的方式呈現(xiàn)蒂秘，這里我是用ggheatmap做的，也可以用pheatmap淘太、complexheatmap或ggplot2做姻僧。
rss_data <- rssPlot$plot$data
devtools::install_github("XiaoLuo-boy/ggheatmap")
library(ggheatmap)
library(reshape2)
rss_data<-dcast(rss_data,
Topic~rss_data$cellType,
value.var = 'Z')
rownames(rss_data) <- rss_data[,1]
rss_data <- rss_data[,-1]
colnames(rss_data)
col_ann <- data.frame(group= c(rep("Neutrophil",1),
rep("Macrophage",1),
rep("mDC",1),
rep("T cell",1),
rep("Mast",1)))#列注釋
rownames(col_ann) <- colnames(rss_data)
groupcol <- c("#D9534F", "#96CEB4", "#CBE86B", "#EDE574", "#0099CC")
names(groupcol) <- c("Neutrophil","Macrophage","mDC", "T cell","Mast")
col <- list(group=groupcol)
text_columns <- sample(colnames(rss_data),0)#不顯示列名
p <- ggheatmap(rss_data,color=colorRampPalette(c('#1A5592','white',"#B83D3D"))(100),
cluster_rows = T,cluster_cols = F,scale = "row",
annotation_cols = col_ann,
annotation_color = col,
legendName="Relative value",
text_show_cols = text_columns)
p
image.png
第二個可視化：
將轉(zhuǎn)錄因子分析結(jié)果與seurat對象結(jié)合，可視化類似于seurat蒲牧！
next_regulonAUC <- regulonAUC[,match(colnames(human_data),colnames(regulonAUC))]
dim(next_regulonAUC)
regulon_AUC <- regulonAUC@NAMES
human_data@meta.data = cbind(human_data@meta.data ,t(assay(next_regulonAUC[regulon_AUC,])))
#自己選定感興趣的或者比較重要的轉(zhuǎn)錄因子撇贺，這里我是隨機(jī)的
TF_plot <- c("ZNF561(+)","FOXP3(+)","YY1(+)","HOXB2(+)",
"TBX21(+)","TCF12(+)","STAT2(+)","SOX21(+)",
"RBBP5(+)","NR2F6(+)","NELFE(+)","MAFG(+)")
DotPlot(human_data, features = TF_plot)+
theme_bw()+
theme(panel.grid = element_blank(),
axis.text.x=element_text(hjust =1,vjust=1, angle = 45))+
labs(x=NULL,y=NULL)+guides(size=guide_legend(order=3))
image.png
上面我們展示的是轉(zhuǎn)錄因子在不同細(xì)胞的評分，按照這個道理冰抢，我們依然可以選定某種細(xì)胞松嘶，看樣本間轉(zhuǎn)錄因子的差別！
DotPlot(human_data, features = TF_plot, group.by = 'group')+
theme_bw()+
theme(panel.grid = element_blank(),
axis.text.x=element_text(hjust =1,vjust=1, angle = 45))+
theme(legend.direction = "horizontal",
legend.position = "bottom")+
labs(x=NULL,y=NULL)
image.png
第三個可視化：
展示轉(zhuǎn)錄因子平均活性挎扰！
cellsPerGroup <- split(rownames(cellTypes),
cellTypes[,selectedResolution])
regulonActivity_byGroup <- sapply(cellsPerGroup,
function(cells)
rowMeans(getAUC(sub_regulonAUC)[,cells]))
regulonActivity_byGroup_Scaled <- t(scale(t(regulonActivity_byGroup),
center = T, scale=T))
regulonActivity_byGroup_Scaled=na.omit(regulonActivity_byGroup_Scaled)
hm <- draw(ComplexHeatmap::Heatmap(regulonActivity_byGroup_Scaled, name="Regulon activity",
row_names_gp=grid::gpar(fontsize=6),
show_row_names = F))
hm #可視化所有的TF
image.png
當(dāng)然了翠订，全部展示沒有啥意義缓升，還是可以提取數(shù)據(jù)，可視化需要的TF蕴轨！
