「生信技能树」三阴性乳腺癌表达矩阵探索

「生信技能树」三阴性乳腺癌表达矩阵探索 by 生信技能树

昨天分享了一个3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，目前免费在B站:

这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz
配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC
TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA
GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA
METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC

然后马上就有了一千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！

扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！

P1 GEO-TNBC-01-文献阅读-1

首先从第一篇文章开始，GSE76275

首先介绍疾病背景
生物信息学最重要的两块应用：表达量 + 突变

要点1：198个TNBC的肿瘤样本
要点2：了解TNBC的定义，其分类是依据 IHC（免疫组化染色）

TNBC这个概念在临床应用很适合，但是在肿瘤分子水平研究这个定义就太过宽泛了
免疫组织化学（IHC）测定显示ER和PR的表达 ≤ 1%，IHC测定HER2结果为0至1+，或IHC2+但荧光原位杂交（FISII）结果阴性（无扩增）

要点3：后面课程的分析结果将验证该结论

==没发现要点4？==
要点5：芯片平台在后续下游分析的探针注释中要用到
要点6：PAM50分类器
要点7：通路注释得出结论

P2 GEO-TNBC-01-文献阅读-2

第二篇文章

在第一篇文章的基础上做了生存分析

细胞系实验不是好的肿瘤学模型，但是胜在比较方便

要点3

生存分析

第三篇文章

P3 GEO-TNBC-02-数据下载及理解-01数据库导读

具体可以看生信技能树公众号推文

重复最后一篇文章的数据和图表 PMID：30175120 数据集GSE76275
视频给的去GEO数据库找文章所用数据集的过程很详细

P4 GEO-TNBC-02-数据下载及理解-02理解exprSet对象

刚开始我觉得分段很麻烦，一直要换脚本，不过后来改一小段就保存，不会改错了影响其他的大部分，还ok。
脚本分段理由

理清思路
明确每一步骤的输入与输出
方便后续单独重复某一步骤（这个确实方便，每一步开始前都是清空环境，载入变量）

