1

由于需要,我最近开始学习 R,到目前为止,我认为非常好。但我仍处于早期阶段。然而,我在 R 中面临着这个重大的紧迫挑战,我将非常感谢一些帮助。我的编程技能显然是业余的,并且肯定会接受我能得到的任何帮助。开始:

  1. 使用 GEOquery 包创建要从 GEO 数据库中检索的数据集列表 (gdslist)
  2. 将 gdslist 项目 (gdsid) 转换为表达式数据。也就是说,可以与我的分析一起使用的数据。为此,GDS2eSet 函数可以正常工作。
  3. 以可以创建类/级别文件 (.cls) 的方式读取此转换后的表达式数据。例如,GDS3715 数据集有 3 个级别——胰岛素抵抗、胰岛素敏感和糖尿病。有时,数据集就这么简单。但在其他时候,例如在这种情况下,出于分析目的,水平将为 6,因为虽然表型上有 3 个水平,但它们已被分为治疗组和未治疗组。在这种情况下,通常会添加一个“代理”列。每个班级/级别都需要分配一个数字编号(0,1,2...)。这几乎是 .cls 文件的一般格式。
  4. 要运行 Siggenes/SAM 分析(也是 R 中的一个包),每个数据集需要两个文件:一个表达式文件(从上面 2 转换的文件)和随附的集群文件(从 3)。
  5. 为了能够在一种循环中为 gdslist 项目运行此过程,并将我的数据存储在指定的目录中。

我目前只能进入第 2 步。我认为第 3 步是挑战的关键……

非常感谢期待。

到目前为止的脚本:

> gdslist = c('GDS3715','GDS3716','GDS3717'...)#up to perhaps 100 datasets
> analysisfunc = function(gdsid) {
    gdsdat = getGEO(gdsid,destdir=".")
    gdseset = GDS2eSet(gdsdat)
    pData(gdseset)$disease.state #Needed assignment, etc...Step 3 stuff ;Siggenes/SAM can perhaps be done here
    return(sprintf("Results from %s should be here",gdsid))
  }
> resultlist = sapply(gdslist,analysisfunc) #loop function 
4

1 回答 1

0

这应该适用于所有 gds 数据集。

 GEOSAM.analysis <- function( gdsid, destdir = getwd() ) {
       require( 'GEOquery' )
       require( 'siggenes' )
       ## test if gdsid is gdsid
       if( length(grep('GDS', gdsid)) == 0 ){
        stop()
       }
       gdsdat = getGEO( gdsid, destdir = destdir )
       gdseset = GDS2eSet( gdsdat )
       gdseset.pData <- pData( gdseset )
       gds.factors <- names( gdseset.pData )
       gds.factors[gds.factors == 'sample'] <- NA
       gds.factors[gds.factors == 'description'] <- NA
       gds.factors <- gds.factors[!is.na( gds.factors )]
       cl.list <- sapply( gdseset.pData[gds.factors], as.character)
       cl.list <- factor( apply( cl.list, 1, function(x){ paste( x , collapse = '-' )} ) )
       if( length( levels ( cl.list ) ) == 2 ){
        levels( cl.list ) <- 0:length( levels( cl.list ) )
       } else {
        levels( cl.list ) <- 1:length( levels( cl.list ) )
       }
       sam.gds <- sam( gdseset, cl.list )
       results.file <- file.path( destdir, paste( gdsid, '.sam.gds.rdata', sep =''  ) )
       save( sam.gds, file = results.file )
       return( sprintf( "Results from %s are saved in '%s'. These can be loaded by 'load('%s')'.",gdsid, results.file, results.file ) )
  }

  gdslist = c('GDS3715', 'GDS3716', 'GDS3717')
  resultlist = sapply(gdslist, GEOSAM.analysis)  
  print(resultlist)
于 2011-12-21T15:30:38.813 回答