從GO數(shù)據(jù)庫找基因集合做AUCell

AUCell

1.提取指定GOterm的基因

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7606976/

在這篇文章里提到，Reactive oxygen species (ROS) 活性氧相關(guān)的17個(gè)GO條目,可以把這些基因提取出來。使用msigdbr包搞定。Data_Sheet_2.XLSX是文章的附件

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7606976/bin/Data_Sheet_2.XLSX

rm(list = ls())
library(tidyverse)
dat = rio::import("Data_Sheet_2.XLSX")
g = dat[-1,1]
g

##  [1] "GO_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"                       
##  [2] "GO_REACTIVE_OXYGEN_SPECIES_METABOLIC_PROCESS"                          
##  [3] "GO_POSITIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"
##  [4] "GO_NEGATIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"
##  [5] "GO_NEGATIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_METABOLIC_PROCESS"   
##  [6] "GO_POSITIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_METABOLIC_PROCESS"   
##  [7] "GO_MITOCHONDRIAL_ELECTRON_TRANSPORT_CYTOCHROME_C_TO_OXYGEN"            
##  [8] "GO_SUPEROXIDE_GENERATING_NADPH_OXIDASE_ACTIVITY"                       
##  [9] "GO_HYDROGEN_PEROXIDE_METABOLIC_PROCESS"                                
## [10] "GO_HYDROGEN_PEROXIDE_CATABOLIC_PROCESS"                                
## [11] "GO_PENTOSE_METABOLIC_PROCESS"                                          
## [12] "GO_RESPONSE_TO_REACTIVE_OXYGEN_SPECIES"                                
## [13] "GO_REGULATION_OF_RESPONSE_TO_REACTIVE_OXYGEN_SPECIES"                  
## [14] "GO_CELLULAR_RESPONSE_TO_REACTIVE_OXYGEN_SPECIES"                       
## [15] "GO_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"         
## [16] "GO_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_METABOLIC_PROCESS"            
## [17] "GO_NEGATIVE_REGULATION_OF_RESPONSE_TO_REACTIVE_OXYGEN_SPECIES"

17條term??

2.從msigdbr查找這些term

library(msigdbr)
GO_df = msigdbr(species = "Homo sapiens",category = "C5") %>% 
  dplyr::select(gene_symbol,gs_name,gs_subcat)

dim(GO_df)

## [1] 1424471       3

table(GO_df$gs_subcat)

## 
##  GO:BP  GO:CC  GO:MF    HPO 
## 721379 115769 122674 464649

GO_df = GO_df[GO_df$gs_subcat!="HPO",]

table(g %in% GO_df$gs_name)

## 
## FALSE 
##    17

全部匹配失敗啦？

其實(shí)是因?yàn)榍熬Y不同

head(GO_df$gs_name,3)

## [1] "GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS"
## [2] "GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS"
## [3] "GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS"

head(g,3)

## [1] "GO_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"                       
## [2] "GO_REACTIVE_OXYGEN_SPECIES_METABOLIC_PROCESS"                          
## [3] "GO_POSITIVE_REGULATION_OF_REACTIVE_OXYGEN_SPECIES_BIOSYNTHETIC_PROCESS"

把前綴去掉來匹配，順便把大寫變成小寫，變得易讀一點(diǎn)。

sn = GO_df$gs_name %>%
  str_to_lower()%>%
  str_replace("_","http:///") %>%
  str_split("http:///",simplify = T)%>%
  .[,2] %>%
  str_replace_all("_"," ")
head(sn)

## [1] "10 formyltetrahydrofolate metabolic process"
## [2] "10 formyltetrahydrofolate metabolic process"
## [3] "10 formyltetrahydrofolate metabolic process"
## [4] "10 formyltetrahydrofolate metabolic process"
## [5] "10 formyltetrahydrofolate metabolic process"
## [6] "10 formyltetrahydrofolate metabolic process"

GO_df = mutate(GO_df,short_name = sn)

g = str_to_lower(g)%>%
  str_replace("_","http:///") %>%
  str_split("http:///",simplify = T)%>%
  .[,2] %>%
  str_replace_all("_"," ")
head(g)

## [1] "reactive oxygen species biosynthetic process"                       
## [2] "reactive oxygen species metabolic process"                          
## [3] "positive regulation of reactive oxygen species biosynthetic process"
## [4] "negative regulation of reactive oxygen species biosynthetic process"
## [5] "negative regulation of reactive oxygen species metabolic process"   
## [6] "positive regulation of reactive oxygen species metabolic process"

