我用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
https://www.jianshu.com/p/5f94ae79f298
简单地说,这个地方需要的是 pData(summaryData.norm)
的一列作为分组变量,你用的 characteristics_ch1
这一列,但是这一列的内容为 “Healthy control of SLE“ 和 ”group: SLE”,这个不是 R 里合法的变量名(里面既有空格也有引号)。所以你应该给数据的 pData 自己加上一列用来在做差异分析的时候分组用。