美文网首页GEO数据挖掘RNA CHIP ATAC -seq
对于多表型数据如何利用limma包进行差异基因分析

对于多表型数据如何利用limma包进行差异基因分析

作者: Kevin_Hhui | 来源:发表于2020-11-20 14:30 被阅读0次

    写在前面
    最近在做GSVA,做到我有点自闭,不是因为这个算法有多难,而是他所耗费的时间有点夸张。加上校园网这两天时不时“调皮”一下,做了两次,等了四天,最后都白白浪费了。
    其实后面去查了挺多资料,其实可以尝试用标准化的数据(data),而不用原始数据counts,虽然官网推荐利用raw count做,但是它耗费的时间太多了。

    吐槽完,来讲一下今天的正文吧

    对于一些多分组数据,我们该如何利用limma包进行差异基因分析

    其实也因为我后面有这个需求,才想着做这个的分享。

    我们都知道

    使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组。这时,如果逐一使用两两比较求差异基因则略显复杂。其实开发limma包的大神们已经替我们考虑到。我自己當下limma包的PDF,仔细研读并将代码运行后分享给有需要的同道,相互学习。

    1. 首先加载前期预处理好的表达矩阵
    library(limma)##加载包
    load("GSE37761exp_groupfile-00.RData")#加载匹配过的基因表达矩阵
    eset <- exp ##将基因表达矩阵赋值给eset
    

    2.创建样本分类信息表

    targets <- read.csv("group_file.csv",row.names = 1)##读入样本数据,包括两列,
    #第一列GSM号,第二列为样本分组#该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则
    targets$Target <- gsub("_",".",targets$Target)##将"_"替换成“.”,也可以不替换
    ##该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则,所以在没个分类前加一个“F”,具体自己定
    targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))
    colnames(targets)=c("FileName","Target")#更改列名,为了和limma包中的一致
    lev<-unique(targets$Target)##使用unique()函数进行去重
    f <- factor(targets$Target, levels=lev) 
    design <- model.matrix(~0+f) #样本矩阵
    colnames(design) <- lev #更改列名为levels名
    

    这个看起来可能有点费力,换一种更清晰的形式

    group_list<-read.table("gene_list.txt") #导入分组数据
    new_group<- group_list[order(group_list[,1]),] #对分组数据进行排序,按照数据框的第一列升序排序
    group<-new_group[,2] #取第二列,这样处理得到的就是每一个位置所对应的样本属性信息
    

    接着,加载limma包,构造如下矩阵。

    suppressMessages(library(limma))
    design <- model.matrix(~0+factor(group))
    colnames(design)=levels(factor(group))
    rownames(design)=colnames(data)
    

    3.构建两组间的比较矩阵
    这段的代码的意思就是比较design矩阵中HC与MDD level(这里为什么只对比了两组,不是说对比多组吗?是的,当然可以对比多组(这个后面会说如何对比多组),但是对比多组后,结果得到的P值只有一个,也就是说结果里的P值是多组对比的综合P值,你不能从这个结果里获取任意两组的P值
    因此,如果大家需要获取任何两组的P值以及FC的话,个人建议每组比较单独做一个topTable

    contrast.matrix<-makeContrasts("HC-MDD",levels=design)
    

    接着跑limma的三个标准步骤,就可以得到所要的结果了

    ##step1
    fit <- lmFit(data,design)
    ##step2
    fit2 <- contrasts.fit(fit, contrast.matrix)
    fit2 <- eBayes(fit2)  
    ##step3
    tempOutput = topTable(fit2, coef=1, n=Inf)
    nrDEG = na.omit(tempOutput) 
    

    4.构建多组间的比较矩阵

    ##两组以上的比较,寻找差异基因,使用topTableF函数
    cont.wt <- makeContrasts("F12h.1-Fcontrol",##举例
                             "F24h.1-F12h.1",
                             "F5d.1-F24h.1",
                             "F10d.1-F5d.1",
                             levels=design) 
    
    fit <- lmFit(eset, design)
    fit2 <- contrasts.fit(fit, cont.wt) 
    fit2 <- eBayes(fit2) 
    tT=topTableF(fit2, adjust="BH",sort.by="F",n=Inf)
    tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) #这一条和下一条代码需要根据自己需要进行更改,topTableF函数得到的结果中包括F值,而logFC则是没两组的比较,均得出了
    colnames(tT)=c("FDR","P.Value","logFC")
    
    1. 保存分析好的差异表达数据。
    write.table(nrDEG,"DEG_GSE98793.txt")
    write.table(tT ,"tT.txt")
    

    两种方法的全部代码:

    两组表型数据差异分析

    data<-read.table("GSE98793.txt")
    group_list<-read.table("gene_list.txt")
    new_group<- group_list[order(group_list[,1]),] 
    group<-new_group[,2]
     
    suppressMessages(library(limma))
     
     
    design <- model.matrix(~0+factor(group))
    colnames(design)=levels(factor(group))
    rownames(design)=colnames(data)
     
    contrast.matrix<-makeContrasts("HC-MDD",levels=design)
     
    ##step1
    fit <- lmFit(data,design)
    ##step2
    fit2 <- contrasts.fit(fit, contrast.matrix)
    fit2 <- eBayes(fit2) 
    ##step3
    tempOutput = topTable(fit2, coef=1, n=Inf)
    nrDEG = na.omit(tempOutput) 
    write.table(nrDEG,"limma_GSE98793.txt")
    

    多组表型数据差异分析

    library(limma)##加载包
    load("GSE37761exp_groupfile-00.RData")#加载匹配过的基因表达矩阵
    eset=exp ##将基因表达矩阵赋值给eset
    targets<-read.csv("group_file.csv",row.names = 1)##读入样本数据,包括两列,
    #第一列GSM号,第二列为样本分组#该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则
    targets$Target=gsub("_",".",targets$Target)##将"_"替换成“.”,也可以不替换
    ##该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则,所以在没个分类前加一个“F”,具体自己定
    targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))
    colnames(targets)=c("FileName","Target")#更改列名,为了和limma包中的一致
    lev<-unique(targets$Target)##使用unique()函数进行去重
    f <- factor(targets$Target, levels=lev) 
    design <- model.matrix(~0+f) #样本矩阵
    colnames(design) <- lev #更改列名为levels名
    ##两组以上的比较,寻找差异基因,使用topTableF函数
    cont.wt <- makeContrasts("F12h.1-Fcontrol",##举例
                             "F24h.1-F12h.1",
                             "F5d.1-F24h.1",
                             "F10d.1-F5d.1",
                             levels=design) 
    
    fit <- lmFit(eset, design)
    fit2 <- contrasts.fit(fit, cont.wt) 
    fit2 <- eBayes(fit2) 
    tT=topTableF(fit2, adjust="BH",sort.by="F",n=Inf)
    ##tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) 这一条和下一条代码需要根据自己需要进行更改,topTableF函数得到的结果中包括F值,而logFC则是没两组的比较,均得出了
    ##colnames(tT)=c("FDR","P.Value","logFC")
    

    用limma包进行多组差异表达分析
    R语言limma包差异基因分析

    相关文章

      网友评论

        本文标题:对于多表型数据如何利用limma包进行差异基因分析

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