Использование GEOquery и SAM / Siggenes в R - PullRequest
1 голос
/ 12 сентября 2011

Я совсем недавно начал изучать R в результате необходимости, и пока, думаю, все хорошо. Но я все еще на самой ранней стадии. Однако я сталкиваюсь с этой серьезной проблемой в R, с которой я буду очень признателен за помощь. Мои навыки программирования, безусловно, любительские и, безусловно, будут принимать любую помощь, которую я могу получить. Здесь идет:

  1. Для создания списка наборов данных (gdslist), которые будут получены из GEO база данных с использованием пакета GEOquery
  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 

1 Ответ

0 голосов
/ 21 декабря 2011

Это должно работать для всех наборов данных 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)
...