多分组芯片数据的差异分析

多分组芯片数据的差异分析 by 生信菜鸟团

继上周二对二分组的芯片进行差异基因分析及富集分析后，今天我们趁热打铁，试试看多分组的芯片能不能行~

1背景知识

数据集在这里：(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE210306)

可以看到是一个很新的数据集，文章还没有出来，先了解一下基本的设计

Cells received 15 sequential fractions of 2 Gy week, allowing irradiated cell populations a period of recovery between exposures.
Non-irradiated controls were handled identically to the irradiated cells without radiation exposure.
All of the experiments were performed within 4–10 passages after the final irradiation

共有8个样品，2个细胞系，那如何进行分组呢？可以这样试试看

MCF-7 vs MCF-7RR （对照/放疗）
MDA-MB-231 vs MDA-MB-231RR （对照/放疗）
MCF-7RR vs MDA-MB-231RR （不同细胞系放疗后）

GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，搞定后只需要一定的生物学背景对数据进行合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

2代码实操

数据集下载

rm(list = ls())  ##清屏

options(stringsAsFactors = F)
getOption('timeout')
options(timeout=10000)

##加载包
library(AnnoProbe)
library(GEOquery)
#install.packages("do")
library(do)
library(stringr) 
##方法1
getGEO("GSE210306")
if(F){
  gset <- geoChina("GSE210306")
  gset
}
##方法2
#浏览器下载 https://ftp.ncbi.nlm.nih.gov/geo/series/GSE210nnn/GSE210306/matrix/GSE210306_series_matrix.txt.gz
gset <- getGEO("GSE210306",
               filename = 'GSE210306_series_matrix.txt.gz',
               AnnotGPL = FALSE, getGPL = FALSE)
gset ##查看得到的芯片数据
a=gset  #
exp=exprs(a) ##提取表达矩阵
dim(exp)
exp[1:4,1:4] #查看dat矩阵的前4行和4列
range(exp) ## 一般在20以内的基本都是log2转化后的
dat=exp
# 发现并不需要log，所以注释了下面的代码
#dat=log2(dat)
boxplot(dat[,1:4],las=2)  
library(limma)
dat=normalizeBetweenArrays(dat)
boxplot(dat[,1:4],las=2)  
pd=pData(a) ##提取临床信息
phe=pd 
phe$title

##生成分组向量
group_list<-str_remove_all(phe$title,"\\d") #\\d表示去除数字
group_list<-str_split(group_list,',',simplify = TRUE)[,1]
group_list<-str_remove_all(group_list,"- parental")
group_list<-gsub('-','',group_list)  ##组名最好不要出现符号"-"
table(group_list)
group_list=factor(group_list,levels = c("MCF","MCFRR","MDAMB","MDAMBRR"))


##探针注释方法1----------
library(data.table)
library(stringr)
# if(!require(tinyarray))install.packages("tinyarray")
# if(!require(devtools))install.packages("devtools")
# if(!require(tinyarray))devtools::install_github("xjsun1221/tinyarray",upgrade = F)
library(tinyarray)
library(tidyverse)
b=fread('GPL6244-17930.txt',data.table = F)
head(b)
ids=data.frame(
  ID=b$ID,  ##注意这里的ID要与dat的行名对应 head(dat)
  symbol=str_split(b$gene_assignment,' // ',simplify = T)[,2]
)
head(ids)
ids=ids[ids$symbol != '',]
head(ids) 
colnames(ids)=c('probe_id','symbol') 
ids$probe_id=as.character(ids$probe_id)  
##将probe_id这一列转为字符，否则在取子集时容易报错：subscript out of bounds

##探针注释方法2--------
##基因注释和整理
# gpl<-gset@annotation
# #http://www.bio-info-trainee.com/1399.html
# #hgu133plus2
# if(F){
#   if(!require(hgu133a.db))BiocManager::install("hgu133a.db")
#   library(hgu133a.db)
#   ls("package:hgu133a.db")
#   ids <- toTable(hgu133aSYMBOL)
#   head(ids)
# }else if(F){
#   getGEO(gpl,destdir = ".")
#   ids = data.table::fread(paste0(gpl,".soft"),header = T,skip = "ID",data.table = F)
#   ids = ids[,c("ID","Gene Symbol")]
#   colnames(ids) = c("probe_id","symbol")
# }else if(T){
#   ids = idmap(gpl,type = "bioc")
# }

