R:limma表达谱分析

添加第一列列名为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:///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)