美文网首页
R:limma表达谱分析

R:limma表达谱分析

作者: 小米羊爱学术 | 来源:发表于2019-03-06 10:00 被阅读0次

    添加第一列列名为id

    清空空字符

    文件保存为csv格式

    #表达矩阵

    >exprSet=read.csv("PC表达差异.csv",header = T,row.names = "id")

    > exprSet[1:5,1:5]

                               N1        N2       N3        N4       N5

    ASCRP000001 7.724977  7.867182 7.716320  7.808999 7.769155

    ASCRP000002 9.358083  9.419830 9.366810 10.023439 9.854187

    ASCRP000004 9.376542 12.278848 9.372987  9.013128 9.864421

    ASCRP000005 8.326384  8.292985 8.341654  8.242107 8.327296

    ASCRP000006 8.056148  9.256802 8.053611  8.311279 8.318445

    #分组

    > group_list = c(rep("normal",6),rep("cancer",6))

    > group_list

    [1] "normal" "normal" "normal" "normal" "normal" "normal" "cancer" "cancer" "cancer" "cancer" "cancer" "cancer"

    #QC检测

    > par(cex=0.7)

    > n.sample = ncol(exprSet)

    > cols=rainbow(n.sample*1.2)

    > boxplot(exprSet, col = cols,main="expression value",las=2)

    #安装limma

    > source("https://bioconductor.org/biocLite.R")

    > biocLite("limma")

    > suppressMessages(library(limma))

    #制作分组矩阵

    >design <- model.matrix(~0+factor(group_list))

    >colnames(design)=levels(factor(group_list))

    >rownames(design)=colnames(exprSet)

    > design

       cancer normal

    N1      0      1

    N2      0      1

    N3      0      1

    N4      0      1

    N5      0      1

    N6      0      1

    C1      1      0

    C2      1      0

    C3      1      0

    C4      1      0

    C5      1      0

    C6      1      0

    #矩阵声明

    > contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

    > contrast.matrix

            Contrasts

    Levels   normal-cancer

    cancer            -1

    normal             1

    #差异分析

    > fit <- lmFit(exprSet,design)

    > fit2 <- contrasts.fit(fit, contrast.matrix)

    > fit2 <- eBayes(fit2)  ## default no trend !!!

    > tempOutput = topTable(fit2, coef=1, n=Inf)

    > nrDEG = na.omit(tempOutput) 

    > head(nrDEG)

                     logFC   AveExpr         t      P.Value  adj.P.Val         B

    ASCRP002607 -0.6921963  8.304751 -6.171921 6.396350e-05 0.06249354 1.9542295

    ASCRP002273 -0.5178892  8.374084 -5.872071 9.897068e-05 0.06249354 1.5815667

    ASCRP004643 -0.3363617  8.378703 -5.805455 1.092339e-04 0.06249354 1.4966747

    ASCRP000624  0.4530361 15.559685  5.634824 1.410412e-04 0.06249354 1.2757264

    ASCRP001621 -0.3981851  8.270907 -5.389916 2.050045e-04 0.06249354 0.9497525

    ASCRP003058 -1.7550218 12.010918 -5.344000 2.201025e-04 0.06249354 0.8874756

    > write.csv(nrDEG,"limma_notrend.results.csv",quote = F)

    最后附上logFC和-log(P.Value)的火山图

    补充:关于limma包差异分析结果的logFC解释

    首先,我们要明白,limma接受的输入参数就是一个表达矩阵,而且是log后的表达矩阵(以2为底)。

    那么最后计算得到的logFC这一列的值,其实就是输入的表达矩阵中case一组的平均表达量减去control一组的平均表达量的值,那么就会有正负之分,代表了case相当于control组来说,该基因是上调还是下调。

    我之前总是有疑问,明明是case一组的平均表达量和control一组的平均表达量差值呀,跟log foldchange没有什么关系呀。

    后来,我终于想通了,因为我们输入的是log后的表达矩阵,那么case一组的平均表达量和control一组的平均表达量都是log了的,那么它们的差值其实就是log的foldchange

    首先,我们要理解foldchange的意义,如果case是平均表达量是8,control是2,那么foldchange就是4,logFC就是2咯

    那么在limma包里面,输入的时候case的平均表达量被log后是3,control是1,那么差值是2,就是说logFC就是2。

    这不是巧合,只是一个很简单的数学公式log(x/y)=log(x)-log(y)

    相关文章

      网友评论

          本文标题:R:limma表达谱分析

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