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 回答