对TNBC做GEO数据挖掘

原文名称：Identification of Key Genes and Pathways in Triple-Negative Breast Cancer by Integrated Bioinformatics Analysis
GEO编号：GSE76275

rm(list = ls()) 
options()$repos #更改镜像
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

## 下载GEO数据
library(GEOquery) #需提前安装
gset <- getGEO('GSE76275', destdir=".",
               AnnotGPL = F,     ## 注释文件
               getGPL = F)       ## 平台文件 由于我这边下载GPL总是报错，故F，后面另下载。

gset <- gset[[1]] #降级处理
exprSet <- exprs(gset) #获取表达矩阵
dim(exprSet)
pdata <- pData(gset) #获取样本信息
colnames(pdata)
group_list <- pdata[,67] #通过观察得知该分组信息在pdata的第67列，提取出来
table(group_list) # 对group_list进行统计，有67个not TN，198个TN
NOT_TN_expr <- exprSet[,grep("not TN",group_list)] #根据group_list的分组把not TN的表达矩阵提取出来
TN_expr     = exprSet[, !(colnames(exprSet) %in% colnames(NOT_TN_expr))] #根据group_list的分组把TN的表达矩阵提取出来
dim(NOT_TN_expr)
dim(TN_expr)
exprSet <- cbind(NOT_TN_expr,TN_expr) #将两者重新合并，得到的是顺序排列的表达矩阵 此表达矩阵原本就按照先not TN，后TN的顺序排列，这一步骤是为了规范化，对无序的表达矩阵同样适用
group_list = c(rep('not_TN',ncol(NOT_TN_expr)), 
               rep('TN', ncol(TN_expr))) #在排好序的表达矩阵的基础上，对group_list进行排序 注意理解 两组无序的数据通过table等使之有序
table(group_list)
save(exprSet,group_list,file = "exprSet_by_group.Rdata") #保存排好序的表达矩阵和分组信息
load("exprSet_by_group.Rdata") #加载

GPL570 <- getGEO("GPL570",destdir = ".") #下载注释信息
GPL570_a <- Table(GPL570) #提取注释信息
GPL570_b <- GPL570_a[,c(1,11)] #通过观察 第一列ID和第11列gene symbol为我们所需 提取出来
ids = GPL570_b[GPL570_b[,2] != '',] # 去除空白的基因名 留下的数据赋值给ids
dim(ids)

# 部分探针对应多个基因，需要将其提取出来
library(plyr)
a <- strsplit(as.character(ids[,2]),"///") #将基因用///分割
tmp <- mapply(cbind, ids[,1],a) #分割后的基因与探针重新合并 成三列
df <- ldply(tmp,data.frame) #将tmp转化为数据框
ID2gene <- df[,2:3] #二选一 已经删除未被注释的ID
save(ID2gene,file="ID2gene.Rdata") #保存
load("ID2gene.Rdata")

## 将表达矩阵中未被注释的探针去除
exprSet <- exprSet[row.names(exprSet)%in%ID2gene[,1],] #对两者进行匹配
ID2gene <- ID2gene[match(rownames(exprSet),ID2gene[,1]),]
colnames(ID2gene) <- c("id","gene") #修改ID2gene的列名
dim(exprSet)
dim(ID2gene) #两者列数相同

tail(sort(table(ID2gene$gene))) # 同一个基因存在不同的表达量，我们需要对其进行筛选，一般只保留最大值
{
  MAX = by(exprSet, ID2gene[,2], 
           function(x) rownames(x)[ which.max(rowMeans(x))])
  MAX = as.character(MAX)    
  exprSet = exprSet[rownames(exprSet) %in% MAX,]  #筛选出最大值的表达矩阵
  rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]
}
dim(exprSet)
exprSet[1:5,1:5]  
save(exprSet, group_list, file = 'finalSet.Rdata')  ## 数据调整全部完成，保存最终数据

## 做PCA分析
load("finalSet.Rdata")
library(ggfortify) #加载R包，需提前安装
data <- as.data.frame(t(exprSet)) #行太大，不要View
dim(data)
data$group <- group_list #给data加列
png("pca_plot.png",res=100)
autoplot(prcomp(data[, 1:(ncol(data)-1)]), data = data, colour = 'group') + theme_bw() 
dev.off() ## 保存png图，完成PCA分析，可以查看一下png

## 用limma包做差异分析
library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
row.names(design) <- colnames(exprSet)
design
contrast.matrix <- makeContrasts("TN-not_TN", levels = design) 
contrast.matrix
load("ID2gene.Rdata")
{
  fit <- lmFit(exprSet, design)  
  fit2 <- contrasts.fit(fit, contrast.matrix)  
  fit2 <- eBayes( fit2 )   ## Empirical Bayes Statistics for Differential Expression
  nrDEG = topTable( fit2, coef = 1, n = Inf ) #产生差异表达矩阵
  write.table( nrDEG, file = "nrDEG.out")
}
head(nrDEG)

## 之后可做热图，可画火山图，可行注释

参考来源：生信技能树

友情链接：

课程分享
生信技能树全球公益巡讲
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小时生信工程师教学视频合辑
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招学徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）

欢迎关注公众号：青岛生信菜鸟团