【原】PySCENIC（三）：pyscenic單細(xì)胞轉(zhuǎn)錄因子分析可視化

TS的美夢(mèng) 2022-11-04 發(fā)布于重慶

展開(kāi)全文

pyscenic分析結(jié)束了，其實(shí)后期最主要的就是可視化了，與其說(shuō)是可視化，其實(shí)核心是對(duì)結(jié)果的提取，因?yàn)槲覀冏詈笾坏玫搅艘粋€(gè)文件：sce_SCENIC.loom,后期我們會(huì)出兩個(gè)三個(gè)文章，用來(lái)可視化結(jié)果，或者一些比較好了可應(yīng)用的思路！本節(jié)詳細(xì)注釋代碼已上傳QQ群文件！

往期內(nèi)容：

PySCENIC（一）：python版單細(xì)胞轉(zhuǎn)錄組轉(zhuǎn)錄因子分析

PySCENIC（二）：pyscenic單細(xì)胞轉(zhuǎn)錄組轉(zhuǎn)錄因子分析

覺(jué)得小編內(nèi)容有用的、有意思的，點(diǎn)贊、分享、關(guān)注一下唄！

?
1、《KS科研分享與服務(wù)》公眾號(hào)有QQ交流群，但是進(jìn)入門檻是20元，請(qǐng)考慮清楚。群里有推文的注釋代碼和示例數(shù)據(jù)，付費(fèi)內(nèi)容半價(jià)，還可以與大家交流。
2、單細(xì)胞轉(zhuǎn)錄組全流程代碼需收費(fèi)，收費(fèi)代碼包含公眾號(hào)付費(fèi)內(nèi)容，也有很多新增加的內(nèi)容。需進(jìn)群或者需單細(xì)胞代碼的小伙伴請(qǐng)?zhí)砑幼髡呶⑿帕私猓?qǐng)備注目的，除此之外請(qǐng)勿添加，謝謝！
3、付費(fèi)文章集合有打包價(jià)哦！
詳情請(qǐng)聯(lián)系作者：
?

先加載需要的R包，都加載了，沒(méi)毛病。

setwd("/home/shpc_100828/Pyscenic/")#加載分析包library(SCopesce_SCENICR)library(AUCell)library(SCENIC)#可視化相關(guān)包，多加載點(diǎn)沒(méi)毛病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_sce_SCENIC("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)

第一個(gè)可視化：

RSS分析，查看細(xì)胞類型特異性轉(zhuǎn)錄因子，需要先加載seurat對(duì)象，提取metadata信息，并進(jìn)行分析！默認(rèn)是點(diǎ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')######計(jì)算細(xì)胞特異性TFcellTypes <-  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

我們也可以提取數(shù)據(jù)，用熱圖的方式呈現(xiàn)，這里我是用ggheatmap做的，也可以用pheatmap、complexheatmap或ggplot2做。

rss_data <- rssPlot$plot$datadevtools::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

第二個(gè)可視化：

將轉(zhuǎn)錄因子分析結(jié)果與seurat對(duì)象結(jié)合，可視化類似于seurat！

next_regulonAUC <- regulonAUC[,match(colnames(human_data),colnames(regulonAUC))]dim(next_regulonAUC)
regulon_AUC <- regulonAUC@NAMEShuman_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))

上面我們展示的是轉(zhuǎn)錄因子在不同細(xì)胞的評(píng)分，按照這個(gè)道理，我們依然可以選定某種細(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)

也可以是是R語(yǔ)言中的SCENIC的展示一樣，用UMAP！

FeaturePlot(human_data, features ="FOXP3(+)")

第三個(gè)可視化：

展示轉(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