我用beadarry和limma包对GSE49454数据集进行差异基因表达分析时,最后卡在这个group name上面了,尝试了好多方法都不行
代码:library(GEOquery)
library(beadarray)
library(illuminaHumanv4.db)
#下载数据
url<-"ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE49nnn/GSE49454/matrix/"
filenm<-"GSE49454_series_matrix.txt.gz"
if(!file.exists("GSE49454_series_matrix.txt.gz")) download.file(paste(url, filenm, sep=""), destfile=filenm)
gse <- getGEO(filename=filenm)
#将下载的数据转换为ExpressionSetIllumina,并注释
summaryData <- as(gse, "ExpressionSetIllumina")
rna <- factor(pData(summaryData)[,"characteristics_ch1"])
#去除非匹配
fData(summaryData)$Status <-
ifelse(fData(summaryData)$PROBEQUALITY=="No match","negative","regular" )
Detection(summaryData) <- calculateDetection(summaryData,
status=fData(summaryData)$Status)
#normalization
summaryData.norm <- normaliseIllumina(summaryData,method="neqc",
status=fData(summaryData)$Status)
group <- pData(summaryData.norm)[ ,"characteristics_ch1"]
limmaRes <- limmaDE(summaryData, SampleGroup="characteristics_ch1")
design <- model.matrix(~0+rna)
design
colnames(design) <- levels(rna)
aw <- arrayWeights(exprs(summaryData.norm), design)
aw
fit <- lmFit(exprs(summaryData.norm), design, weights=aw)
contrasts <- makeContrasts(group: SLE-group: Healthy control of SLE, levels=design)
contr.fit <- eBayes(contrasts.fit(fit, contrasts))
topTable(contr.fit, coef=1)
报错:
contrasts <- makeContrasts(SLE-control, levels=design)
Error in makeContrasts(SLE - control, levels = design) :
The levels must by syntactically valid names in R, see help(make.names). Non-valid names: rnagroup: Healthy control of SLE,rnagroup: SLE
limmaRes <- limmaDE(summaryData, SampleGroup="characteristics_ch1")
Calculating array weights
Array weights
Error in makeContrasts(contrasts = contrast, levels = design) :
The levels must by syntactically valid names in R, see help(make.names). Non-valid names: group: Healthy control of SLE,group: SLE