Thanks James, Valerie, and Steve.

I was tripped up by the distance() function in BioC 2.11, which requires
both arguments to have same length, but upgrading to 2.12 solved that.

Cheers,
Gad




On 14 May 2013 02:23, Steve Lianoglou <lianoglou.steve@gene.com> wrote:

> Hi,
>
> On Mon, May 13, 2013 at 8:48 AM, James W. MacDonald <jmacdon@uw.edu>
> wrote:
> > Hi Gad,
> >
> >
> > On 5/12/2013 11:41 PM, Gad Abraham wrote:
> >>
> >> Hi,
> >>
> >> I'm using VariantAnnotation, which is a very nice package, to annotate
> >> some
> >> SNPs, based on
> >>
> >>
> http://adairama.wordpress.com/2013/02/15/functionally-annotate-snps-and-indels-in-bioconductor
> >> .
> >>
> >> For SNPs annotated as intergenic and therefore falling between two genes
> >> (PRECEDEID and FOLLOWID), is there a way of getting the actual distance
> >> between the SNPs and these two genes?
> >>
> >> I'm using:
> >>
> >> head(m, 3)
> >>           RS  Chr      Pos
> >> 1 rs2187668 chr6 32713862
> >> 2 rs9357152 chr6 32772938
> >> 3 rs3099844 chr6 31556955
> >>
> >> target<- with(m,
> >>        GRanges(seqnames=Rle(Chr),
> >>           ranges=IRanges(Pos, end=Pos, names=RS),
> >>           strand=Rle(strand("*"))
> >>        )
> >> )
> >>
> >> loc<- locateVariants(target, TxDb.Hsapiens.UCSC.hg18.knownGene,
> >> AllVariants())
> >> names(loc)<- NULL
> >> out<- as.data.frame(loc)
> >>
> >> head(out, 3)
> >>        names seqnames    start      end   LOCATION GENEID PRECEDEID
> >> FOLLOWID
> >> GeneSymbol PrecedeSymbol FollowSymbol
> >> 1 rs2187668     chr6 32713862 32713862     intron   3117<NA>
> >> <NA>    HLA-DQA1<NA>          <NA>
> >> 4 rs9357152     chr6 32772938 32772938 intergenic<NA>       3119
> >> 3117<NA>       HLA-DQB1     HLA-DQA1
> >> 5 rs3099844     chr6 31556955 31556955 intergenic<NA>       4277
> >> 352961<NA>           MICB        HCG26
> >>
> >> I'd like to get the distance from rs9357152 to HLA-DQB1 and to HLA-DQA1.
> >
> >
> > There might be a far superior way to do this, but one way is to assume
> that
> > cds == genes and go from there:
> >
> >> gns <- cds(TxDb.Hsapiens.UCSC.hg19.knownGene,
> >> columns=c("cds_id","gene_id"))
> >> ind <- sapply(values(gns)$gene_id, function(x) any(x %in%
> >> c("3117","3119")))
> >> sum(ind)
> > [1] 28
> >> gns[ind,]
> > GRanges with 28 ranges and 2 metadata columns:
> >             seqnames               ranges strand   |    cds_id
> > gene_id
> > <Rle> <IRanges> <Rle>   | <integer> <CharacterList>
> >    [1]          chr6 [32605236, 32605317]      +   |     73476
> > 3117
> >    [2]          chr6 [32609087, 32609335]      +   |     73477
> > 3117
> >    [3]          chr6 [32609749, 32610030]      +   |     73478
> > 3117
> >    [4]          chr6 [32609749, 32610164]      +   |     73479
> > 3117
> >    [5]          chr6 [32610387, 32610541]      +   |     73480
> > 3117
> >    ...           ...                  ...    ... ...       ...
> > ...
> >   [24] chr6_qbl_hap6   [3860872, 3860982]      -   |    233860
> > 3119
> >   [25] chr6_qbl_hap6   [3861407, 3861780]      -   |    233861
> > 3119
> >   [26] chr6_qbl_hap6   [3861499, 3861780]      -   |    233862
> > 3119
> >   [27] chr6_qbl_hap6   [3864670, 3864939]      -   |    233863
> > 3119
> >   [28] chr6_qbl_hap6   [3866378, 3866486]      -   |    233864
> > 3119
> >
> > Note that these two genes also map to the extra Chr6 haplotypes.
> >
> >> z <- cbind(as.data.frame(gns[ind,]), distance(GRanges("chr6",
> >> IRanges(start=32772938, end=32772938)), gns[ind,]))[,c(1,2,4,6,7,8)]
> > Warning message:
> > In .local(x, y, ...) :
> >   The behavior of distance() has changed in Bioconductor 2.12. See
> ?distance
> > for details.
> >> colnames(z) <- c(colnames(z)[1:5], "Distance")
>
> Saving one line of code makes it twice as nice for half the price:
>
> R> z <- cbind(as.data.frame(gns[ind,]),
> Distance=distance(GRanges("chr6", IRanges(start=32772938,
> end=32772938)), gns[ind,]))[,c(1,2,4,6,7,8)]
>
> R> head(z)
>    seqnames    start width cds_id gene_id Distance
> 1|     chr6 32605236    82  73476    3117   167620
> 2|     chr6 32609087   249  73477    3117   163602
> 3|     chr6 32609749   282  73478    3117   162907
> 4|     chr6 32609749   416  73479    3117   162773
> 5|     chr6 32610387   155  73480    3117   162396
> 6|     chr6 32628013    14  78928    3119   144911
>
> This time, it's free ... this time.
>
> -steve
>
> --
> Steve Lianoglou
> Computational Biologist
> Department of Bioinformatics and Computational Biology
> Genentech
>

	[[alternative HTML version deleted]]

