美文网首页R语言做生信生信分析流程GEO数据挖掘
R语言GEO数据挖掘:步骤三:进行基因差异分析

R语言GEO数据挖掘:步骤三:进行基因差异分析

作者: mayoneday | 来源:发表于2019-03-26 21:55 被阅读71次

    进行基因差异分析

    1.读取数据

    rm(list = ls())  ## 魔幻操作,一键清空~
    options(stringsAsFactors = F)#在调用as.data.frame的时,将stringsAsFactors设置为FALSE可以避免character(字符)类型自动转化为factor(因子,如0,1)类型
    # 注意查看下载文件的大小,检查数据 
    load(file = 'step1-output.Rdata')#加载第一步保存的文件,此文件包括dat和group_list,加载后得到此两数据
    
    dat group_list

    2.检测数据

    # 每次都要检测数据
    dat[1:4,1:4] #查看1到4行和1到4列的dat文件
    table(group_list) #table函数,查看group_list中的分组个数
    
    dat[1:4,1:4] table(group_list)

    3.为每个数据集绘制箱形图,比较数据集中的数据分布(单一某个基因分析)

    boxplot(dat[1,]~group_list) #按照group_list分组画箱线图,定义分组,(此步骤可以省略)
    
    bp=function(g){         #定义一个函数g,函数为{}里的内容
      library(ggpubr)
      df=data.frame(gene=g,stage=group_list)#把分组和数据糅合在一起
      p <- ggboxplot(df, x = "stage", y = "gene",#定义图形的横纵坐标
                     color = "stage", palette = "jco",
                     add = "jitter")#
        p + stat_compare_means()#  Add p-value值
    }
    bp(dat[1,]) ## 调用上面定义好的函数,避免同样的绘图代码重复多次敲。
    bp(dat[2,])
    bp(dat[3,])
    bp(dat[4,])
    dim(dat)
    

    4.对某种基因在样本中的表达情况根据分组统一分析(所有基因一起分析)

    4.1数据表达形式

    用limma包,这里注意,limma包是对基因芯片表达矩阵的分析,不能对逆转录RNAseq表达矩阵进行分析(因为数据特征不同),RNAseq需要用另一种方法

    library(limma)
    design=model.matrix(~factor( group_list ))#把分组信息变化成因子
    
    design
    fit=lmFit(dat,design)#在给定一系列阵列的情况下,拟合每个基因的线性模型
    
    QQ截图20190326210520.jpg
    fit=eBayes(fit)#eBayes函数指在微阵列线性模型拟合下,通过经验Bayes对标准误差向一个共同值的缓和,计算缓和t-统计、缓和f-统计和微分表达式的对数概率。 用法举列:EBAYES(配合,比例=0.01,标准偏差,系数,极限值=C(0.1,4) 
    ## 上面是limma包用法的一种方式 
    
    #把上面计算出的差异再进行处理
    options(digits = 4) #设置全局的数字有效位数为4
    #topTable(fit,coef=2,adjust='BH') 
    topTable(fit,coef=2,adjust='BH') #提取一个表的列上的基因从一个线性模型的拟合
    #例子:topTable(fit, coef=NULL, number=10, genelist=fit$genes, adjust.method="BH",
             #sort.by="B", resort.by=NULL, p.value=1, lfc=0, confint=FALSE)  
    
    所有基因差异总数据表

    解读此表


    表中数据意义

    但是上面的用法做不到随心所欲的指定任意两组进行比较,所有还有下一种方法

    design <- model.matrix(~0+factor(group_list))#把分组变成因子形式
    
    
    以分组信息为因子得到design
    colnames(design)=levels(factor(group_list))#改design列名为分组信息
    
    改之后列名变了
    head(design)#查看分组文件
    exprSet=dat#exprSet 为表达矩阵
    rownames(design)=colnames(exprSet)#更改分组信息design行名为分组名,也就是样本名字
    design
    #构造了design矩阵,得出每一个样本属于哪一个组,属于为1,不属于为0
    
    
    QQ截图20190326213059.jpg

    处理好了分组信息,再自定义比较元素

    #自定义比较元素
    contrast.matrix<-makeContrasts("noTNBC-TNBC",
                                   levels = design)
    #也可以用下面的,和上面等同意思
    #contrast.matrix<-makeContrasts("paste0(unique(group_list),collapse = "-"),levels = design)
    contrast.matrix ##这个矩阵声明,我们要把 Tumor 组跟 Normal 进行差异分析比较
    

    自定义函数进行比较

    deg = function(exprSet,design,contrast.matrix){#定义一个叫deg的函数
      ##step1#logFC后是前的多少倍,进行差异化比较
      fit <- lmFit(exprSet,design)
      ##step2
      fit2 <- contrasts.fit(fit, contrast.matrix) 
      ##这一步很重要,大家可以自行看看效果
      
      fit2 <- eBayes(fit2)  ## default no trend !!!
      ##eBayes() with trend=TRUE
      ##step3
      tempOutput = topTable(fit2, coef=1, n=Inf)
      nrDEG = na.omit(tempOutput) 
      #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
      head(nrDEG)
      return(nrDEG)
    }
    deg = deg(exprSet,design,contrast.matrix)
    head(deg)#查看前几行文件
    save(deg,file = 'deg.Rdata')#存储文件
    #得出了两两差异化比较的图
    
    得出所有基因差异化数据表

    4.2分析基因表达差异的可视化表示:火山图,热图

    for volcano 火山图

    if(T){
      nrDEG=deg
      head(nrDEG)
      attach(nrDEG)
      plot(logFC,-log10(P.Value))
      library(ggpubr)
      df=nrDEG
      df$v= -log10(P.Value) #df新增加一列'v',值为-log10(P.Value)
      ggscatter(df, x = "logFC", y = "v",size=0.5)
      
      df$g=ifelse(df$P.Value>0.01,'stable
                  ', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
                  ifelse( df$logFC >2,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
                          ifelse( df$logFC < -2,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1.5 的基因,再if 判断:如果logFC <1.5,则为down(下调)基因,否则为stable基因
      )
      table(df$g)
      df$name=rownames(df)
      head(df)
      ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
      ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
                label = "name", repel = T,
                #label.select = rownames(df)[df$g != 'stable'] ,
                label.select = c('TTC9', 'AQP3', 'CXCL11','PTGS2'), #挑选一些基因在图中显示出来
                palette = c("#00AFBB", "#E7B800", "#FC4E07") )
      ggsave('volcano.png')
      
      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") )
      ggsave('MA.png')
      
      
    }
    

    for heatmap

    if(T){ 
      load(file = 'step1-output.Rdata')
      # 每次都要检测数据
      dat[1:4,1:4]
      table(group_list)
      x=deg$logFC #deg取logFC这列并将其重新赋值给x
      names(x)=rownames(deg) #deg取probe_id这列,并将其作为名字给x
      cg=c(names(head(sort(x),100)),#对x进行从小到大排列,取前100及后100,并取其对应的探针名,作为向量赋值给cg
           names(tail(sort(x),100)))
      library(pheatmap)
      pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行,所得到的矩阵来画热图
      n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于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) #将ac的行名也就分组信息(是‘no TNBC’还是‘TNBC’)给到n的列名,即热图中位于上方的分组信息
      pheatmap(n,show_colnames =F,
               show_rownames = F,
              cluster_cols = F, 
               annotation_col=ac,filename = 'heatmap_top200_DEG.png') #列名注释信息为ac即分组信息
      
      
    }
    
    write.csv(deg,file = 'deg.csv')
    dev.new()
    

    热土和火山图都是傻瓜式的,只要的前面得出的deg数据(也就是基因差异表达数据)是正确的

    最后

    感谢jimmy的生信技能树团队!

    感谢导师岑洪老师!

    感谢健明、孙小洁,慧美等生信技能树团队的老师一路以来的指导和鼓励!

    特别注明:此文中编码来自生信技能树健明老师

    相关文章

      网友评论

        本文标题:R语言GEO数据挖掘:步骤三:进行基因差异分析

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