基因集变异分析（GSVA）
软件：
- GSEAbase
- GSVA
- limma
- Seurat
首先，针对某个细胞类型继续细分其细胞亚型：
library(Seurat)library(GSEABase)library(Biobase)library(genefilter)library(limma)library(RColorBrewer)library(GSVA)library(dplyr)###读入Seurat对象Seurat.obj<-readRDS("seurat.rds")DefaultAssay(Seurat.obj) <- "integrated"#提取第三个cluster或者加入你定义了细胞类型，则为该细胞名称sub <- subset(Seurat.obj, idents = "3")##为该细胞组定义身份(stim为case-control分组依据，名称来自于官网）Idents(sub) <- "stim"##再跑一次高解析度的聚类，分出细胞亚类sub<-ScaleData(sub, verbose = FALSE)sub<-RunPCA(sub,npcs=30,verbose = FALSE)sub<-RunUMAP(sub, reduction = "pca", dims = 1:20)sub<-FindNeighbors(sub, reduction = "pca", dims = 1:20)sub <- FindClusters(sub, resolution = 0.2)###保存该对象saveRDS(sub,"sub_cluster.rds")##可以把身份加上case-control和聚类双重信息sub$celltype.case<-paste(Idents(sub),sub$stim,sep="_")sub$celltype<-Idents(sub)Idents(sub)<-"celltype.case"
接着，做 GSVA 分析，获取 GSVA 矩阵，做 GSVA 分析可以参考传统 bulk-RNA 分析的方法
这里还需要注意几个问题：
- 处理普通 RNA 数据需要预先过滤，但是单细胞数据取自 Seurat 对象，已经预先过滤好了
- 如果输入是原始 counts 值，需要设置参数
kcdf="Possion"，但如果是 TPM 值，默认就好，因为我们输入是标准化后的数据，所以用默认参数 - 默认参数
mx.diff=TURE，结果是一个类似于负二项分布，因为后面要做差异分析，所以需要使用该参数，如果设置mx.diff=FALSE，则为高斯分布
##读取Geneset#（http://software.broadinstitute.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.0/h.all.v7.0.symbols.gmt)geneSets <- getGmt("h.all.v7.0.symbols.gmt")##获取表达矩阵，这里选择标准化后但没有均一化的数据expMat<-GetAssayData(sub,slot="data")mydata<-as.matrix(expMat)##提供列注释annotation_col<-data.frame(Type=Idents(sub))##运行GSVA，返回二项式分布的结果res_es <- gsva(mydata, geneSets, min.sz=10, max.sz=500, verbose=FALSE, parallel.sz=8)
然后，使用 limma 线性模型分析新表达矩阵的差异，寻找 case-control 差异通路
##设置分组grouP<-as.factor(annotation_col$Type)desigN <- model.matrix(~ grouP + 0)rownames(desigN)<-colnames(mydata)#head(desigN)#这里需要手动构建差异矩阵，我这里用每个亚类的case和control比较，寻找其中的差异comparE <- makeContrasts(clu1=grouP0_case-grouP0_control,clu2=grouP1_case-grouP1_control,clu3=grouP2_case-grouP2_control,levels=desigN)#构建线性模型fiT <- lmFit(res_es, desigN)fiT2 <- contrasts.fit(fiT, comparE)fiT3 <- eBayes(fiT2)#这里获取一步ANNOVA检验获取的不同分组的log2FC,可以用来构造热图Diff<-topTableF(fiT3,adjust="BH",p.value=0.05,num=50)write.table(Diff,"diff_path.txt",sep="\t",quote=FALSE)
然后绘图
library(pheatmap)library(dplyr)data<-read.table("diff_path.txt",sep="\t",stringsAsFactors = FALSE,header = TRUE,check.names = FALSE)data$name<-rownames(data)data<-tbl_df(data)##提取cluster对应的log2FCdata<-select(data,name,0,1,2)data<-as.data.frame(data)nn<-gsub("HALLMARK_","",data$name)nn<-tolower(gsub("_"," ",nn))rownames(data)<-nndata<-data[,-1]pheatmap(haha,color = colorRampPalette(c("navy", "white", "firebrick3"))(100),angle_col = 45)
然后，使用 limma 构建不同分组，分别对比不同 cluster 之间的差异通路
##因为要使用cluster区别差异，所以先恢复细胞类型定义Idents(sub)<-"celltype"annotation_col<-data.frame(Type=Idents(sub))some_cluster<-annotation_colsome_cluster$Type<-as.character(some_cluster$Type)#把0-clu的定义一组，不是0-clu的定义为另外的组some_cluster$Type[which(some_cluster$Type!="0")]<-"control"some_cluster$Type[which(some_cluster$Type="0")]<-"case"#构建分组矩阵grouP<-as.factor(some_cluster$Type)desigN <- model.matrix(~ grouP + 0)rownames(desigN)<-colnames(mydata)comparE <- makeContrasts(grouPcase-grouPcontrol,levels=desigN)#构建线性模型fiT <- lmFit(res_es, desigN)fiT2 <- contrasts.fit(fiT, comparE)fiT3 <- eBayes(fiT2)#提取差异通路Diff<-topTable(fiT3,p.value=0.05,num=Inf)some_diff<-Diffsome_diff<-Diff['t']colnames(some_diff)<-"0"##这里因为有多个分组互相比较，最后需要整合不同分组的t-value表，使用mergesome_diff$name<-rownames(some_diff)fin<-merge(some_diff,fin,by="name")write.table(fin,"diff_path.txt",sep="\t",quote=FALSE,row.names=FALSE)
以上完成后，获取的差异通路信息可以通过画热图的方式展示，画图和前面一致
TF 转录因子分析
按照 SCENIC 的官方教程，跑完整个流程（这里耗时很长，所以不能直接跑，要提交作业活着 screen,避免中间间断
SCENIC 基因互作数据库需要提前下载，人类的一共两个，每个 1.2G
hg19-500bp-upstream-7species.mc9nr.feather
hg19-tss-centered-10kb-7species.mc9nr.feather
library(Seurat)library(SCENIC)#读取对象endo<-readRDS("sub_cluster.rds")DefaultAssay(endo) <- "integrated"endo$celltype.case<-paste(Idents(endo),endo$stim,sep="_")endo$celltype<-Idents(endo)Idents(endo)<-"celltype.case"##生成亚型分类信息cellInfo<-data.frame(CellType=Idents(endo))#提取表达矩阵exprMat<-GetAssayData(endo,slot="data")exprMat<-as.matrix(exprMat)##很关键，很多中间数据都保存在这里面，不要乱改名字dir.create("int")saveRDS(cellInfo, file="int/cellInfo.Rds")org="hgnc" # or hgnc, or dmel##这个数据库需要从官网下载dbDir="/SingleCell/SCENIC/cisTarget_databases"myDatasetTitle="someproject"data(defaultDbNames)dbs <- defaultDbNames[[org]]scenicOptions <- initializeScenic(org=org, dbDir=dbDir, dbs=dbs, datasetTitle=myDatasetTitle, nCores=4)scenicOptions@inputDatasetInfo$cellInfo <- "int/cellInfo.Rds"saveRDS(scenicOptions, file="int/scenicOptions.Rds")###过滤参数需要根据你输入的数据格式有所调整genesKept <- geneFiltering(exprMat, scenicOptions=scenicOptions,minCountsPerGene=3*.01*ncol(exprMat),minSamples=ncol(exprMat)*.01)exprMat_filtered <- exprMat[genesKept, ]rm(endo)runCorrelation(exprMat_filtered, scenicOptions)runGenie3(exprMat_filtered, scenicOptions)scenicOptions@settings$verbose <- TRUEscenicOptions@settings$nCores <- 4scenicOptions@settings$seed <- 123runSCENIC_1_coexNetwork2modules(scenicOptions)runSCENIC_2_createRegulons(scenicOptions)runSCENIC_3_scoreCells(scenicOptions, exprMat)
接下来，我们已经得到了新的转录调控打分矩阵，需要做的就是继续进行差异分析，我们接下来进行不同细胞亚类型之间的差异分析，从而确认到底是哪些互作专一性影响某个细胞亚类。
library(SCENIC)library(AUCell)library(Seurat)library(limma)##这里需要在/int的上级目录scenicOptions <- readRDS("int/scenicOptions.Rds")cellInfo<-readRDS("int/cellInfo.Rds")regulons <- loadInt(scenicOptions, "regulons")#获取互作基因关系regulons2 <- loadInt(scenicOptions, "aucell_regulons")#获取打分矩阵regulonAUC <- loadInt(scenicOptions, "aucell_regulonAUC")regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]endo<-readRDS("sub_cluster.rds")###差异分析,这里使用的方法和前面按照不同cluster分组的方法一致annotation_col<-data.frame(Type=Idents(endo))some_clu<-annotation_colsome_clu$Type<-as.character(some_clu$Type)some_clu$Type[which(some_clu$Type!="0")]<-"control"some_clu$Type[which(some_clu$Type=="0")]<-"case"grouP<-as.factor(some_clu$Type)desigN <- model.matrix(~ grouP + 0)rownames(desigN)<-colnames(getAUC(regulonAUC))comparE <- makeContrasts(grouPcase-grouPcontrol,levels=desigN)fiT <- lmFit(as.data.frame(getAUC(regulonAUC)), desigN)fiT2 <- contrasts.fit(fiT, comparE)fiT3 <- eBayes(fiT2)Diff<-topTable(fiT3,p.value=0.05,num=50)some_diff<-Diffsome_diff<-Diff['t']colnames(some_diff)<-"0"###这里依旧按照上面的方法采用merge把不同分组的合并一下然后保存
接下来，我们需要分析下显著差异的某个转录因子所 target 的基因所对应的功能，我们可以使用 Y 叔的clusterprofiler
library(clusterProfiler)library(org.Hs.eg.db)library(DOSE)##选择junB基因对应的通路gene<-regulons2$`JUNB (52g)`id<-bitr(gene, fromType="SYMBOL", toType="ENTREZID",OrgDb="org.Hs.eg.db")ego<-enrichGO(OrgDb="org.Hs.eg.db", gene = id$ENTREZID,pvalueCutoff = 0.05,readable=TRUE)ekk<-enrichKEGG(gene=id$ENTREZID,organism='hsa',pvalueCutoff=0.05)plotdata<-ego@resultplotdata$GeneRatio<-sapply(plotdata$GeneRatio,function(x) eval(parse(text = x)))write.table(plotdata,"JUNB_ego.txt",sep="\t",row.names=FALSE,quote=FALSE)plotdata<-ekk@resultplotdata$GeneRatio<-sapply(plotdata$GeneRatio,function(x) eval(parse(text = x)))write.table(plotdata,"JUNB_ekk.txt",sep="\t",row.names=FALSE,quote=FALSE)
这里可以绘制条形图
接下来，我们还需要把对应的转录因子以及其 targets 对应的 AUCell score 投影到 umap 图上
library(Seurat)library(ggplot2)endo<-readRDS("sub_cluster.rds")##因为要将转录调控矩阵在聚类图上显示，所以需要将表达矩阵存储到Seurat对象中#但是因为features数目不对，所以只是替换一下对应的数值#feature提取依旧使用基因的名字，按照顺序对应就是了fin<-GetAssayData(endo,slot="data")fin[1:length(rownames(getAUC(regulonAUC))),]<-getAUC(regulonAUC)new_endo<-SetAssayData(endo,slot="data",new.data=fin,assay="RNA")pdf("test3.pdf",width=4,height=4)FeaturePlot(new_endo,feature="RP4-669L17.10",cols = c("grey", "red"),min.cutoff="q9",slot="data")+theme(title = element_text(face="bold"),text = element_text(face="bold"))+theme(panel.grid.major = element_line(linetype = "blank"), legend.text = element_text(size = 10),panel.grid.minor = element_line(linetype = "blank"), panel.background = element_rect(colour = "black", size = 2, linetype = "solid"))dev.off()
若有收获，就点个赞吧
0 人点赞
