limma package计算差异表达

test.txt:文件格式：

gene	TCGA-3M-AB46-01A-11R-A414-31	TCGA-3M-AB47-01A-22R-A414-31	TCGA-B7-A5TK-01A-12R-A36D-31	TCGA-B7-A5TN-01A-21R-A31P-31
A2LD1	431	606	2731	194
A2M	35864	15942	22086	109363
AADACL2	0	2	10	0
AADAT	643	388	317	241

pheno.txt文件格式：

TCGA-3M-AB46-01A-11R-A414-31 a
TCGA-3M-AB47-01A-22R-A414-31 b
TCGA-B7-A5TK-01A-12R-A36D-31 b
TCGA-B7-A5TN-01A-21R-A31P-31 a
TCGA-BR-6457-01A-21R-1802-13 a
TCGA-BR-6458-01A-11R-1802-13 a
TCGA-BR-6706-01A-11R-1884-13 b

## Librarys
library(limma)
## Matrix File
raw.data <-read.delim("test.txt")
attach(raw.data)
names(raw.data)
d <- raw.data[, 2:130]
rownames(d) <- raw.data[, 1]
# Pheno data file
pheno <- read.table("pheno.txt",header = F)
colnames(pheno) <- c("sample","group")
Group<-factor(pheno$group,levels=levels(pheno$group))
##To design matrix---
design<-model.matrix(~0+Group)
##Normalisation 
y <- voom(d,design,plot=TRUE)
colnames(design)
fit <-lmFit(y,design)
##Designing Contrast Matrix for group Differentiation
cont.wt<-makeContrasts(Groupb-Groupa,levels=design)
fit2 <-contrasts.fit(fit,cont.wt)
fit3<-eBayes(fit2)
DE<-topTable(fit3,number = 'all')
write.csv(DE, "limma_genes_test.csv")
selected = rownames(DE[DE$P.Value <0.05,])
# 默认是holm方法控制FWER
esetSel = y[selected, ]
heatmap(esetSel@.Data[[1]])
#p<- DE[DE$P.Value <1,]$P.Value 
#esetSela <- y[]
#p.adjust(c(0.001,0.02,0.03),"BH")

 

转载于:https://www.cnblogs.com/qiniqnyang/p/6564054.html