文章Tinagl1 Suppresses Triple-negative Breast Cancer Progression by Simultaneously Targeting Integrin/FAK and EGFR Signaling Pathways，根据文章找到GSE号: GSE99507,
实验设计：Tinagl1过表达三个样本+正常对照三个样本。
画出来的图如下：
heatmap.png
1. 安装包，并准备工作
install.packages("GEOquery")
Installing package into ‘/usr/local/lib/R/3.5/site-library’
(as ‘lib’ is unspecified)
Warning in install.packages :
package ‘GEOquery’ is not available (for R version 3.5.0)
- 更新R,https://www.r-project.org, 升级到3.5.2版本
- google后，找解决方案https://blog.csdn.net/yepeng2007fei/article/details/78112557
- 使用如下代码，安装成功
source("[http://bioconductor.org/biocLite.R](http://bioconductor.org/biocLite.R)")
biocLite("GEOquery")
library(GEOquery)
此步骤耗费我1h，本来满满的成就感，直到我看到这篇文章的开头，http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html，平时不看文，实操两行泪
2. 下载数据
- 第三种方法下载数据
gset <- getGEO('GSE99507', destdir=".")
- 第二种下载
a = read.table('文件名', sep = '\t',quote = "", fill = T, comment.char = "!", header = T) - jimmy 现成版，将下载GEO数据集的表达矩阵封装成函数。后续操作过程中，发现将过程封装成函数后，在后续还是要将exprSet调出来，索性不封装了。
3. ID转换
GPL中的探针和基因名称 与 GSE中探针对应的表达关系 相match。
- 下载并转换表达矩阵
rm(list = ls())
options(stringsAsFactors = F) #读表的时候，不要把字符读成factor
library(GEOquery)
eSet <- getGEO('GSE99507', destdir = '.')
exprSet <- exprs(eSet[[1]]) #exprs用来提取表达矩阵
pdata <- pData(eSet[[1]]) #pData用来提取样本信息
View(pdata) #瞅一眼看分组信息
#以下进行分组，分别是Tinagl1过表达组和正常组
library(stringr)
group_list <- trimws(str_split(pdata$title,' ',simplify = T)[,5])
group_list <- str_replace(group_list,'LM2','Vector')
table(group_list)
group_list
- 对芯片进行注释
gpl <- getGEO("GPL17077", destdir = ".") #此GPL没有注释的R包，只能手动找
colnames(Table(gpl))
head(Table(gpl)[,c(1,7)])
write.csv(Table(gpl)[,c(1,7)],"GPL17077.csv")
ids <- read.csv("GPL17077.csv") #获得探针对应基因的表达矩阵
ids = ids[,-1]
View(ids)
- 将表达矩阵中的探针用对应的基因名ID转换
length(unique(ids$GENE_SYMBOL))
tail(sort(table(ids$GENE_SYMBOL)))
table(sort(table(ids$GENE_SYMBOL)))
table(rownames(exprSet) %in% ids$ID)
dim(exprSet)
exprSet <- exprSet[rownames(exprSet) %in% ids$ID,]
dim(exprSet)
ids <- ids[match(rownames(exprSet),ids$ID),] #改ID顺序
View(ids)
head(ids)
exprSet[1:5,1:5] #瞅一眼表达矩阵的开始和ids的是否相同
#ids的GENE_SYMBOL对表达矩阵进行分类
tmp <- by(exprSet,
ids$GENE_SYMBOL,
function(x)
rownames(x)[which.max(rowMeans(x))] )
#分类，如果有一个基因对应多个探针,只保留在所有样本里面平均表达量最大的那个探针
probes <- as.character(tmp)
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% probes ,] #过滤出有多个探针的基因
dim(exprSet)
save(exprSet,group_list,file = "step1.Rdata")
rownames(exprSet)=ids[match(rownames(exprSet),ids$ID),2]
exprSet[1:5,1:5]
View(exprSet)
library(reshape2)
exprSet_L<- melt(exprSet)
colnames(exprSet_L) <- c('probe','sample','value')
exprSet_L$group <- rep(group_list,each=nrow(exprSet))
head(exprSet_L)
View(exprSet_L)
save(exprSet_L,group_list,file = "step2.Rdata")
#看表达矩阵样本分组信息是否合理
exprSet["ACTB",] #关键蛋白常见蛋白β-actin表达量是否合理
library("ggplot2")
p <- ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
4. 画图
- PCA图
install.packages("ggfortify")
library("ggfortify")
df <- as.data.frame(t(exprSet))
df$group <- group_list
autoplot(prcomp(df[,1:(ncol(df)-1)]),data = df, colour = "group")
- 差异分析
library( "limma" )
dim(exprSet)
group_list
design <- model.matrix( ~0 + factor( group_list ) ) #根据分组变成矩阵
colnames( design ) <- levels( factor( group_list ) )
rownames( design ) <- colnames( exprSet )
design
#下面做比较矩阵
contrast.matrix <- makeContrasts(paste0(unique(group_list),
collapse = "-"), levels = design )
contrast.matrix
fit <- lmFit(exprSet, design )
fit2 <- contrasts.fit( fit, contrast.matrix )
fit2 <- eBayes( fit2 )
nrDEG <- topTable( fit2, coef = 1, n = Inf )
write.table( nrDEG, file = "nrDEG.out")
head(nrDEG)
View(nrDEG)
- 画热图
install.packages("pheatmap")
library( "pheatmap" )
choose_gene <- head( rownames( nrDEG ), 20 ) #取了前20个差异最大的基因，这数可以自己定
choose_matrix <- exprSet[ choose_gene, ]
annotation_col = data.frame( CellType = factor( group_list ) )
rownames( annotation_col ) = colnames( exprSet )
pheatmap( fontsize = 5, choose_matrix, annotation_col = annotation_col, show_rownames = T,annotation_legend = T, filename = "heatmap_1.png")
跌跌撞撞，终于花了8h作了一张图出来，使用的都是别人写的代码，主要参考的代码有：
http://www.bio-info-trainee.com/3409.html
https://www.jianshu.com/p/0dcc6030343e
https://mp.weixin.qq.com/s/nVnL_uc6GEkD_cLpMHCYXg
目前还有部分代码不是很理解，还要继续努力啊。
