Hi Dario,

Yep, I did that with snp135common; the really annoying part is that you
need to look in the .sql file for the appropriate table's schema to find
the column names, and (at least in my case) I had to use zcat and 'cut' to
get the .txt columns that were present for all entries so that I could
convert it to a GRanges.  Buit after that it was fine.  read.table,
names(), and a df2GR() and resize() later, we seem to be OK:

$ wget
ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp135Common.txt.gz
$ zcat snp135Common.txt.gz | cut -f2,3,4,5,7,14 | gzip >
snp135Common.BioC.txt.gz
$ R
## Welcome to R version 2.15.0 (2012-03-30)
## BiocInstaller version 1.4.3, ?biocLite for help
snp135common.df <- read.table("snp135Common.BioC.txt.gz")
## time passes...
names(snp135common.df) <-
c('chrom','chromStart','chromEnd','name','strand','score')
source("df2GR.R")
snp135common.GR <- df2GR(snp135common.df, keepColumns=TRUE)
show(snp135common.GR)
## GRanges with 11525489 ranges and 1 elementMetadata col:
##              seqnames               ranges strand   |     score
##                 <Rle>            <IRanges>  <Rle>   | <numeric>
##   rs55998931     chr1       [10491, 10492]      +   |  0.115787
##   rs58108140     chr1       [10582, 10583]      +   |  0.394721
##   rs75454623     chr1       [14929, 14930]      +   |      0.48
##   rs71262674     chr1       [20244, 20245]      -   |       0.5
##   rs71262673     chr1       [20303, 20304]      -   |       0.5
## ...
## 'score' is average heterozygosity, in case some with MAF < 0.01 sneak in
## loading snp135.txt (cut properly) instead of snp135common would bring
those in, along with indels, etc.
## ...
genome(snp135common.GR)<-'hg19'

## Now compare to Herve's ranges:
##
library(SNPlocs.Hsapiens.dbSNP.20111119)
rsidsToGRanges(c('rs55998931','rs71262673'))
## GRanges with 2 ranges and 2 elementMetadata cols:
##      seqnames         ranges strand |   RefSNP_id alleles_as_ambig
##         <Rle>      <IRanges>  <Rle> | <character>      <character>
##  [1]      ch1 [10492, 10492]      * |    55998931                Y
##  [2]      ch1 [20304, 20304]      * |    71262673                S
## ...

snp135common.GR[c('rs55998931','rs71262673')]
## GRanges with 2 ranges and 1 elementMetadata col:
##              seqnames         ranges strand |     score
##                 <Rle>      <IRanges>  <Rle> | <numeric>
##   rs55998931     chr1 [10491, 10492]      + |  0.115787
##   rs71262673     chr1 [20303, 20304]      - |       0.5
## ...
ranges(snp135common.GR[c('rs55998931','rs71262673')])==ranges(rsidsToGRanges(c('rs55998931','rs71262673')))
## [1] FALSE FALSE

## UCSC ranges are too wide!  And the stranding is funny.
##
strand(snp135common.GR) <- '*'
snp135common.GR <- resize(snp135common.GR, width=1, fix='end')
ranges(snp135common.GR[c('rs55998931','rs71262673')])==ranges(rsidsToGRanges(c('rs55998931','rs71262673')))
## [1] TRUE TRUE

save(snp135common.GR, file="snp135common.GR.rda")
## in lieu of a proper FeatureDb, for now
##
quit('n')

