美文网首页生信GEO数据挖掘
GEO实战演练(一)

GEO实战演练(一)

作者: ShanSly | 来源:发表于2019-07-15 00:21 被阅读153次

    《Epidermal Growth Factor Receptor Activation: An Upstream Signal for Transition of Quiescent Astrocytes into Reactive Astrocytes after Neural Injury》

    首先ncbi中检索该文章,得到数据GSE5282

    第一步:Start:获取数据并得到表达矩阵

    rm(list = ls())
    options(stringsAsFactors = F) 
    library(GEOquery)
    eSet <- getGEO("GSE5282", 
                   destdir = '.',   
                   getGPL = F)
    class(eSet)
    length(eSet)
    class(eSet[[1]]) 
    exp <- exprs(eSet[[1]])
    exp[1:4,1:4] 
    exp1 <- exp[,c(4,5,7,1,2)]
    exp1[1:4,1:3]
    pd <- pData(eSet[[1]])
    pd
    save(pd,exp,exp1,file = "step1output.Rdata")
    eSet[[1]]@annotation
    

    第二步:设置分组信息

    ##control:EGF4h
    m(list = ls())  
    options(stringsAsFactors = F)
    load(file = "step1output.Rdata")
    class(pd)
    dim(pd)
    colnames(pd)
    group_list1=rep(c('control','EGF_4h'),c(3,2)) 
    group_list1 
    
    save(group_list1,file = "step2output.Rdata")
    

    第三步:检验组间数据的显著性(PCA)

    rm(list = ls())  ## 魔幻操作,一键清空~
    options(stringsAsFactors = F)
    load(file = "step1output.Rdata")
    load(file = "step2output.Rdata")
    exp1[1:4,1:4]
    dat1=as.data.frame(t(exp1))
    dim(exp1)
    dim(dat1)
    dat1=cbind(dat1,group_list1)
    
    library(FactoMineR)
    library(factoextra)
    dat.pca <- PCA(dat1[,-ncol(dat1)], graph = FALSE)
    fviz_pca_ind(dat.pca,
                 geom.ind = "point", 
                 col.ind = dat1$group_list1, 
                 palette = c("#00AFBB", "#E7B800"),
                 addEllipses = TRUE, 
                 legend.title = "Groups") 
    cg=names(tail(sort(apply(exp,1,sd)),1000)) 
    n=t(scale(t(exp[cg,]))) 
    thr=2
    n[n>thr]=thr    
    n[n< -thr]= -thr  
    n[1:4,1:4]
    pheatmap(n,show_colnames =F,show_rownames = F)
    

    Result:


    PCA

    通过图示,很明显发现两者的组间差异明显,数据合理。

    第四步:筛选差异基因

    rm(list = ls())  
    options(stringsAsFactors = F)
    options(digits = 4)   
    load(file = "step1output.Rdata")
    load(file = "step2output.Rdata")
    exp1[1:4,1:4]
    table(group_list1) 
    boxplot(exp1[1,]~group_list1)
    
    bp=function(g){         
      library(ggpubr)
      df=data.frame(gene=g,group=group_list1)
      p <- ggboxplot(df, x = "group", y = "gene",
                     color = "group", palette = "jco",
                     add = "jitter")
      p + stat_compare_means(label.y = 8)
    }
    bp(exp1[1,])
    
    library(limma)
    design=model.matrix(~factor(group_list1))
    fit=lmFit(exp1,design)
    fit=eBayes(fit)
    deg1=topTable(fit,coef=2,number = Inf)
    head(deg1)
    deg1$logFC <- (deg1$logFC)
    bp(exp1[rownames(deg1)[1],])
    bp(exp1[rownames(deg1)[2],])
    library(dplyr)
    deg1 <- mutate(deg1,probe_id=rownames(deg1))  
    
    head(deg1)
    library(rat2302.db)
    ls("package:rat2302.db") 
    ids <- toTable(rat2302SYMBOL)  
    head(ids)
    
    deg1 <- inner_join(deg1,ids,by="probe_id")
    head(deg1)
    logFC_t=1
    change=ifelse(deg1$P.Value>0.01,'stable', 
                  ifelse( deg1$logFC >logFC_t,'up', 
                          ifelse( deg1$logFC < -logFC_t,'down','stable') )
    )
    
    deg1 <- mutate(deg1,change)
    head(deg1)
    table(deg1$change)
    
    
    library(ggplot2)
    library(clusterProfiler)
    library(org.Hs.eg.db)
    #s2e <- bitr(unique(deg1$symbol), fromType = "SYMBOL",
                toType = c( "ENTREZID"),
                OrgDb = org.Hs.eg.db)
    
    
    #######################################################
    head(s2e)
    head(deg1)
    deg1 <- inner_join(deg1,s2e,by=c("symbol"="SYMBOL"))
    
    head(deg1)
    save(group_list1,deg1,file = "step4output.Rdata")
    
    

    Q:s2e <- bitr(unique(deg1$symbol), fromType = "SYMBOL", toType = c( "ENTREZID"), OrgDb = org.Hs.eg.db)deg1 <- inner_join(deg1,s2e,by=c("symbol"="SYMBOL"))
    得到的很尴尬,全部成了stable基因,如果用这个deg1去画热图 就天下大乱了。

    第五步:火山图

    rm(list = ls())  
    options(stringsAsFactors = F)
    load(file = "step4output.Rdata")
    
    head(deg1)
    
    plot(deg1$logFC,-log10(deg1$P.Value))
    
    library(ggpubr)
    
    table(dat1$change)
    ggscatter(dat1, x = "logFC", y = "v",size=0.5,color = "change")
    
    
    ggscatter(dat1, x = "logFC", y = "v", color = "change",size = 0.5,
              label = "symbol", repel = T,
              label.select = dat1$symbol[1:30] ,
              #label.select = c('CD36','DUSP6'), #挑选一些基因在图中显示出来
              palette = c("#00AFBB", "#999999", "#FC4E07") )
    

    Result:


    volcano

    第五步:heatmap

    rm(list = ls())  ## 魔幻操作,一键清空~
    options(stringsAsFactors = F)
    load(file = 'step1output.Rdata')
    load(file = 'step2output.Rdata')
    load(file = 'step4output.Rdata')
    
    
    x=deg1$logFC 
    names(x)=deg1$probe_id 
    cg=c(names(head(sort(x),100)),
         names(tail(sort(x),100)))
    
    library(pheatmap)
    pheatmap(exp1[cg,],show_colnames =F,show_rownames = F)
    
    n=t(scale(t(exp1[cg,])))
    
    thr=2
    n[n>thr]=thr 
    n[n< -thr]= -thr
    
    n[1:4,1:4]
    pheatmap(n,show_colnames =F,show_rownames = F)
    #显示分组
    ac=data.frame(group=group_list1)
    rownames(ac)=colnames(n) 
    
    pheatmap(n,show_colnames =T,
             show_rownames = F,
             annotation_col = ac,color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
    #保存
    pdf(file = "heatmap.pdf")
    pheatmap(n,show_colnames =F,
             show_rownames = F,
             #cluster_cols = F,
             annotation_col=ac
    )
    dev.off()
    
    

    Result:


    heatmap

    相关文章

      网友评论

        本文标题:GEO实战演练(一)

        本文链接:https://www.haomeiwen.com/subject/llbwkctx.html