Using GEOquery and SAM/Siggenes in R

I very recently started learning R as a result of a need and so far, so good, I开发者_如何学Python think. But I'm still in the very early stages. I am however faced with this major urgent challenge in R that I will greatly appreciate some help with. My programming skills are quite obviously amateur and will most definitely accept any help I can get. Here goes:

1. To create a list of datasets (gdslist) to be retrieved from the GEO database using the GEOquery package
2. To convert the gdslist items (gdsid) to expression data. That is, data that can work with my analysis. For this, a GDS2eSet function does the work fine.
3. To read in this converted expression data in a way that a class/levels file (.cls) can be created. The GDS3715 dataset for example, has 3 levels – insulin resistant, insulin sensitive and diabetic. Sometimes, the datasets are as straightforward as that. But other times, like in this case, the levels will be 6 for analysis purposes, because although there are phenotypically 3 levels, they have been divided into treated and untreated groups. There is often an added “agent” column in such cases. Each class/level needs to be assigned a numerical number (0,1,2...). That’s pretty much the general format of a .cls file.
4. To run a Siggenes/SAM analysis (also a package in R), two files are needed for each dataset: an expression file (the converted file from 2 above) and the accompanying cluster file (from 3).
5. To be able to run this process for the gdslist items in a kind of loop and have my data stored in a specified directory.

I can currently only get to step 2. I think step 3 is the crux of the challenge...

Many thanks in anticipation.

Script so far:

> 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 

---

This should work for all the gds datasets.

 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)