TNBC数据分析-GSE27447-GPL6244

GSM678364 non-triple negative breast cancer tumor B2
GSM678365 non-triple negative breast cancer tumor B4
GSM678366 non-triple negative breast cancer tumor B5
GSM678367 non-triple negative breast cancer tumor B8
GSM678368 non-triple negative breast cancer tumor B9
GSM678369 triple negative breast cancer tumor B12
GSM678370 non-triple negative breast cancer tumor B14
GSM678371 non-triple negative breast cancer tumor B17
GSM678372 non-triple negative breast cancer tumor B20
GSM678373 triple negative breast cancer tumor B23
...

rm(list = ls())
library(AnnoProbe)
library(GEOquery)  
library(ggplot2) 
library(ggstatsplot) 
library(reshape2)
library(patchwork)

gse_number <- 'GSE27447'
gpl_number = "GPL10558"
gset <- geoChina(gse_number)
a=gset[[1]] 
dat=exprs(a) 
dim(dat)

# 检查，判断需不需要取log
dat[1:4,1:4] 
dat = log2(dat)
dat[1:4,1:4] 
library(limma)
dat=normalizeBetweenArrays(dat)

# 画图，使用ggplot需宽数据变长数据
class(dat)
data <- as.data.frame(dat)
data <- melt(data)
head(data)
title <- paste (gse_number, "/", a@annotation, sep ="")
p1 <- ggplot(data,aes(x=variable,y=value))+
      geom_boxplot()+
      theme_ggstatsplot()+
      theme(panel.grid = element_blank(),
            axis.text=element_text(size=10,face = 'bold'),
            axis.text.x=element_text(angle=90),
            plot.title = element_text(hjust = 0.5,size =15))+
      xlab('')+
      ylab('')+
      ggtitle(title)
p1

pd = pData(a)
dim(pd)
library(stringr)
group_list=ifelse(str_detect(pd$`disease state:ch1`,"non-triple"),"control","TNBC")
group_list = factor(group_list,
               levels = c("control","TNBC"))
table(group_list)

library(stringr)
ids=idmap('GPL6244')
#超级好用的函数，首选，如果不行再尝试其他
dim(ids)

> head(ids)
           probe_id   symbol
586103 ILMN_1343291   EEF1A1
586104 ILMN_1343295    GAPDH
586105 ILMN_1651209 SLC35E2A
586106 ILMN_1651221   EFCAB1
586107 ILMN_1651228    RPS28
586108 ILMN_1651229    IPO13

library(tidyr)
library(dplyr)
#接下来，使探针与基因symbol一一对应
ids=as.data.frame(ids)
table(rownames(dat) %in% ids$probe_id)
dat=dat[rownames(dat) %in% ids$probe_id,]
ids=ids[match(rownames(dat),ids$probe_id),]
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
ids=ids[ids$probe_id %in%  rownames(dat),]
dat=dat[ids$probe_id,]

#下一步是：基因symbol去冗余-按照表达量最大值筛选
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/exp
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
#把ids的symbol这一列中的每一行给dat作为dat的行名
rownames(dat)=ids$symbol
dat[1:4,1:4]
table(group_list)
save(dat,pd,ids,group_list,gse_number,gpl_number,file = "step1-output.Rdata")

> dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
       GSM678364 GSM678365 GSM678366 GSM678367
ZZZ3    8.662448  9.804468  9.144349  9.707848
ZZEF1   8.788682  9.189457  8.981535  8.599375
ZYX    10.648260 10.034978 10.622474 10.430494
ZYG11B  8.769574  8.913443  8.958082  8.715722

>table(group_list)
group_list
control    TNBC 
     14       5  

## 下面是画PCA的必须操作，需要看说明书。
exp <- dat
exp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame 
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra") 
library(ggstatsplot)
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
dim(exp)
exp[,12488]
# 画图，主成分分析图p2
this_title <- paste0(gse_number,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group_list, # color by groups
                         palette = "Dark2",
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups")+
                   ggtitle(this_title)+
                   theme_ggstatsplot()+
                   theme(plot.title = element_text(size=15,hjust = 0.5))

p2

# 下面是1000_sd热图
library(pheatmap)
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
n=t(scale(t(dat[cg,]))) 
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
ac=data.frame(Group=group_list)
rownames(ac)=colnames(n)
# 画图，高变基因的表达量热图p3
p3 <- pheatmap::pheatmap(n,
               show_colnames =F,
               show_rownames = F,
               main = gse_number,
               annotation_col=ac,
               breaks = seq(-3,3,length.out = 100))#因为已经手动设置了表达量最大值，所以，可以不用设置break
p3

# 画图，样品相关性热图p4
colD=data.frame(Group=group_list)
exprSet=t(exp)
rownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的exp
p4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列 操作
                   annotation_col = colD,
                   show_rownames = F,
                   main = gse_number
                   )
p4

