Valerie Obenchain | 19 Jun 00:11 2013

Re: SNP annotation in R

Hi Ina,

It looks like the chromosome names in 'target' don't match those in the 
TranscriptDb. There is an example on the locateVariants man page of how 
to inspect the seqlevels (chromosome names) with seqlevels() and change 
them with renameSeqlevels().

Once you've made the change, confirm that the names match with the 
intersect() function. This is also demonstrated on the man page.

Valerie

On 06/18/2013 01:47 PM, Ina Hoeschele wrote:
> Hi,
>    I am interested in annotating lists of SNPs in R, mostly I am interested in finding the gene a SNP is located
in and for intergenic SNPs which genes are close by. I found this code that wsa developed with the help of
Valerie Obenchain:
> http://adairama.wordpress.com/2013/02/15/functionally-annotate-snps-and-indels-in-bioconductor/
> However, when I adapted this code to my problem I got the error message below. Can someone please give me a hint?
> Thanks, Ina
>
> input <- SNP.map.retain[SNPs,1:3]
> rownames(input) <- NULL
> colnames(input) <- c("rsid", "chr", "pos")
> input
> #       rsid chr      pos
> #1  rs661761   1 30618672
> #2 rs3754508   1 17255474
> #3 rs6426403   1  4169394
> #4  rs228727   1  7770423
> #5 rs1635592   1 17535994
> target <- with(input,
>                  GRanges( seqnames = Rle(chr),
>                           ranges   = IRanges(pos, end=pos, names=rsid),
>                           strand   = Rle(strand("*")) ) )
> target
> #GRanges with 5 ranges and 0 metadata columns:
> #            seqnames               ranges strand
> #               <Rle>            <IRanges>  <Rle>
> #   rs661761        1 [30618672, 30618672]      *
> #  rs3754508        1 [17255474, 17255474]      *
> #  rs6426403        1 [ 4169394,  4169394]      *
> #   rs228727        1 [ 7770423,  7770423]      *
> #  rs1635592        1 [17535994, 17535994]      *
> #  ---
> #  seqlengths:
> #    1
> #   NA
>
> loc <- locateVariants(target, TxDb.Hsapiens.UCSC.hg19.knownGene, AllVariants())
> Error in .mk_isActiveSeqReplacementValue(x, value) :
>    the names of the supplied 'isActiveSeq' must match the names of the current 'isActiveSeq'
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@...
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>

_______________________________________________
Bioconductor mailing list
Bioconductor@...
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor


Gmane