table(g %in% GO_df$short_name)

## 
## FALSE  TRUE 
##     1    16

經(jīng)過數(shù)據(jù)框搜索，發(fā)現(xiàn)匹配不上的那個(gè)其實(shí)是因?yàn)樯倭藗z空格。。。

所以直接改一下就好。

“superoxide generating nadph oxidase activity”

“superoxide generating nad p h oxidase activity”

g[!(g %in% GO_df$short_name)] = "superoxide generating nad p h oxidase activity"
table(g %in% GO_df$short_name)

## 
## TRUE 
##   17

rosgene = GO_df %>%
  filter(short_name %in% g) %>%
  pull(gene_symbol) %>%
  unique()
length(rosgene)

## [1] 402

2.AUCell

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9897923/
用這個(gè)文章里的方法，將單細(xì)胞亞群的marker基因與ros相關(guān)基因取交集，用作AUCell的基因集

The intersection of marker genes was selected based on strong population specificity (adj_p < 0.05 & |avg_log2FoldChange| > 1.5 & pct.1 > 0.5 & pct.2 < 0.5) from each cell subgroup and factors related to OS responses.

sce.all是seurat+singleR得出的對(duì)象，markers是findallmarkers得到的數(shù)據(jù)框。

load("sce.all.Rdata")
load("makers.Rdata")
k = markers$p_val_adj < 0.05 & 
  abs(markers$avg_log2FC) > 1.5 & 
  markers$pct.1 > 0.5 & 
  markers$pct.2 < 0.5
table(k)

## k
## FALSE  TRUE 
## 13100  1483

gi = intersect(markers[k,"gene"],rosgene) 
gi = factor(gi,levels = gi)
library(Seurat)
DotPlot(sce.all,features = gi,cols = "RdYlBu",
        group.by = "RNA_snn_res.0.5")+RotatedAxis()

AUCell用于計(jì)算每個(gè)細(xì)胞中特定基因集的活性程度。用上面的gi作為基因集來計(jì)算。

AUCell的三個(gè)步驟:

Build the rankings：矩陣中的每個(gè)細(xì)胞里，給基因進(jìn)行排序。 Calculate the Area Under the Curve (AUC)：計(jì)算每個(gè)細(xì)胞的AUC值 Set the assignment thresholds：計(jì)算活性區(qū)分的閾值

library(GSEABase)
geneSets <- GeneSet(as.character(gi), setName="ROS")
geneSets

## setName: ROS 
## geneIds: ACP5, SOD3, ..., PMAIP1 (total: 32)
## geneIdType: Null
## collectionType: Null 
## details: use 'details(object)'

library(AUCell)
exprMatrix = sce.all@assays$RNA@data
cells_rankings <- AUCell_buildRankings(exprMatrix)

##  min   1%   5%  10%  50% 100% 
##  474  716  794  859 1409 5996

cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings)
set.seed(333)
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist=TRUE, assign=TRUE) 

auc_thr = cells_assignment$ROS$aucThr$selected
auc_thr

##       L_k2 
## 0.07917315

3.FeaturePlot可視化

標(biāo)準(zhǔn)版：

library(Seurat)
a = as.numeric(getAUC(cells_AUC))
sce.all$auc_score = a
FeaturePlot(sce.all,features = "auc_score")

定制版：

dat<- data.frame(sce.all@meta.data, 
                sce.all@reductions$umap@cell.embeddings,
                sce.all@active.ident)
colnames(dat)[ncol(dat)] = "seurat_annotation"
class_avg <- dat %>%
  group_by(seurat_annotation) %>%
  summarise(
    UMAP_1 = median(UMAP_1),
    UMAP_2 = median(UMAP_2)
  )

library(ggpubr)
ggplot(dat, aes(UMAP_1, UMAP_2))  +
  geom_point(aes(colour  = auc_score)) + viridis::scale_color_viridis(option="A") +
  ggrepel::geom_label_repel(aes(label = seurat_annotation),
                            data = class_avg,
                            label.size = 0,
                            segment.color = NA)+
  theme(legend.position = "none") + theme_bw()

dat$auc_group = ifelse(dat$auc_score>auc_thr,"active","non_active")

4.活躍和不活躍組的細(xì)胞比例條形圖

ggplot(dat,aes(auc_group,fill = seurat_annotation))+
  geom_bar(position = "fill", alpha = 0.9)+
  scale_fill_brewer(palette = "Set1")+
  theme_classic()

http://www.itdecent.cn/p/2fb20f44da67

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容