library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4) #设置全局的数字有效位数为4
deg = topTable(fit,coef=2,adjust='BH', n=Inf) 
#结果检验
boxplot(dat[rownames(deg)[1],]~group_list)
deg[1,]
#核对无误

> deg[1:10,]
         logFC AveExpr      t   P.Value adj.P.Val     B  symbol
NEK11   -1.650   6.932 -5.925 8.582e-06   0.07263 2.905   NEK11
EGR2     2.171   6.356  5.836 1.044e-05   0.07263 2.761    EGR2
HOXB3   -1.468   7.031 -5.790 1.157e-05   0.07263 2.685   HOXB3
ARFGEF3 -1.904   9.214 -5.546 1.995e-05   0.09393 2.281 ARFGEF3
NUDT16  -1.015   6.484 -5.427 2.609e-05   0.09828 2.080  NUDT16
CHRM3    2.228   5.439  5.233 4.045e-05   0.11969 1.749   CHRM3
C3orf14 -2.128   6.837 -5.192 4.448e-05   0.11969 1.677 C3orf14
ANKEF1  -1.333   7.065 -5.123 5.199e-05   0.12008 1.558  ANKEF1
TSPAN1  -4.524   7.740 -5.080 5.737e-05   0.12008 1.483  TSPAN1
GFRA3    4.226   3.723  4.940 7.916e-05   0.14912 1.237   GFRA3

nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))#简单画图看一下
df=nrDEG
df$v= -log10(P.Value) #df新增加一列'v',作为新的绘图参数，值为-log10(P.Value) 
#设定上下调基因
df$g=ifelse(df$P.Value>0.05,'stable', 
              ifelse( df$logFC >2,'up', 
                      ifelse( df$logFC < -2,'down','stable') )
)
#统计上下调基因数量
table(df$g)
#给绘制火山图用的数据新增一列symbol
df$name=rownames(df)
head(df)
logFC_t = 2
#设置可循环使用的plot标题
this_tile <- paste0('Cutoff for logFC is ',round(logFC_t,3),
                      '\nThe number of up gene is ',nrow(df[df$g == 'up',]) ,
                      '\nThe number of down gene is ',nrow(df[df$g == 'down',])
)
#画图，火山图p5
p5 <- ggplot(data = df, 
              aes(x = logFC, 
                  y = -log10(P.Value))) +
    geom_point(alpha=0.6, size=1.5, 
               aes(color=g)) +
    ylab("-log10(Pvalue)")+
    scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+
    geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) +
    xlim(-3, 3)+
    theme_classic()+
    ggtitle(this_tile )+
    theme(plot.title = element_text(size=12,hjust = 0.5),
          legend.title = element_blank(),
          )
p5

#热图
library(pheatmap)
x=deg$logFC 
names(x)=rownames(deg) 
cg=c(names(head(sort(x),100)),
       names(tail(sort(x),100)))#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n) 
# 画图，上下调的差异基因热图p6
p6 <- pheatmap(n,show_colnames =F,
           show_rownames = F,
           cluster_cols = T, 
           main = gse_number,
           annotation_col=ac) 
p6

# 加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
deg$symbol=rownames(deg)
library(org.Hs.eg.db)
library(clusterProfiler)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人
#bitr()用于SYMBOL转ENTREZID
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
dim(deg)
dim(s2e)
setdiff(deg$symbol,s2e$SYMBOL)
DEG <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))

library(dplyr)
library(ggplot2)  
geneList=DEG$logFC
names(geneList)=DEG$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)
library(clusterProfiler)
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',#按需替换
                  #nPerm        = 1000,
                  minGSSize    = 10,
                  pvalueCutoff = 0.9,
                  verbose      = FALSE)
tmp=kk_gse@result
dim(tmp)
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#按需替换
#DOSE::setReadable():mapping geneID to gene Symbol
tmp=kk@result
dim(tmp)
pro='comp1'
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
save(kk,file = 'gsea_kk.Rdata')

# 展现前6个上调通路和6个下调通路
down_k <- kk_gse[tail(order(kk_gse$enrichmentScore,decreasing = F)),];down_k$group=-1
up_k <- kk_gse[head(order(kk_gse$enrichmentScore,decreasing = F)),];up_k$group=1

dat=rbind(up_k,down_k)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]
# gsea分析结果p7
p7<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) + 
  scale_x_discrete(name ="Pathway names") +
  scale_y_continuous(name ="log10P-value") +
  coord_flip() + 
  theme_ggstatsplot()+
  theme(plot.title = element_text(size = 15,hjust = 0.5),  
        axis.text = element_text(size = 12,face = 'bold'),
        panel.grid = element_blank())+
  ggtitle("Pathway Enrichment") 
p7
#具体看上面条形图里面的每个通路的gsea分布情况p8
library(enrichplot)
p8 <- gseaplot2(kk, geneSetID = rownames(down_k))+
      gseaplot2(kk, geneSetID = rownames(up_k))
p8