$ cat df2GR.R
##
require(GenomicRanges)
##
## df2GR: the stupid man's FeatureDb!
##
df2GR <- function(df, keepColumns=FALSE, ignoreStrand=FALSE) { # orig. from
KDH
  stopifnot(class(df) == "data.frame")
  subs <- list(chromStart='start', chromEnd='end', chrom='chr',
seqnames='chr')
  for(s in names(subs)) names(df) = gsub(s, subs[[s]], names(df),
ignore=TRUE)
  stopifnot(all(c("start", "end") %in% names(df)))
  if('genome' %in% names(attributes(df))) g <- attr(df, 'genome') else g <-
NULL
  if(substr(df$chr, 1, 3)[1] != 'chr') df$chr <- paste('chr', df$chr,
sep='')
  df <- subset(df, !is.na(start) & !is.na(end))
  if(!ignoreStrand && ("strand" %in% names(df))) {
    if(is.numeric(df$strand)) df$strand <- strandMe(df$strand)
    GR <- with(df, GRanges(chr, IRanges(start=start, end=end),
strand=strand))
  } else {
    GR <- with(df, GRanges(chr, IRanges(start=start, end=end)))
  }
  if('name' %in% names(df)) {
    names(GR) <- df$name
    df$name <- NULL
  } else {
    names(GR) <- rownames(df)
  }
  if(keepColumns) {
    skipped = c("rangename","chr","start","end","width","strand")
    elementMetadata(GR) <- as(df[, setdiff(names(df), skipped), drop=F],
                              "DataFrame")
  }
  if('X' %in% names(elementMetadata(GR))) {
    if(all(is.na(GR$X))) {
      GR$X <- NULL
    } else {

names(elementMetadata(GR))[which(names(elementMetadata(GR))=='X')]='score'
    }
  }
  if(!is.null(g)) genome(GR) <- g
  return(GR)
}

Hope this helps someone (particularly Moiz ;-))



On Mon, May 7, 2012 at 5:00 PM, Dario Strbenac <D.Strbenac@garvan.org.au>wrote:

> I e-mailed UCSC and they said the preferred way is to download by FTP.
> Which means more lines of code to parse the text file into columns, then
> split the exons and widths columns up to be able to make GRanges.
>
> ---- Original message ----
> >Date: Sun, 6 May 2012 22:36:53 -0700
> >From: "Tim Triche, Jr." <tim.triche@gmail.com>
> >Subject: Re: [Bioc-devel] GenomicFeatures FeatureDB EST Table Error
> >To: D.Strbenac@garvan.org.au
> >Cc: bioc-devel@r-project.org
> >
> >   Actually, try downloading the same thing from the
> >   Table Browser and see if there isn't something like
> >   the following at the tail of the file:
> >   843  chr11  33910774  33910775  rs4756078  0  +  CC
> >   C/G/T  genomic  single
> >   by-cluster,by-frequency,by-2hit-2allele,by-hapmap,by-1000genomes
> >   0.361204  0.223906  intron  exact  1
> >
> SingleClassTriA---------------------------------------------------------------------------
> >   procedures have exceeded timeout: 1200 seconds,
> >   function has ended.
> >
> ---------------------------------------------------------------------------
> >   (this is from my attempted download of the
> >   snp135common track, but it appears to be happening
> >   to you as well)
> >   It would appear that we're being throttled.
> >   On Sun, May 6, 2012 at 10:00 PM, Dario Strbenac
> >   <D.Strbenac@garvan.org.au> wrote:
> >
> >     Hi,
> >
> >     It seems some data has been added to the EST data
> >     table in UCSC that GenomicFeatures cannot parse.
> >
> >     > ESTs <- makeFeatureDbFromUCSC(genome = "hg19",
> >     track = "est", tablename = "all_est")
> >     Download the all_est table ... Error in scan(file,
> >     what, nmax, sep, dec, quote, skip, nlines,
> >     na.strings,  :
> >      line 4587019 did not have 22 elements
> >
> >     The relevant lines of the session information are
> >     :
> >
> >     R version 2.15.0 (2012-03-30)
> >     Platform: x86_64-unknown-linux-gnu (64-bit)
> >
> >     other attached packages:
> >     [1] GenomicFeatures_1.8.1
> >
> >     - Dario.
> >
> >     _______________________________________________
> >     Bioc-devel@r-project.org mailing list
> >     https://stat.ethz.ch/mailman/listinfo/bioc-devel
> >
> >   --
> >   A model is a lie that helps you see the truth.
> >   Howard Skipper
>
>
> --------------------------------------
> Dario Strbenac
> Research Assistant
> Cancer Epigenetics
> Garvan Institute of Medical Research
> Darlinghurst NSW 2010
> Australia
>
> _______________________________________________
> Bioc-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/bioc-devel
>



-- 
*A model is a lie that helps you see the truth.*
*
*
Howard Skipper<http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>

	[[alternative HTML version deleted]]

