添加第一列列名为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)
网友评论