setwd("E:/GSE25066")#環(huán)境設(shè)置
library(limma)#加載差異分析包limma
#將分組文件加載到環(huán)境中,分組信息第一列為樣本名横朋,第二列為分組信息如“high”“l(fā)ow”
targets<-read.csv("group.csv")
#將表達(dá)矩陣加載到環(huán)境中,行為基因炕置,列為樣本迫摔,這里應(yīng)該注意去除重復(fù)項(xiàng)。
eset<-read.csv("expreset-basal1.csv",row.names = "symbol")
targets$Target=gsub("_",".",targets$Target)##若數(shù)據(jù)中存在特殊符號症脂,將"_"替換成“.”,也可以不替換
##該數(shù)據(jù)集中實(shí)際存在不符合R的命名原則淫僻,所以在沒個分類前加一個“F”诱篷,具體自己定
targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))
colnames(targets)=c("FileName","Target")#更改列名,為了和limma包中的一致
lev<-unique(targets$Target)##使用unique()函數(shù)進(jìn)行去重
f <- factor(targets$Target, levels=lev)
design <- model.matrix(~0+f)
colnames(design) <- lev
cont.wt <- makeContrasts("high-low",
? ? ? ? ? ? ? ? ? ? ? +? levels=design)
fit <- lmFit(eset, design)#前面矩陣的row.name=“symbol”
fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
tT=topTable(fit2, adjust="BH",sort.by="logFC",n=Inf)
tT = subset(tT, select=c("adj.P.Val","P.Value","logFC"))
colnames(tT)=c("FDR","P.Value","logFC")
write.csv(tT,"DEGbasal.csv")
#最后的tT就是得到的差異基因雳灵,其中包含基因棕所,P.Value和logFC