写在开头
前面我们已经完成了两次内容的分享，今天今天写第三部分。
3.PROGENy--转录组通路活性评分
4.DoRothEA--转录组转录因子评分
正文部分
通常我们都是利用基因集合去做单细胞评分，那么是否可以通过有权重的基因集来计算样本的通路活性？这就是PROGENy开发出来的意义。之前我们分析了PROGENy在单细胞里的应用，今天我们看看在bulk-seq转录组中，如何使用PROGENy。
加载r包和示例数据
rm(list=ls())library(airway)library(DESeq2)data(airway)library(data.table)library(airway,quietly = T)data(airway)
获取表达矩阵和分组信息
#1首先我们获取airway这个包里面的airway表达量矩阵，它有8个样品，分成了两组，我们命名为'control','case' ----mat <- assay(airway) ;mat[1:4,1:8];colnames(mat)mat=mat[,c('SRR1039513',"SRR1039517","SRR1039521",'SRR1039509','SRR1039508', 'SRR1039520', 'SRR1039512', 'SRR1039516')]keep_feature <- rowSums (mat > 1) > 1ensembl_matrix <- mat[keep_feature, ]library(AnnoProbe)ids=annoGene(rownames(ensembl_matrix),'ENSEMBL','human')ids=ids[!duplicated(ids$SYMBOL),]ids=ids[!duplicated(ids$ENSEMBL),]symbol_matrix= ensembl_matrix[match(ids$ENSEMBL,rownames(ensembl_matrix)),]rownames(symbol_matrix) = ids$SYMBOLdata.frame(row.names = colnames(airway),grouplist=as.character(airway@colData$dex))group_list = as.character(airway@colData$dex);group_listgroup_list=ifelse(group_list=='untrt','control','case' )group_list = factor(group_list,levels = c('control','case' ))#其实可以使用下述的对应关系进行基因名称转换，防止AnnoProbe的滞后性head( as.data.frame(airway@rowRanges@elementMetadata@listData) )group_listsymbol_matrix[1:4,1:6]
表达矩阵，行为基因，列为样本
DEseq2差异分析（可选）
#2如果做差异分析的话------if (F) {(colData <- data.frame(row.names=colnames(symbol_matrix),group_list=group_list) )dds <- DESeqDataSetFromMatrix(countData = symbol_matrix,colData = colData,design = ~ group_list)dds <- DESeq(dds)res <- results(dds,contrast=c("group_list",levels(group_list)[2],levels(group_list)[1]))resOrdered <- res[order(res$padj),]DEG =as.data.frame(resOrdered)DEG_deseq2 = na.omit(DEG)head(DEG_deseq2)}
差异分析结果：
progeny通路分析
#3progeny通路分析--library(progeny)pathways <- progeny(symbol_matrix,scale = TRUE, organism = "Human", top = 100,perm = 1, verbose = FALSE)pathways;colnames(airway)
获取的通路活性矩阵。注意：矩阵的行为样本，列为通路：
得到通路得分矩阵之后，后续的可视化就各显神通了！
通路可视化
library(pheatmap)myColor = colorRampPalette(c("Darkblue", "white","red"))(100)pheatmap(pathways,fontsize=14, show_rownames = F,color=myColor, main = "PROGENy", angle_col = 45, treeheight_col = 0,border_color = NA)
下面热图的行为样本，不符合常规，所以需要转置下矩阵
转置矩阵之后，再可视化
pheatmap(t(pathways),fontsize=14, show_rownames = T,scale = 'row',annotation_col = data.frame(row.names = colnames(airway),grouplist=as.character(airway@colData$dex)) ,color=myColor, main = "PROGENy", angle_col = 45, treeheight_col = 0,border_color = NA)
这个热图横轴就是PROGENy分析的14个与肿瘤发生发展的信号了，纵轴是各个样本。在实际应用中，只需要提供自己数据的基因表达矩阵即可。
加上组别信息
pheatmap(t(pathways)[,c('SRR1039513',"SRR1039517","SRR1039521",'SRR1039509','SRR1039508', 'SRR1039520', 'SRR1039512', 'SRR1039516')],fontsize=14, show_rownames = T,scale = 'row', cluster_rows = T,cluster_cols = FALSE,annotation_col = data.frame(row.names = colnames(airway),grouplist=as.character(airway@colData$dex)) ,color=myColor, main = "PROGENy", angle_col = 45, treeheight_col = 0,border_color = NA)
结尾
从我的示例数据的热图结果上看。progeny通路分析对此数据集的效果不是很好，control 和treat组别之间的区别不是非常明显。
但是！是否可以把control 和treated组别的值取个平均值，再进行可视化，效果是否会好些呢？这个问题就留给读者吧。
参考：https://saezlab.github.io/progeny/articles/ProgenySingleCell.htmlhttps://www.jianshu.com/p/4058050d546ehttps://saezlab.github.io/progeny/articles/progeny.html
