小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看

搜索
分享

【原】萌新學完GEO課程復現(xiàn)SCI文章差異基因的熱圖

 健明 2021-07-14

文章標題是：Tinagl1 Suppresses Triple-Negative Breast Cancer Progression and Metastasis by Simultaneously Inhibiting Integrin/FAK and EGFR Signaling .

菜鳥一枚，所以現(xiàn)在只是復現(xiàn)這篇文章里的B圖；C圖里的通路GO/KEGG注釋及A圖中的GSEA下回分解?? （感覺任重而道遠）

從文章中，我們找到GSE號

All microarray data generated in this study have been deposited as a superseries at the NCBI Gene Expression Omnibus with the accession code GSE99507, and can be accessed at https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE99507.

點開，看到如下分組，這里面的分組信息先留個懸念

<從下面開始，跟著老大的嗶哩嗶哩GEO挖掘視頻一步一步走>

1. 現(xiàn)在我們需要找到這個GSE號對應的GPL平臺，然后找到這個GPL平臺所對應的注釋R包

老大GEO視頻中一閃而過的文章是:從GEO數(shù)據(jù)庫下載得到表達矩陣 一文就夠

現(xiàn)在我們需要通這個GPL平臺找到對應的注釋R包然后找到下面這個包裝函數(shù)的板塊，截圖如下：

圖片中是老大已經(jīng)包裝好的函數(shù)，我們這次要下載的GSE號是99507，因此在R里代碼框輸入下面即可，如圖：

得到當前GSE號所對應的芯片平臺，如下圖：即GPL17077.

這個時候，我們需要非常熟悉這個對象哦

我們也可以在GEO網(wǎng)站上看到平臺信息是這樣的，截圖如下

接下來，因此我們可以去老大的搜集貼找對應關系，參考老大的這篇文章：(16)芯片探針與基因的對應關系-生信菜鳥團博客2周年精選文章集，截圖如下：

從以上列表中，用GPL17077去找對應的注釋包，發(fā)現(xiàn)沒有。（老大已經(jīng)總結地很全了，但是依然有些GPL平臺是沒有注釋包或者是有我們沒有找到。）

而老大GEO視頻里的GPL平臺其實是有注釋R包的，就是下圖中的這個hgu95av2.db，視頻中老大從R語言20練習題里有關于ID轉換的代碼，如下圖：

總之，現(xiàn)在的GPL17077平臺現(xiàn)在并沒有對應的注釋R包，當然更不是視頻中的這個hgu95av2.db了。為了能繼續(xù)跟著按照老大視頻中講的繼續(xù)做，我去Google搜索，用的是前面在GEO網(wǎng)站里GPL平臺對應的平臺信息，前面也有截圖，如下

Google找的說是這個hgug4110b.db的注釋R包，其實在老大總結里這個包對應的是GPL887平臺，但是我想還是試試，于是我嘗試了下載下，代碼如下啦

# 首先我下載了一個錯誤的包
BiocManager::install('hgug4110b.db')
library('hgug4110b.db')
?hgug4110b
ls("package:hgug4110b.db")

#下面這段很眼熟，就是我上面在R語言20題里的截圖，我把視頻中的hgu95av2.db包換成我在谷歌搜索到的hgug4110b.db包
library('hgug4110b.db')
ids=toTable(hgug4110bSYMBOL)
length(unique(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))


exprSet<-dat
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(exprSet)


ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp = by(exprSet,ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)
#寫一個函數(shù)，可以之間過濾掉其他探針并且進行ID轉換，
#此時的id轉換是用bioconductor包來轉換的，如果下載的matrix，則不用轉換
jimmy <- function(exprSet,ids){
                             tmp = by(exprSet,
                             ids$symbol,
                              function(x) rownames(x)[which.max(rowMeans(x))] )
                probes = as.character(tmp)
                print(dim(exprSet))
                exprSet=exprSet[rownames(exprSet) %in% probes ,]
                print(dim(exprSet))
                return(exprSet)
}
new.exprSet <- jimmy(exprSet,ids)


#此時會報錯，因為現(xiàn)在的exprSet和ids的長度不一樣

