### R code from vignette source 'SRAdb.Rnw' ################################################### ### code chunk number 1: init ################################################### options(width=50) ################################################### ### code chunk number 2: SRAdb.Rnw:56-59 ################################################### library(SRAdb) sqlfile <- 'SRAmetadb.sqlite' if(!file.exists('SRAmetadb.sqlite')) sqlfile <<- getSRAdbFile() ################################################### ### code chunk number 3: SRAdb.Rnw:64-67 (eval = FALSE) ################################################### ## timeStart <- proc.time() ## sqlfile <- getSRAdbFile() ## proc.time() - timeStart ################################################### ### code chunk number 4: SRAdb.Rnw:72-73 ################################################### file.info('SRAmetadb.sqlite') ################################################### ### code chunk number 5: SRAdb.Rnw:78-79 ################################################### sra_con <- dbConnect(SQLite(),sqlfile) ################################################### ### code chunk number 6: SRAdb.Rnw:89-91 ################################################### sra_tables <- dbListTables(sra_con) sra_tables ################################################### ### code chunk number 7: SRAdb.Rnw:95-96 ################################################### dbListFields(sra_con,"study") ################################################### ### code chunk number 8: SRAdb.Rnw:100-101 ################################################### dbGetQuery(sra_con,'PRAGMA TABLE_INFO(study)') ################################################### ### code chunk number 9: SRAdb.Rnw:105-107 ################################################### colDesc <- colDescriptions(sra_con=sra_con)[1:5,] colDesc[, 1:4] ################################################### ### code chunk number 10: j1 ################################################### rs <- dbGetQuery(sra_con,"select * from study limit 3") rs[, 1:3] ################################################### ### code chunk number 11: 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 12: SRAdb.Rnw:128-134 ################################################### getTableCounts <- function(tableName,conn) { sql <- sprintf("select count(*) from %s",tableName) return(dbGetQuery(conn,sql)[1,1]) } do.call(rbind,sapply(sra_tables[c(2,4,5,11,12)], getTableCounts, sra_con, simplify=FALSE)) ################################################### ### code chunk number 13: SRAdb.Rnw:139-143 ################################################### rs <- dbGetQuery(sra_con, paste( "SELECT study_type AS StudyType, count( * ) AS Number FROM `study` GROUP BY study_type order by Number DESC ", sep="")) rs ################################################### ### code chunk number 14: SRAdb.Rnw:148-152 ################################################### rs <- dbGetQuery(sra_con, paste( "SELECT instrument_model AS 'Instrument Model', count( * ) AS Experiments FROM `experiment` GROUP BY instrument_model order by Experiments DESC", sep="")) rs ################################################### ### code chunk number 15: SRAdb.Rnw:156-160 ################################################### rs <- dbGetQuery(sra_con, paste( "SELECT library_strategy AS 'Library Strategy', count( * ) AS Runs FROM `experiment` GROUP BY library_strategy order by Runs DESC", sep="")) rs ################################################### ### code chunk number 16: SRAdb.Rnw:168-170 ################################################### conversion <- sraConvert( c('SRP001007','SRP000931'), sra_con = sra_con ) conversion[1:3,] ################################################### ### code chunk number 17: SRAdb.Rnw:174-175 ################################################### apply(conversion, 2, unique) ################################################### ### code chunk number 18: SRAdb.Rnw:183-196 ################################################### rs <- getSRA( search_terms = "breast cancer", out_types = c('run','study'), sra_con ) dim(rs) rs <- getSRA( search_terms = "breast cancer", out_types = c("submission", "study", "sample", "experiment", "run"), sra_con ) # get counts for some information interested apply( rs[, c('run','sample','study_type','platform', 'instrument_model')], 2, function(x) {length(unique(x))} ) ################################################### ### code chunk number 19: SRAdb.Rnw:200-203 ################################################### rs <- getSRA (search_terms ='"breast cancer"', out_types=c('run','study'), sra_con) dim(rs) ################################################### ### code chunk number 20: SRAdb.Rnw:207-210 ################################################### rs <- getSRA( search_terms ='MCF7 OR "MCF-7"', out_types = c('sample'), sra_con ) dim(rs) ################################################### ### code chunk number 21: SRAdb.Rnw:214-217 ################################################### rs <- getSRA( search_terms ='submission_center: GEO', out_types = c('submission'), sra_con ) dim(rs) ################################################### ### code chunk number 22: SRAdb.Rnw:221-224 ################################################### rs <- getSRA( search_terms ='Carcino*', out_types = c('study'), sra_con=sra_con ) dim(rs) ################################################### ### code chunk number 23: SRAdb.Rnw:231-232 ################################################### rs = listSRAfile( c("SRX000122"), sra_con, fileType = 'sra' ) ################################################### ### code chunk number 24: SRAdb.Rnw:237-239 ################################################### rs = getSRAinfo ( c("SRX000122"), sra_con, sraType = "sra" ) rs[1:3,] ################################################### ### code chunk number 25: SRAdb.Rnw:244-245 ################################################### getSRAfile( c("SRR000648","SRR000657"), sra_con, fileType = 'sra' ) ################################################### ### code chunk number 26: SRAdb.Rnw:251-252 ################################################### ## system ("fastq-dump SRR000648.lite.sra") ################################################### ### code chunk number 27: SRAdb.Rnw:256-259 (eval = FALSE) ################################################### ## getFASTQinfo( c("SRR000648","SRR000657"), sra_con, srcType = 'ftp' ) ## ## getSRAfile( c("SRR000648","SRR000657"), sra_con, fileType = 'fastq' ) ################################################### ### code chunk number 28: SRAdb.Rnw:266-286 (eval = FALSE) ################################################### ## ## List fasp addresses for associated fastq files: ## listSRAfile ( c("SRX000122"), sra_con, fileType = 'fastq', srcType='fasp') ## ## ## get fasp addresses for associated fastq files: ## getFASTQinfo( c("SRX000122"), sra_con, srcType = 'fasp' ) ## ## ## download fastq files using fasp protocol: ## # the following ascpCMD needs to be constructed according custom ## # system configuration ## # common ascp installation in a Linux system: ## ascpCMD <- 'ascp -QT -l 300m -i ## /usr/local/aspera/connect/etc/asperaweb_id_dsa.putty' ## ## ## common ascpCMD for a Mac OS X system: ## # ascpCMD <- "'/Applications/Aspera Connect.app/Contents/ ## # Resources/ascp' -QT -l 300m -i '/Applications/ ## # Aspera Connect.app/Contents/Resources/asperaweb_id_dsa.putty'" ## ## getSRAfile( c("SRX000122"), sra_con, fileType = 'fastq', ## srcType = 'fasp', ascpCMD = ascpCMD ) ################################################### ### code chunk number 29: SRAdb.Rnw:290-296 (eval = FALSE) ################################################### ## ## List fasp addresses of sra files associated with "SRX000122" ## listSRAfile( c("SRX000122"), sra_con, fileType = 'sra', srcType='fasp') ## ## ## download sra files using fasp protocol ## getSRAfile( c("SRX000122"), sra_con, fileType = 'sra', ## srcType = 'fasp', ascpCMD = ascpCMD ) ################################################### ### code chunk number 30: SRAdb.Rnw:311-312 (eval = FALSE) ################################################### ## startIGV("mm") ################################################### ### code chunk number 31: SRAdb.Rnw:316-323 (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 32: SRAdb.Rnw:338-350 (eval = FALSE) ################################################### ## library(SRAdb) ## library(Rgraphviz) ## ## g <- sraGraph('primary thyroid cell line', sra_con) ## attrs <- getDefaultAttrs(list(node=list( ## fillcolor='lightblue', shape='ellipse'))) ## plot(g, attrs=attrs) ## ## ## similiar search as the above, returned much larger data.frame and graph is too clouded ## g <- sraGraph('Ewing Sarcoma', sra_con) ## plot(g) ## ################################################### ### code chunk number 33: SRAdb.Rnw:357-358 ################################################### dbDisconnect(sra_con) ################################################### ### code chunk number 34: SRAdb.Rnw:370-389 (eval = FALSE) ################################################### ## ## library(SRAdb) ## ## setwd('1000g') ## if( ! file.exists('SRAmetadb.sqlite') ) { ## sqlfile <- getSRAdbFile() ## } else { ## sqlfile <- 'SRAmetadb.sqlite' ## } ## sra_con <- dbConnect(SQLite(),sqlfile) ## ## ## get all related accessions ## rs <- getSRA( search_terms = '"1000 Genomes Project"', ## sra_con=sra_con, acc_only=TRUE) ## dim(rs) ## head(rs) ## ## ## get counts for each data types ## apply( rs, 2, function(x) {length(unique(x))} ) ################################################### ### code chunk number 35: SRAdb.Rnw:394-396 (eval = FALSE) ################################################### ## runs <- tail(rs$run) ## fs <- getSRAinfo( runs, sra_con, sraType = "sra" ) ################################################### ### code chunk number 36: SRAdb.Rnw:400-401 (eval = FALSE) ################################################### ## getSRAfile( runs, sra_con, fileType ='sra', srcType = "ftp" ) ################################################### ### code chunk number 37: SRAdb.Rnw:405-408 (eval = FALSE) ################################################### ## ascpCMD <- "'/Applications/Aspera Connect.app/Contents/Resources/ascp' -QT -l 300m -i '/Applications/Aspera Connect.app/Contents/Resources/asperaweb_id_dsa.putty'" ## ## sra_files = getSRAfile( runs, sra_con, fileType ='sra', srcType = "fasp", ascpCMD = ascpCMD ) ################################################### ### code chunk number 38: SRAdb.Rnw:412-415 (eval = FALSE) ################################################### ## for( fq in basename(sra_files$fasp) ) { ## system ("fastq-dump SRR000648.lite.sra") ## } ################################################### ### code chunk number 39: SRAdb.Rnw:423-424 ################################################### toLatex(sessionInfo())