### R code from vignette source '/Users/sdavis/Documents/git/publicDataTutorial/inst/doc/publicDataTutorial.Rnw' ################################################### ### code chunk number 1: loadgeoquery ################################################### library(GEOquery) ################################################### ### code chunk number 2: testgeoquery ################################################### # If you have network access, the more typical way to do this # would be to use this: # gds <- getGEO("GDS507") gds <- getGEO(filename=system.file("extdata/GDS507.soft.gz",package="GEOquery")) ################################################### ### code chunk number 3: publicDataTutorial.Rnw:76-80 ################################################### # If you have network access, the more typical way to do this # would be to use this: # gds <- getGEO("GSM11805") gsm <- getGEO(filename=system.file("extdata/GSM11805.txt.gz",package="GEOquery")) ################################################### ### code chunk number 4: metagsm ################################################### # Look at gsm metadata: Meta(gsm) ################################################### ### code chunk number 5: tablegsm ################################################### # Look at data associated with the GSM: # but restrict to only first 5 rows, for brevity Table(gsm)[1:5,] ################################################### ### code chunk number 6: columngsm ################################################### # Look at Column descriptions: Columns(gsm) ################################################### ### code chunk number 7: columnsgds ################################################### Columns(gds) ################################################### ### code chunk number 8: publicDataTutorial.Rnw:128-139 ################################################### # Again, with good network access, one would do: # gse <- getGEO("GSE781",GSEMatrix=FALSE) gse <- getGEO(filename=system.file("extdata/GSE781_family.soft.gz",package="GEOquery")) head(Meta(gse)) # names of all the GSM objects contained in the GSE names(GSMList(gse)) # and get the first GSM object on the list class(GSMList(gse)[[1]]) head(Meta(GSMList(gse)[[1]])) # and the names of the GPLs represented names(GPLList(gse)) ################################################### ### code chunk number 9: gsematrix ################################################### # Note that GSEMatrix=TRUE is the default gse2553 <- getGEO('GSE2553',GSEMatrix=TRUE) show(gse2553) show(pData(phenoData(gse2553[[1]]))[1:5,c(1,6,8)]) ################################################### ### code chunk number 10: gdseset ################################################### eset <- GDS2eSet(gds,do.log2=TRUE) ################################################### ### code chunk number 11: eset ################################################### eset pData(eset) ################################################### ### code chunk number 12: gdsmalist1 ################################################### #get the platform from the GDS metadata Meta(gds)$platform #So use this information in a call to getGEO gpl <- getGEO(filename=system.file("extdata/GPL97.annot.gz",package="GEOquery")) ################################################### ### code chunk number 13: gdsmalist2 ################################################### MA <- GDS2MA(gds,GPL=gpl) MA ################################################### ### code chunk number 14: gseconvert1 ################################################### gsmplatforms <- lapply(GSMList(gse),function(x) {Meta(x)$platform}) gsmplatforms ################################################### ### code chunk number 15: gseconvert2 ################################################### Table(GSMList(gse)[[1]])[1:5,] # and get the column descriptions Columns(GSMList(gse)[[1]])[1:5,] ################################################### ### code chunk number 16: gseconvert3 ################################################### # get the probeset ordering probesets <- Table(GPLList(gse)[[1]])$ID # make the data matrix from the VALUE columns from each GSM # being careful to match the order of the probesets in the platform # with those in the GSMs data.matrix <- do.call('cbind',lapply(GSMList(gse),function(x) {tab <- Table(x) mymatch <- match(probesets,tab$ID_REF) return(tab$VALUE[mymatch]) })) data.matrix <- apply(data.matrix,2,function(x) {as.numeric(as.character(x))}) data.matrix <- log2(data.matrix) data.matrix[1:5,] ################################################### ### code chunk number 17: gseconvert4 ################################################### require(Biobase) # go through the necessary steps to make a compliant ExpressionSet rownames(data.matrix) <- probesets colnames(data.matrix) <- names(GSMList(gse)) pdata <- data.frame(samples=names(GSMList(gse))) rownames(pdata) <- names(GSMList(gse)) pheno <- as(pdata,"AnnotatedDataFrame") eset2 <- new('ExpressionSet',exprs=data.matrix,phenoData=pheno) eset2 ################################################### ### code chunk number 18: getgeosuppfiles ################################################### df = getGEOSuppFiles('GSM3922') ################################################### ### code chunk number 19: getgeosuppfilesreturn ################################################### df ################################################### ### code chunk number 20: publicDataTutorial.Rnw:279-285 ################################################### gpl97 <- getGEO('GPL97') Meta(gpl97)$title head(Meta(gpl97)$series_id) length(Meta(gpl97)$series_id) head(Meta(gpl97)$sample_id) length(Meta(gpl97)$sample_id) ################################################### ### code chunk number 21: publicDataTutorial.Rnw:290-295 ################################################### gsmids <- Meta(gpl97)$sample_id # Feel free to run the next two lines, but I leave them out # here to cut down on processing time # gsmlist <- sapply(gsmids[1:5],getGEO) # names(gsmlist) ################################################### ### code chunk number 22: geoquerymirror ################################################### destdir = tempdir() # this will be downloaded x = getGEO('GDS507',destdir=destdir) ################################################### ### code chunk number 23: geoquerymirror2 ################################################### # this will NOT be downloaded # local copy will be used instead y = getGEO('GDS507',destdir=destdir) ################################################### ### code chunk number 24: publicDataTutorial.Rnw:363-364 ################################################### options(width=55) ################################################### ### code chunk number 25: publicDataTutorial.Rnw:371-372 ################################################### library(GEOmetadb) ################################################### ### code chunk number 26: publicDataTutorial.Rnw:376-379 ################################################### if(!file.exists('GEOmetadb.sqlite')) { getSQLiteFile() } ################################################### ### code chunk number 27: publicDataTutorial.Rnw:386-387 ################################################### file.info('GEOmetadb.sqlite') ################################################### ### code chunk number 28: publicDataTutorial.Rnw:391-393 ################################################### con <- dbConnect(SQLite(),'GEOmetadb.sqlite') dbDisconnect(con) ################################################### ### code chunk number 29: publicDataTutorial.Rnw:406-407 ################################################### con <- dbConnect(SQLite(),'GEOmetadb.sqlite') ################################################### ### code chunk number 30: publicDataTutorial.Rnw:411-413 ################################################### geo_tables <- dbListTables(con) geo_tables ################################################### ### code chunk number 31: publicDataTutorial.Rnw:417-418 ################################################### dbListFields(con,'gse') ################################################### ### code chunk number 32: publicDataTutorial.Rnw:422-423 ################################################### sqliteQuickSQL(con,'PRAGMA TABLE_INFO(gpl)') ################################################### ### code chunk number 33: publicDataTutorial.Rnw:429-431 ################################################### rs <- dbGetQuery(con,'select * from gse limit 5') rs[,1:7] ################################################### ### code chunk number 34: publicDataTutorial.Rnw:435-438 ################################################### rs <- dbGetQuery(con,paste("select gse,title from gse where", "contributor like '%Sean%Davis%'",sep=" ")) rs ################################################### ### code chunk number 35: publicDataTutorial.Rnw:443-447 ################################################### rs <- dbGetQuery(con,paste("select gsm,supplementary_file", "from gsm where gpl='GPL96'", "and supplementary_file like '%CEL.gz'")) dim(rs) ################################################### ### code chunk number 36: publicDataTutorial.Rnw:452-458 ################################################### rs <- dbGetQuery(con,paste("select gpl.bioc_package,gsm.gpl,", "gsm,gsm.supplementary_file", "from gsm join gpl on gsm.gpl=gpl.gpl", "where gpl.manufacturer='Affymetrix'", "and gsm.supplementary_file like '%CEL.gz' ")) rs[1:5,] ################################################### ### code chunk number 37: publicDataTutorial.Rnw:462-467 ################################################### getTableCounts <- function(tableName,conn) { sql <- sprintf("select count(*) from %s",tableName) return(dbGetQuery(conn,sql)[1,1]) } do.call(rbind,sapply(geo_tables,getTableCounts,con,simplify=FALSE)) ################################################### ### code chunk number 38: publicDataTutorial.Rnw:476-477 ################################################### conversion <- geoConvert('GPL96') ################################################### ### code chunk number 39: publicDataTutorial.Rnw:482-487 ################################################### lapply(conversion, dim) conversion$gse[1:5,] conversion$gsm[1:5,] conversion$gds[1:5,] conversion$sMatrix[1:5,] ################################################### ### code chunk number 40: publicDataTutorial.Rnw:494-495 ################################################### getBiocPlatformMap(con, bioc=c('hgu133a','hgu95av2')) ################################################### ### code chunk number 41: publicDataTutorial.Rnw:501-514 ################################################### sql <- paste("SELECT DISTINCT gse.title,gse.gse", "FROM", " gsm JOIN gse_gsm ON gsm.gsm=gse_gsm.gsm", " JOIN gse ON gse_gsm.gse=gse.gse", " JOIN gse_gpl ON gse_gpl.gse=gse.gse", " JOIN gpl ON gse_gpl.gpl=gpl.gpl", "WHERE", " gsm.molecule_ch1 like '%total RNA%' AND", " gse.title LIKE '%breast cancer%' AND", " gpl.organism LIKE '%Homo sapiens%'",sep=" ") rs <- dbGetQuery(con,sql) dim(rs) print(rs[1:5,],right=FALSE) ################################################### ### code chunk number 42: publicDataTutorial.Rnw:519-520 ################################################### dbDisconnect(con) ################################################### ### code chunk number 43: publicDataTutorial.Rnw:524-525 (eval = FALSE) ################################################### ## file.remove('GEOmetadb.sqlite') ################################################### ### code chunk number 44: publicDataTutorial.Rnw:540-544 ################################################### library(SRAdb) if(!file.exists('SRAdb.sqlite')) { sqlfile <- getSRAdbFile() } ################################################### ### code chunk number 45: publicDataTutorial.Rnw:549-550 ################################################### file.info('SRAmetadb.sqlite') ################################################### ### code chunk number 46: publicDataTutorial.Rnw:555-558 ################################################### #open connection sra_con <- dbConnect(SQLite(),'SRAdb.sqlite') sra_con ################################################### ### code chunk number 47: publicDataTutorial.Rnw:568-570 ################################################### sra_tables <- dbListTables(sra_con) sra_tables ################################################### ### code chunk number 48: publicDataTutorial.Rnw:574-575 ################################################### dbListFields(sra_con,'study') ################################################### ### code chunk number 49: publicDataTutorial.Rnw:579-580 ################################################### sqliteQuickSQL(sra_con,'PRAGMA TABLE_INFO(study)') ################################################### ### code chunk number 50: j1 ################################################### rs <- dbGetQuery(sra_con,'select * from study limit 3') rs[,1:5] ################################################### ### code chunk number 51: j2 ################################################### rs <- dbGetQuery(sra_con,paste("select study_accession,study_title from study where", "study_description like 'Transcriptome%'",sep=" ")) rs[1:3,] ################################################### ### code chunk number 52: publicDataTutorial.Rnw:600-605 ################################################### getTableCounts <- function(tableName,conn) { sql <- sprintf("select count(*) from %s",tableName) return(dbGetQuery(conn,sql)[1,1]) } do.call(rbind,sapply(sra_tables,getTableCounts,sra_con,simplify=FALSE)) ################################################### ### code chunk number 53: publicDataTutorial.Rnw:614-616 ################################################### conversion <- sraConvert(c('SRP001007','SRP000931'), sra_con= sra_con) conversion[1:3,] ################################################### ### code chunk number 54: publicDataTutorial.Rnw:620-621 ################################################### apply(conversion, 2, unique) ################################################### ### code chunk number 55: publicDataTutorial.Rnw:629-631 ################################################### rs <- getSRA (search_terms ='breast cancer', out_types=c('run','study'), sra_con=sra_con) dim(rs) ################################################### ### code chunk number 56: publicDataTutorial.Rnw:635-637 ################################################### rs <- getSRA (search_terms ='"breast cancer"', out_types=c('run','study'), sra_con=sra_con) dim(rs) ################################################### ### code chunk number 57: publicDataTutorial.Rnw:641-643 ################################################### rs <- getSRA (search_terms ='MCF7 OR "MCF-7"', out_types=c('sample'), sra_con=sra_con) dim(rs) ################################################### ### code chunk number 58: publicDataTutorial.Rnw:647-649 ################################################### rs <- getSRA (search_terms ='submission_center: GEO', out_types=c('submission'), sra_con=sra_con) dim(rs) ################################################### ### code chunk number 59: publicDataTutorial.Rnw:653-655 ################################################### rs <- getSRA (search_terms ='Carcino*', out_types=c('study'), sra_con=sra_con) dim(rs) ################################################### ### code chunk number 60: publicDataTutorial.Rnw:661-662 ################################################### listSRAfile ("SRX000122", sra_con=sra_con, sraType='litesra') ################################################### ### code chunk number 61: publicDataTutorial.Rnw:666-668 ################################################### rs <- getSRAinfo (in_acc=c("SRX000122"), sra_con=sra_con) rs[1:3,] ################################################### ### code chunk number 62: publicDataTutorial.Rnw:672-673 (eval = FALSE) ################################################### ## getSRAfile (in_acc=c("SRR000648","SRR000657"), sra_con=sra_con, destdir=getwd(), sraType='litesra') ################################################### ### code chunk number 63: publicDataTutorial.Rnw:684-685 (eval = FALSE) ################################################### ## startIGV("mm") ################################################### ### code chunk number 64: publicDataTutorial.Rnw:689-696 (eval = FALSE) ################################################### ## exampleBams = file.path(system.file('extdata',package='SRAdb'), ## dir(system.file('extdata',package='SRAdb'),pattern='bam$')) ## sock <- IGVsocket() ## IGVgenome(sock, 'hg18') ## IGVload(sock, exampleBams) ## IGVgoto(sock, 'chr1:1-1000') ## IGVsnapshot(sock) ################################################### ### code chunk number 65: publicDataTutorial.Rnw:708-713 ################################################### library(Rgraphviz) acc <- getSRA (search_terms ='colon cancer', out_types=c('sra'), sra_con=sra_con, acc_only=TRUE) g <- entityGraph(acc) attrs <- getDefaultAttrs(list(node=list(fillcolor='lightblue', shape='ellipse'))) plot(g, attrs= attrs) ################################################### ### code chunk number 66: publicDataTutorial.Rnw:717-720 ################################################### g <- sraGraph('colon cancer', sra_con) attrs <- getDefaultAttrs(list(node=list(fillcolor='lightblue', shape='ellipse'))) plot(g, attrs=attrs) ################################################### ### code chunk number 67: publicDataTutorial.Rnw:725-726 ################################################### dbDisconnect(sra_con) ################################################### ### code chunk number 68: publicDataTutorial.Rnw:732-733 ################################################### toLatex(sessionInfo())