f='GSE76275_eSet.Rdata'
if(!file.exists(f)){  gset <- getGEO('GSE76275', destdir=".",                 AnnotGPL = F,     ## 注释文件                 getGPL = F)       ## 平台文件  save(gset,file=f)   ## 保存到本地}load('GSE76275_eSet.Rdata')  ## 载入数据

中间插播解释了什么是R语言里的函数
要了解表达矩阵的含义，查看表达矩阵的结构

P5 GEO-TNBC-02-数据下载及理解-03了解实验设计

分组信息 pData(矩阵) -----group
ifelse

最后

P6 GEO-TNBC-02-数据下载及理解-04检查矩阵之PCA

清空环境，加载PCA需要的两个包，载入数据~~表达矩阵 + 分组
整理PCA需要的数据：t + as.data.frame()

P7 GEO-TNBC-02-数据下载及理解-05检查矩阵之Heatmap

重新清空环境，载入数据。用第一步的数据就好，上一步的数据变换加了group不需要
apply（清晰明了，赞）：1行2列（1，2）

取最显著的前1000个基因，之前学习时，尝试过下面这两种

## 这里为什么会用到tail函数呢？sort排序后从小到大，tail就是从最后面开始取，相当于从大到小# 但是有个sort不是有个参数是按从大到小么decreasing = Tcg=names(tail(sort(apply(dat,1,sd)),1000)) #apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
### 试错dat = rnorm(1000)dim(dat) <- c(20,50)apply(dat,1,sd)sort(apply(dat,1,sd))##a = (sort(apply(dat,1,sd),decreasing = T))[1:10]a## tail更简单一点，厉害b=tail(sort(apply(dat,1,sd)),10)b

取出来的数据是一样的，顺序不一样

热图：为了突显样本之间的差异，可以先把数据标准化。后续加注释等

P8 GEO-TNBC-03-主要流程之差异分析-01差异分析导读

差异分析

P9 GEO-TNBC-03-主要流程之差异分析-02 Limma差异分析

design是limma做差异分析需要的一个输入值，了解它可由group_list构建即可

P10 GEO-TNBC-03-主要流程之差异分析-03基因ID转换

基因注释 Totable + merge

P11 GEO-TNBC-03-主要流程之差异分析-04绘制火山图

火山图

标记重要基因

##1.有指定的两个基因for_label <- dat%>%   filter(symbol %in% c("TRPM3","SFRP1")) 
#%>%# 2.head(10)，如果没有指向性的基因，就可以取前十个最差异的基因标出基因名for_label = dat %>% head(10)

##3.可以挑出上调和下调的前三个差异基因x1 = head(deg[deg$change=='up',],3)x2 = head(deg[deg$change=='down',],3)
for_label = rbind(x1,x2)

P12 GEO-TNBC-03-主要流程之差异分析-05绘制热图

绘制火山图不需要表达矩阵，只要差异基因结果的表格就可以。其中包含gene symbol和logFC和P.Value
热图：需要前面的表达矩阵，可以自定义取多少个基因来绘制

P13 GEO-TNBC-03-主要流程之差异分析-06 GO和KEGG注释分析

ID转换

library(clusterProfiler)library(org.Hs.eg.db)df <- bitr(unique(a1$ensembl_id),                   fromType = "SYMBOL",                   toType = "ENSEMBL",                   OrgDb = "org.Hs.eg.db")

merge函数可以根据相同的列名将两个矩阵合并，也可以指定两个不同的列名进行合并

g_kegg = kegg_plot(up_kegg,down_kegg)print(g_kegg)

kegg_plot这个绘图函方法使用较多，Jimmy老师将它包装成一个函数使用，在代码中可以直接找到这个函数. source('function.R')
下面是function.R

### ---------------###### Create: Jianming Zeng### Date: 2018-08-10 17:07:49### Email: jmzeng1314@163.com### Blog: http://www.bio-info-trainee.com/### Forum:  http://www.biotrainee.com/thread-1376-1-1.html### CAFS/SUSTC/Eli Lilly/University of Macau### Update Log: 2018-08-10  First version###### ---------------
draw_h_v <- function(exprSet,need_DEG,n='DEseq2',group_list,logFC_cutoff){  ## we only need two columns of DEG, which are log2FoldChange and pvalue  ## heatmap    library(pheatmap)  choose_gene=head(rownames(need_DEG),50) ## 50 maybe better  choose_matrix=exprSet[choose_gene,]  choose_matrix[1:4,1:4]  choose_matrix=t(scale(t(log2(choose_matrix+1))))   ## http://www.bio-info-trainee.com/1980.html  annotation_col = data.frame( group_list=group_list  )  rownames(annotation_col)=colnames(exprSet)  pheatmap(choose_matrix,show_colnames = F,annotation_col = annotation_col,           filename = paste0(n,'_need_DEG_top50_heatmap.png'))      library(ggfortify)  df=as.data.frame(t(choose_matrix))  df$group=group_list  png(paste0(n,'_DEG_top50_pca.png'),res=120)  p=autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')+theme_bw()  print(p)  dev.off()      if(! logFC_cutoff){    logFC_cutoff <- with(need_DEG,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) )      }  # logFC_cutoff=1    need_DEG$change = as.factor(ifelse(need_DEG$pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,                                     ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')  )  this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),                      '\nThe number of up gene is ',nrow(need_DEG[need_DEG$change =='UP',]) ,                      '\nThe number of down gene is ',nrow(need_DEG[need_DEG$change =='DOWN',])  )  library(ggplot2)  g = ggplot(data=need_DEG,              aes(x=log2FoldChange, y=-log10(pvalue),                  color=change)) +    geom_point(alpha=0.4, size=1.75) +    theme_set(theme_set(theme_bw(base_size=20)))+    xlab("log2 fold change") + ylab("-log10 p-value") +    ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+    scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)  print(g)  ggsave(g,filename = paste0(n,'_volcano.png'))  dev.off()}

注释流程

首先limma分析表达矩阵得到差异分析结果
通过差异分析结果的logFC以及Pvalue挑选出上调基因和下调基因
通过与KEGG通路数据库比对，判断是否包含这些挑选出的上下调基因，判断是否在该通路存在富集现象

P14 GEO-TNBC-04-IHC-01三阴性乳腺癌定义

对象储存了各种小变量，要了解结构，
直接载入第一步得到的表达矩阵
去genecards上查找

得到了对应的别名之后

A in%in B：A是否包含于B中，为真输出TURE，为否输出FALSE

上面尝试着看了一个所需要的基因对应的探针为3个，下面查看我们所需要的4个基因，20个探针

我们要处理的每一个数据要对它每一行每一列是什么内容都了解，再做处理

热图（改变了探针名）

可以按照基因来zscore

P15 GEO-TNBC-04-IHC-02确定基因表达为阴性

P16 GEO-TNBC-05-PAM50-01heatmap背后的生物学意义

上面把TNBC最显著的三个marker：ER,PR和HER2作图，每个marker基因对应多个探针

可以看到，这3个探针表达得很差异
之前基因名对应多个探针时直接会取最大表达量的探针，这个图完美的解释了这一点：最大表达量的探针才有足够的信息来符合生物学规律在TNBC不表达，在非TNBC高表达
看不到其他那些探针差异表达并不代表没有差异，极大值和极小值的存在，导致一些变化被颜色掩盖了，衡量变化的参数应该通过计算，是否具有统计学差异。可以通过zscore探索一下每个基因的高低，再画热图，就可以看出这些低表达的基因是否有隐藏的差异

x=t(scale(t(x)))x[x>2]=2x[x< -2]= -2

P17 GEO-TNBC-05-PAM50-02PAM50

通过分类基因把样本分为了4个组

class
dim
str

可以查看一下管家基因 GAPDH 的表达

通过分类，得到了5个分组，4个肿瘤组，一个正常组

if(T){  ddata=t(dat)  ddata[1:4,1:4]  s=colnames(ddata);head(s)  library(org.Hs.eg.db)  s2g=toTable(org.Hs.egSYMBOL)  g=s2g[match(s,s2g$symbol),1];head(g)  #  probe Gene.symbol Gene.ID  dannot=data.frame(probe=s,                    "Gene.Symbol" =s,                     "EntrezGene.ID"=g)  ddata=ddata[,!is.na(dannot$EntrezGene.ID)]#去除列的NA值  dannot=dannot[!is.na(dannot$EntrezGene.ID),] #去除行的NA值

模型有很多种，但是pam50模型最出名

P18 GEO-TNBC-05-PAM50-03两个热图的区别

看看这5个分类是否有实际意义，可以画热图看看

前面用1000个差异表达的基因绘制了热图，接下来用PAM50分类器包含的基因来绘制热图。一样是从原始表达矩阵中提取出这50个基因在所有样本中的表达量来绘制热图。

library(genefu)pam50genes=pam50$centroids.map[c(1,3)]pam50genes[pam50genes$probe=='CDCA1',1]='NUF2'pam50genes[pam50genes$probe=='KNTC2',1]='NDC80'pam50genes[pam50genes$probe=='ORC6L',1]='ORC6'x=dat###判断一下50个基因是否在原来的表达矩阵中x=x[pam50genes$probe[pam50genes$probe %in% rownames(x)] ,]

归一化后再看,只看基因内部的表达量

x=t(scale(t(x)))x[x>1.6]=1.6x[x< -1.6]= -1.6

P19 GEO-TNBC-06-GSEA-01由超几何分布到GSEA

回到前面的差异分析，使用合适的阈值挑选出的基因在250-500之间都ok
enrichKEGG函数需要的输入数据为基因的entreID

可以看出这个通路很显著

下面就是GSEA

GSEA需要一个包含基因名gene symbol以及logFC值的表格。

  geneList=DEG$logFC  names(geneList)=DEG$symbol  geneList=sort(geneList,decreasing = T)

KEGG和GO的注释数据库被做成了一个可操作的R包，这里的GSEA分析需要用到的数据库需要自己下载

这里的循环是针对下载的每一个注释文件

 ## 下面使用lapply循环读取每个gmt文件，并且进行GSEA分析  gsea_results <- lapply(gmts, function(gmtfile){    # gmtfile=gmts[2]    # 如果不确定循环里面的代码，就尝试赋值，测试下面的代码。    geneset <- read.gmt(file.path(d,gmtfile))     print(paste0('Now process the ',gmtfile))    egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)

说明这个基因很显著，
下面这个就一般，不怎么显著

P20 GEO-TNBC-06-GSEA-02运行GSEA代码

最后得到了一系列基因集
就可以拿ES比较高的GSEA图

P21 GEO-TNBC-06-GSEA-03GSEA原理

去MSigDB下载基因集数据

可以下载ALL，把下面的几个都包装在一起了。也可以一个个的下载好分开走流程

这里字错了，应该是“以及”吧。

GSEA的输入输出

第一：要下载相应的GSEA基因集

第二：关于GSEA 理解输入输出

1.需要输入一个geneList选取差异表达的基因
2.需要输入gmts数据集（从GSEA官网下载）：
d='./MsigDB/symbols'gmts <- list.files(d,pattern = 'all')gmts
3.lapply函数是对gmts基因集的循环，如果下载全基因集或者只关注一个基因集，那就不需要循环：如注释信息 # gmtfile=gmts[2]:相当于只取第二个基因集
geneset <- read.gmt(file.path(d,gmtfile)) ：用read.gmt 这个函数读取gmt文档信息
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE) head(egmt)：相当于还是主函数。主要的数据分析
gseaplot(egmt, geneSetID = rownames(egmt[1,])) ：选取了一条基因集画图：选取的基因集的名字是rownames(egmt[1,])
7：gsea_results_list<- lapply(gsea_results, function(x){ cat(paste(dim(x@result)),'\n') x@result}) ##########这个的意思应该是吧gsea_results进行了操作列出 gsea_results_list
8：gsea_results_df <- do.call(rbind, gsea_results_list) ###########把gsea_results_list做转化成data.frame。do.call神奇操作
9：gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE') 。选择有差异的基因集进行画图。此处要注意gsea_results[[2]]其中的数字2 要和KEGG_CELL_CYCLE是对应的。所以画图前一定要用notepad++查看下载下来的基因集和数字的对应关系。

P22 GEO-TNBC-07-GSVA-01GSVA与GSEA的比较

还有ssGSEA，PGSEA等

GSVA也需要的是基因集啊和表达矩阵
GSEA需要的是基因的排序 + 基因集
GSVA在每个样本中计算得到一个值，可以类比于ES

之前进行差异分析的数据是==基因和样本==，现在做差异分析的数据是==通路和样本==

之前得到基因的差异分析结果，后续绘制了热图喝火山图。同样，通路的差异分析结果也可以做这些后续分析

P23 GEO-TNBC-07-GSVA-02GSVA代码及结果解释

GSVA的输入输出

计算每个样本在基因集中的富集程度，然后计算Fold change 和Pvalue。做类似于热图的图。所以后面的很多算法和制作热图一致。
关注input：函数需要输入X=dat这个expression data，row为基因名，col为sample 。另外还需要gmt这个数据集。
d='./MSigDB/symbols/'gmts=list.files(d,pattern = 'all') 通过list.files这个可以列出当前路径下的所有gmt文件
es_max <- lapply(gmts, function(gmtfile) ；使用lapply函数进行循环，对内一个下载的基因集进行操作；如果我们只关心一个基因集，就不需要这个循环了比如 #gmtfile=gmts[8]
geneset <- getGmt(file.path(d,gmtfile)) es.max <- gsva(X, geneset, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)
#######这一步很重要，是函数的最重要步骤。geneset <- getGmt(file.path(d,gmtfile)) 通过这个getGmt这个函数可以从gmt 文件中获得geneset；gsva这个主函数通过输入X这个表达矩阵和geneset数据集进行计算Fold change 和P value。
5adjPvalueCutoff <- 0.001 logFCcutoff <- log2(2)####进行阈值调整
6.帅选差异表达的基因集##类似做热图和火山图（因为实验设计是本身存在不同的分组信息）；
6.1
es_deg <- lapply(es_max, function(es.max)这个函数对所有的基因集得到的es.max进行循环操作
6.2
分组内容如下：design <- model.matrix(~0+factor(group_list)) colnames(design)=levels(factor(group_list)) rownames(design)=colnames(es.max)
6.3 deg = function(es.max,design,contrast.matrix) 这个函数对es.max按照分组信息进行基因集富集情况进行差异分析（类似于寻找差异表达基因）
7.最终会输出ex.max 和ex.deg这两个数据
8.最后作图，做差异富集的图，（类似差异表达的基因）

### 对 MigDB中的全部基因集做GSVA分析。## 还有ssGSEA, PGSEA{  load(file = 'step1-output.Rdata')  # 每次都要检测数据  dat[1:4,1:4]    library(hgu133plus2.db)  ids=toTable(hgu133plus2SYMBOL) #通过看hgu133plus2.db这个包的说明书知道提取probe_id（探针名）和symbol（基因名）的对应关系的表达矩阵的函数为toTable  head(ids)  dat=dat[ids$probe_id,]  dat[1:4,1:4]   ids$median=apply(dat,1,median)  ids=ids[order(ids$symbol,ids$median,decreasing = T),]  ids=ids[!duplicated(ids$symbol),]  dat=dat[ids$probe_id,]  rownames(dat)=ids$symbol  dat[1:4,1:4]      X=dat  table(group_list)  ## Molecular Signatures Database (MSigDb)   d='./MSigDB/symbols/'  gmts=list.files(d,pattern = 'all')  gmts  library(ggplot2)  library(clusterProfiler)  library(org.Hs.eg.db)  library(GSVA) # BiocManager::install('GSVA')    if(T){    es_max <- lapply(gmts, function(gmtfile){       #gmtfile=gmts[8];gmtfile      geneset <- getGmt(file.path(d,gmtfile))        es.max <- gsva(X, geneset,                      mx.diff=FALSE, verbose=FALSE,                      parallel.sz=1)      return(es.max)    })    adjPvalueCutoff <- 0.001    logFCcutoff <- log2(2)    es_deg <- lapply(es_max, function(es.max){      table(group_list)      dim(es.max)      design <- model.matrix(~0+factor(group_list))      colnames(design)=levels(factor(group_list))      rownames(design)=colnames(es.max)      design      library(limma)      contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),                                     levels = design)      contrast.matrix<-makeContrasts("TNBC-noTNBC",                                     levels = design)            contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较            deg = function(es.max,design,contrast.matrix){        ##step1        fit <- lmFit(es.max,design)        ##step2        fit2 <- contrasts.fit(fit, contrast.matrix)         ##这一步很重要，大家可以自行看看效果                fit2 <- eBayes(fit2)  ## default no trend !!!        ##eBayes() with trend=TRUE        ##step3        res <- decideTests(fit2, p.value=adjPvalueCutoff)        summary(res)        tempOutput = topTable(fit2, coef=1, n=Inf)        nrDEG = na.omit(tempOutput)         #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)        head(nrDEG)        return(nrDEG)      }            re = deg(es.max,design,contrast.matrix)      nrDEG=re      head(nrDEG)       return(nrDEG)    })  }     gmts    save(es_max,es_deg,file='gsva_msigdb.Rdata')    load(file='gsva_msigdb.Rdata')    library(pheatmap)  lapply(1:length(es_deg), function(i){    # i=2    print(i)    dat=es_max[[i]]    df=es_deg[[i]]    df=df[df$P.Value<0.01 & abs(df$logFC) > 0.3,]   #######在这里调阈值，保证输出    print(dim(df))    if(nrow(df)>5){    #######nrow(df)>5:只有差异富集的基因集大于五个以上才画图      n=rownames(df)      dat=dat[match(n,rownames(dat)),]  ######n=rownames(df),筛选表达矩阵数据      ac=data.frame(g=group_list)   ########对应后面的图中的annotation_col      rownames(ac)=colnames(dat)      rownames(dat)=substring(rownames(dat),1,50)      pheatmap::pheatmap(dat,                          fontsize_row = 8,height = 11,                         annotation_col = ac,show_colnames = F,                         filename = paste0('gsva_',strsplit(gmts[i],'[.]')[[1]][1],'.pdf'))          }})    adjPvalueCutoff <- 0.001  logFCcutoff <- log2(2)  df=do.call(rbind ,es_deg)  es_matrix=do.call(rbind ,es_max)  df=df[df$P.Value<0.01 & abs(df$logFC) > 0.5,]  write.csv(df,file = 'GSVA_DEG.csv')}

整个课程清晰明了，代码也都给了，只要跟着视频一步步跑下来，理解一下代码的含义，应该都能完整的得到数据图。

但是最主要的是要理解生物学背景。多爬楼看看生信技能树以及生信菜鸟团，会得到意想不到的收获。特别是早期的文章，虽然不是最新进展，但是很基础，现在的推文就有些复杂了哈哈哈，变高级了。

这些代码大家都可以测试一下，

也欢迎加入我们的其它分门别类的学习交流群，如下：

ixxmu / mp_duty

「生信技能树」三阴性乳腺癌表达矩阵探索 #321