文章Tinagl1 Suppresses Triple-negative Breast Cancer Progression by Simultaneously Targeting Integrin/FAK and EGFR Signaling Pathways，根據(jù)文章找到GSE號: GSE99507,
實(shí)驗(yàn)設(shè)計(jì)：Tinagl1過表達(dá)三個樣本+正常對照三個樣本。
畫出來的圖如下：

heatmap.png

1. 安裝包，并準(zhǔn)備工作

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)

此步驟耗費(fèi)我1h，本來滿滿的成就感，直到我看到這篇文章的開頭，http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html，平時不看文，實(shí)操兩行淚

2. 下載數(shù)據(jù)

• 第三種方法下載數(shù)據(jù)

gset <- getGEO('GSE99507', destdir=".")

• 第二種下載
  a = read.table('文件名', sep = '\t',quote = "", fill = T, comment.char = "!", header = T)
• jimmy 現(xiàn)成版，將下載GEO數(shù)據(jù)集的表達(dá)矩陣封裝成函數(shù)。后續(xù)操作過程中，發(fā)現(xiàn)將過程封裝成函數(shù)后，在后續(xù)還是要將exprSet調(diào)出來，索性不封裝了。

3. ID轉(zhuǎn)換

GPL中的探針和基因名稱 與 GSE中探針對應(yīng)的表達(dá)關(guān)系 相match。

• 下載并轉(zhuǎn)換表達(dá)矩陣

rm(list = ls())
options(stringsAsFactors = F) #讀表的時候，不要把字符讀成factor
library(GEOquery)
eSet <- getGEO('GSE99507', destdir = '.')
exprSet <-  exprs(eSet[[1]]) #exprs用來提取表達(dá)矩陣
pdata <- pData(eSet[[1]]) #pData用來提取樣本信息
View(pdata) #瞅一眼看分組信息

#以下進(jìn)行分組，分別是Tinagl1過表達(dá)組和正常組
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

• 對芯片進(jìn)行注釋

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") #獲得探針對應(yīng)基因的表達(dá)矩陣
ids = ids[,-1]
View(ids)

• 將表達(dá)矩陣中的探針用對應(yīng)的基因名ID轉(zhuǎn)換

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] #瞅一眼表達(dá)矩陣的開始和ids的是否相同

#ids的GENE_SYMBOL對表達(dá)矩陣進(jìn)行分類
tmp <- by(exprSet,
         ids$GENE_SYMBOL,
         function(x)
         rownames(x)[which.max(rowMeans(x))] ) 
#分類，如果有一個基因?qū)?yīng)多個探針,只保留在所有樣本里面平均表達(dá)量最大的那個探針
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")

#看表達(dá)矩陣樣本分組信息是否合理
exprSet["ACTB",] #關(guān)鍵蛋白常見蛋白β-actin表達(dá)量是否合理
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 ) ) #根據(jù)分組變成矩陣
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個差異最大的基因，這數(shù)可以自己定
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
http://www.itdecent.cn/p/0dcc6030343e
https://mp.weixin.qq.com/s/nVnL_uc6GEkD_cLpMHCYXg
目前還有部分代碼不是很理解，還要繼續(xù)努力啊。