[BioC] SNP annotation in R

[Valerie Obenchain]

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 at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>