文章Tinagl1 Suppresses Triple-negative Breast Cancer Progression by Simultaneously Targeting Integrin/FAK and EGFR Signaling Pathways悯嗓，根據(jù)文章找到GSE號: GSE99507,
實驗設(shè)計：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)
此步驟耗費我1h，本來滿滿的成就感，直到我看到這篇文章的開頭跛十，http://www.bio-info-trainee.com/bioconductor_China/software/GEOquery.html唐含，平時不看文浅浮，實操兩行淚
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.reibang.com/p/0dcc6030343e
https://mp.weixin.qq.com/s/nVnL_uc6GEkD_cLpMHCYXg
目前還有部分代碼不是很理解，還要繼續(xù)努力啊蹬碧。