#####取出dat矩阵中与探针对应的probe_id
ids=ids[ids$probe_id %in%  rownames(dat),]
dat[1:4,1:4]  
dat=dat[ids$probe_id,] 


##计算基因的表达量
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  

dat['ACTB',]
dat['GAPDH',]
save(gset,exp,ids,dat,group_list,phe,file = 'step1-output.Rdata')

差异基因分析

rm(list = ls())
load('step1-output.Rdata')

# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("Glimma")
library(limma)
library(Glimma)

#创建一个desigen矩阵
design<-model.matrix(~0+group_list)
colnames(design)<-gsub("group_list","",colnames(design))
#用makeContrasts函数设定组间对比
px = levels(group_list)
contr.matrix<-makeContrasts(
  paste0(px[2],"-",px[1]), 
  paste0(px[4],"-",px[3]), 
  paste0(px[4],"-",px[2]),
    levels = colnames(design))  #对照组放在levels的第一个位置

par (mfrow=c (1,2)) ##实现一页多图

##方法1---------
#从表达计数数据中删除异方差
#原文见于CSDN:使用limma、Glimma和edgeR,RNA-seq数据分析
v<-voom(dat,design,plot=TRUE)
vfit<-lmFit(v,design)
vfit<-contrasts.fit(vfit,contrasts = contr.matrix)
efit<-eBayes(vfit)
plotSA(efit,main="Final model:Mean-variance trend")

# #检查差异基因数量
dt=decideTests(efit,p.value = 0.05,lfc = 0.5)
summary(dt)

# MCFRR-MCF MDAMBRR-MDAMB MDAMBRR-MCFRR
Down          19             7           249
NotSig     23274         23295         22577
Up            14             5           481
vennDiagram(dt[,1:3],circle.col =c("turquoise","salmon","cornflowerblue"))

##方法2--------
##1:3是因为共有3次组间比较
deg = lapply(1:3, function(i){
  topTable(efit,coef = i,number = Inf)
})
names(deg) = colnames(contr.matrix)

###选择任意两个分组进行差异表达可视化
plotMD(efit,column = 1,
       status = decideTests(efit)[,3],
       main=colnames(efit)[1] )

##给deg增加上下调信息------
library(dplyr)
library(clusterProfiler)
library(org.Hs.eg.db)

head(deg[[1]]) ##随时查看数据

for(i in 1:3){
  #1.加基因列
  deg[[i]] <- mutate(deg[[i]],symbol=rownames(deg[[i]]))
  deg[[i]] <- inner_join(deg[[i]],ids,by="symbol")
  #去重
  deg[[i]] <- deg[[i]][!duplicated(deg[[i]]$symbol),]
  logFC_t=with(deg[[i]],mean(abs(logFC)))
  P.Value_t=0.05
  #3.change
  {logFC_t=0.5
    P.Value_t = 0.05
    change=ifelse(deg[[i]] $P.Value>P.Value_t,'stable', 
                  ifelse( deg[[i]]$logFC >logFC_t,'up', 
                          ifelse( deg[[i]]$logFC < -logFC_t,'down','stable') )
    )
    deg[[i]] <- mutate(deg[[i]],change)
    print(table(deg[[i]]$change))
    deg[[i]] <- mutate(deg[[i]],v = -log10(P.Value))
  }
  #4.加ENTREZID列，后面富集分析要用
  s2e <- bitr(unique(deg[[i]]$symbol), fromType = "SYMBOL",
              toType = c( "ENTREZID"),
              OrgDb = org.Hs.eg.db)
  deg[[i]] <- inner_join(deg[[i]],s2e,by=c("symbol"="SYMBOL"))
  
}

#down stable     up 
    19  23274     14