由于expeSet里的ID能在ids里對應的很少，大部分都是對應不上的，如下圖

對應上的基因數(shù)TRUE只有8190，因此我認為我谷歌到的這個hgug4110b.db并不對。

不過==沒關系==鼓起勇氣，繼續(xù)探索

而且老大視頻里又給了找不到對應的注釋R包時，通過下面的代碼來實現(xiàn)探針名和基因名的對接，代碼網(wǎng)址：https://github.com/jmzeng1314/my-R/blob/master/9-microarray-examples/E-MTAB-3017.R

#代碼如下
gpl <- getGEO('GPL17077', destdir=".")
colnames(Table(gpl)) 
head(Table(gpl)[,c(1,6,7)])  ## you need to check this , which column do you need 
write.csv(Table(gpl)[,c(1,6,7)],"GPL17077.csv")
ids=read.csv('GPL17077.csv')
#這樣我們就同樣獲得了探針對應基因的表達矩陣啦

前面我所認為的第一步剛剛結束了，獲得了探針對應基因的表達矩陣

2.下載數(shù)據(jù)，獲得分組信息

下載數(shù)據(jù)：老大在視頻里講了三種找到數(shù)據(jù)集后的數(shù)據(jù)下載方式，我簡單羅列如下（目的得到表達矩陣）：

1.直接下載raw data，但不推薦大家用，原始數(shù)據(jù)

2.下載表達矩陣 series matrix file(s)，下載后可讀到R里面，考驗網(wǎng)速

a=read.table('GSE42872_series_matrix.txt.gz')
> class(a)
[1] "data.frame"
> str(gset)

3.在R里面讀取GSE號.

gset <- getGEO("GSE42589")

加載GEO包下載

library(GEOquery)
gset <- getGEO('GSE42872',destdir=".",AnnotGPL = F,getGPL = F) 

老大推薦第三種方式，下面是代碼

# 就這么來就對了，沒毛病
###step1-downloadR
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù) 
f='GSE99507_eSet.Rdata'

library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
  gset <- getGEO('GSE99507', destdir=".",
                 AnnotGPL = F,     ## 注釋文件
                 getGPL = F)       ## 平臺文件
  save(gset,file=f)   ## 保存到本地
}
load('GSE99507_eSet.Rdata')  ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs
dim(dat)
#  GPL3921  [HT_HG-U133A] Affymetrix HT Human Genome U133A Array
# Bioconductor - hgu133a.db
dat[1:4,1:4]


pd=pData(a) #獲得分組信息，就得到了前面一張截圖，即前面懸念的解答處-GSE號對應的Samples分組信息，主要是根據(jù)文章里的那個截圖的分組信息來對字符串進行拆分（str_split）和替換（str_replace）
library(stringr)
group_list = trimws(str_split(pd$title,' ',simplify = T)[,5])#拆分
group_list = str_replace(group_list,'LM2','Vector')
table(group_list)
save(dat,group_list,file = 'step1-output.Rdata')
#此時，要注意保存的這個文件step1-output.Rdata，我們要知道此時的dat表達矩陣是什么樣的

3.ID轉換

###step2-ids-filter.R

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load('step1-output.Rdata')
gpl <- getGEO('GPL17077', destdir=".")
colnames(Table(gpl)) ## [1] 41108    17
head(Table(gpl)[,c(1,7)])  ## you need to check this , which column do you need 
write.csv(Table(gpl)[,c(1,7)],"GPL17077.csv")
ids=read.csv('GPL17077.csv')
ids=ids[,-1]
dat[1:4,1:4] 
rownames(dat)

#如前所述，現(xiàn)在沒有這個平臺對應的注視R包，所以下面是根據(jù)下載的dat和ids來進行的ID轉換步驟
ids[ids$GENE_SYMBOL!='',]
ids=ids[ids$GENE_SYMBOL!='',]
match(rownames(dat),ids$ID)
ids[match(rownames(dat),ids$ID),]
ids=ids[match(rownames(dat),ids$ID),]
ids=ids[complete.cases(ids),]
dat=dat[ids$ID,]

