(五)04_DEG差异基因

作者: 养猪场小老板 | 来源:发表于2020-01-24 23:15 被阅读0次

    第一,同样的,清除变量+加载之前的数据

    > rm(list = ls()) 
    > load(file = "step2output.Rdata")
    > group_list
    [1] control control control treat   treat   treat  
    Levels: control treat
    

    差异分析,用limma包来做

    需要表达矩阵exp和group_list分组对象,不需要改

    {
    > library(limma)#加载limma包
    #需要自己做好三个数据(表达矩阵,分组矩阵,差异比较矩阵)
    #总共三个步骤(lmFit,eBayes,topTable)
    +  design=model.matrix(~group_list)#把group_list设置成一个model.matrix
    #> design
    #  (Intercept) group_listtreat
    #1           1               0
    #2           1               0
    #3           1               0
    #4           1               1
    #5           1               1
    #6           1               1
    #attr(,"assign")
    #[1] 0 1
    #attr(,"contrasts")
    #attr(,"contrasts")$group_list
    #[1] "contr.treatment"
    +  fit=lmFit(exp,design)#lmFit用于线性拟合;至少两个输入,一个是表达矩阵,一个是分组对象。
    #表达矩阵必须是matrix类数据结构,每一列都是存放一个样本,每一行是一个探针信息或者是注释后的基因名
    +  fit=eBayes(fit)#根据lmFit的拟合结果进行统计推断
    +  deg=topTable(fit,coef=2,number = Inf)#得出差异比较矩阵
    }
    #view(deg)
    
    差异表达矩阵deg.png image.png
    image.png

    为deg数据框添加几列

    {
      #1.加probe_id列:把deg的行名(探针名)赋值给probe_id,并加一列
      library(dplyr)
      deg <- mutate(deg,probe_id=rownames(deg))#mutate增加一列
      head(deg)
      #2.加symbol列,火山图要用
      deg <- inner_join(deg,ids,by="probe_id")
      head(deg)
      #按照symbol列去重复#因为存在一个探针映射多个基因的情况
      deg <- deg[!duplicated(deg$symbol),]
    }
    
    #3.加change列,下面两行是阈值,可修改
    logFC_t=1 
    logP_t = 0.01
    
    test1 = deg$P.Value < logP_t
    test2 = deg$logFC < -logFC_t
    test3 = deg$logFC > logFC_t
    
    {
      change = ifelse(test1 & test2 , 
                   "down" , 
                   ifelse(test1 & test3 ,
                          "up", 
                          "stable"))#不满足down和up时,就是stable;分成三种
      deg <- mutate(deg,change)#加一列有down和up和stable
      #4.加ENTREZID列,用于富集分析(symbol转entrezid,然后inner_join)
      library(ggplot2)
      library(clusterProfiler)
      library(org.Hs.eg.db)
      s2e <- bitr(deg$symbol, fromType = "SYMBOL",
                  toType = "ENTREZID",
                  OrgDb = org.Hs.eg.db)#人类基因ID转换库
      #其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
      deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
    }
    
    save(group_list,deg,file = "step4output.Rdata")
    
    最后数据如图: image.png

    相关文章

      网友评论

        本文标题:(五)04_DEG差异基因

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