添加第一列列名為id

清空空字符

文件保存為csv格式

#表達(dá)矩陣

>exprSet=read.csv("PC表達(dá)差異.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）的火山圖

補(bǔ)充:關(guān)于limma包差異分析結(jié)果的logFC解釋

首先，我們要明白，limma接受的輸入?yún)?shù)就是一個表達(dá)矩陣，而且是log后的表達(dá)矩陣（以2為底）。

那么最后計算得到的logFC這一列的值，其實就是輸入的表達(dá)矩陣中case一組的平均表達(dá)量減去control一組的平均表達(dá)量的值，那么就會有正負(fù)之分，代表了case相當(dāng)于control組來說，該基因是上調(diào)還是下調(diào)。

我之前總是有疑問，明明是case一組的平均表達(dá)量和control一組的平均表達(dá)量差值呀，跟log foldchange沒有什么關(guān)系呀。

后來，我終于想通了，因為我們輸入的是log后的表達(dá)矩陣，那么case一組的平均表達(dá)量和control一組的平均表達(dá)量都是log了的，那么它們的差值其實就是log的foldchange

首先，我們要理解foldchange的意義，如果case是平均表達(dá)量是8，control是2，那么foldchange就是4，logFC就是2咯

那么在limma包里面，輸入的時候case的平均表達(dá)量被log后是3，control是1，那么差值是2，就是說logFC就是2。

這不是巧合，只是一個很簡單的數(shù)學(xué)公式log(x/y)=log(x)-log(y)