#下面是針對那些有多個基因去重的步驟，注意理解為什么要取median（中位數(shù)），巧妙
ids$median=apply(dat,1,median)
ids=ids[order(ids$GENE_SYMBOL,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$GENE_SYMBOL),]
dat=dat[ids$ID,]
rownames(dat)=ids$GENE_SYMBOL
dat[1:4,1:4]  
dim(dat)

save(dat,group_list,file = 'step2-output.Rdata')

#此時此時，要注意保存的這個文件step2-output.Rdata，我們要知道此時的dat表達矩陣是什么樣的

4.PCA和Heatmap

###step3-pca-heatmap

##PCA主成分分析

(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
## 下面是畫PCA的必須操作，需要看說明書。
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,group_list)
library("FactoMineR")
library("factoextra") 
# The variable group_list (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,repel =T,
             #geom.ind = "point", # show points only (nbut not "text")
             col.ind = dat$group_list, # color by groups
             palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)



##熱圖：heatmap-sd

rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
cg=names(tail(sort(apply(dat, 1, sd)),1000))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) 
cg1=names(tail(sort(rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) ),500))
cg2=names(head(sort(rowMeans(dat[,1:3]) -  rowMeans(dat[,4:6]) ),500))
cg=c(cg1,cg2)
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,
         show_rownames = F,
         cluster_cols = F,
         annotation_col=ac)

5.差異分析

##step4-DEG
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
load(file = 'step2-output.Rdata')######要知道此時是step2中的dat哦
# 每次都要檢測數(shù)據(jù)
dat[1:4,1:4]
table(group_list)

##boxplot圖
boxplot(dat[1,]~group_list)
t.test(dat[1,]~group_list)
bp=function(g){
  library(ggpubr)
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
bp(dat[1,])
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])

##DEG-limma包

library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
#topTable(fit,coef=2,adjust='BH') 
topTable(fit,coef=2,adjust='BH')
bp(dat['KLK10',])
bp(dat['FXYD3',])
deg=topTable(fit,coef=2,adjust='BH',number = Inf)

head(deg) 
save(deg,file = 'deg.Rdata')


## for volcano 
if(T){
  nrDEG=deg
  head(nrDEG)
  attach(nrDEG)
  plot(logFC,-log10(P.Value))
  library(ggpubr)
  df=nrDEG
  df$v= -log10(P.Value)
  ggscatter(df, x = "logFC", y = "v",size=0.5)

  df$g=ifelse(df$P.Value>0.01,'stable',
              ifelse( df$logFC >1.5,'up',
                      ifelse( df$logFC < -1.5,'down','stable') )
  )
  table(df$g)#此時獲得的上調基因為20，下調基因為145
  df$symbol=rownames(df)
  ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
  ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
            label = "symbol", repel = T,
            #label.select = rownames(df)[df$g != 'stable'] ,
            label.select = rownames(head(head(deg))),
            palette = c("#00AFBB", "#E7B800", "#FC4E07") )

  ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
  df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
                  ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
  table(df$p_c )
  ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2, 
            palette = c("green", "red", "black") )


}


## 熱圖：for heatmap 
rm(list = ls())  ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
if(T){ 
  load(file = 'step2-output.Rdata')
  # 每次都要檢測數(shù)據(jù)
  dat[1:4,1:4]
  table(group_list)
  load(file = 'deg.Rdata')
  x=deg$logFC
  names(x)=rownames(deg)
  cg=c(names(head(sort(x),20)),#老大的代碼中是100、100，由于我獲得上下調基因少，就改了 
       names(tail(sort(x),145)))
  library(pheatmap)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
  n=t(scale(t(dat[cg,])))
  n[n>2]=2
  n[n< -2]= -2
  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(g=group_list)
  rownames(ac)=colnames(n)
  pheatmap(n,show_colnames =F,
           show_rownames = F,
           cluster_cols = F,
           annotation_col=ac)


}

到這里就差不多結束了，基本得到和文章中差不多的熱圖。還有很多做圖技巧需要慢慢學習。其中的代碼均從老大那copy過來的，需要邊做邊理解，注意每一步變量的變化，以不變應萬變，這點很重要哇。

十分<感謝老大>的指導和幫助，第一次做到這里，很開心哇～??

    轉藏 分享 獻花（0）

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多