添加第一列列名為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("
> 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)
