|
原文名称：Identification of Key Genes and Pathways in Triple-Negative Breast Cancer by Integrated Bioinformatics Analysis rm(list = ls()) options()$repos #更改镜像options('repos' = c(CRAN='https://mirrors.tuna./CRAN/'))options()$BioC_mirroroptions(BioC_mirror='https://mirrors.ustc.edu.cn/bioc/')options('repos' = c(CRAN='https://mirrors.tuna./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个TNNOT_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] != '',] # 去除空白的基因名 留下的数据赋值给idsdim(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] #二选一 已经删除未被注释的IDsave(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)) #行太大，不要Viewdim(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)designcontrast.matrix <- makeContrasts('TN-not_TN', levels = design) contrast.matrixload('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)## 之后可做热图，可画火山图，可行注释参考来源：生信技能树 |
|
|