#down stable     up 
     7  23294      6
#down stable     up 
   250  22576    481

save(deg,contr.matrix,file = "step2output.Rdata")

上下调基因并不是很多，这是为什么呢？答案可能在后续的paper中，大家可以期待一下

可视化

rm(list = ls())

##加载包
library(tinyarray)
library(dplyr)

load("step1-output.Rdata")
load( "step2output.Rdata")


pca_plot <- draw_pca(exp,group_list)
volc = lapply(1:length(deg),function(k){
  draw_volcano(deg[[k]],
               lab = colnames(contr.matrix)[k],
               pkg = 2,  ##需查看参数的意义
               logFC_cutoff  = 0.5,
               pvalue_cutoff = 0.05,
               symmetry = T)
})

dev.off()
volc[[3]] ##挑选一个看看效果


cgs = list()
for(k in 1:3){
  cgs[[k]] = filter(deg[[k]],abs(logFC) > logFC & P.Value < 0.05) %>% 
    dplyr::select(probe_id)
}
cg = intersect(intersect(cgs[[1]]$probe_id,cgs[[2]]$probe_id),cgs[[3]]$probe_id)
cg = union(union(cgs[[1]]$probe_id,cgs[[2]]$probe_id),cgs[[3]]$probe_id)
cg = c(head(sort(cg),200),tail(sort(cg),200))  
h = draw_heatmap(exp[cg,],group_list,cluster_cols = F)  ##exp行名要求是探针名


library(patchwork)
plotlist = c(list(h,pca_plot),volc)
layout <- '
AAABBB
AAABBB
CCDDEE
'
wrap_plots(plotlist) + plot_layout(design = layout)+ plot_layout(guides = 'collect')

library(ggplot2)
ggsave("patch_work.pdf")

仔细观察，可以发现dat和exp这两个矩阵是有区别的，可以利用dat进行后续的go和kegg富集分析，有兴趣的同学可以接着往下试试看

这里是粗心的小编打的补丁，补上上周两分组分析中缺失的"kegg_plot_function.R"

### ---------------
###
### Create: Jianming Zeng
### Date: 2018-07-09 20:11:07
### 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-07-09  First version
###

### Update: Xiaojie Liang
### Date: 2021-07-12 23:38:07
### Email: liangxiaojiecom@163.com
### ---------------
library(ggthemes)
library(ggplot2)
kegg_plot <- function(up.data,down.data){
  dat=rbind(up.data,down.data)
  colnames(dat)
  dat$pvalue = -log10(dat$pvalue)
  dat$pvalue=dat$pvalue*dat$group 
  
  dat=dat[order(dat$pvalue,decreasing = F),]
  
  gk_plot <- ggplot(dat,aes(reorder(Description, pvalue), y=pvalue)) +
    geom_bar(aes(fill=factor((pvalue>0)+1)),stat="identity", width=0.7, position=position_dodge(0.7)) +
    coord_flip() +
    scale_fill_manual(values=c("#0072B2", "#B20072"), guide="none") +
    labs(x="", y="" ) +
    theme_pander()  +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          #axis.ticks.x = element_blank(),
          axis.line.x = element_line(size = 0.3, colour = "black"),#x轴连线
          axis.ticks.length.x = unit(-0.20, "cm"),#修改x轴刻度的高度，负号表示向上
          axis.text.x = element_text(margin = margin(t = 0.3, unit = "cm")),##线与数字不要重叠
          axis.ticks.x = element_line(colour = "black",size = 0.3) ,#修改x轴刻度的线                         
          axis.ticks.y = element_blank(),
          axis.text.y  = element_text(hjust=0),
          panel.background = element_rect(fill=NULL, colour = 'white')
    )
}

以上就是内容的全部，其实还有另外一个包可以更快捷地完成多分组的差异分析，下次给大家演示看吧

ixxmu / mp_duty

多分组芯片数据的差异分析 #2671

多分组芯片数据的差异分析 by 生信菜鸟团

1背景知识

2代码实操

数据集下载

差异基因分析

可视化