Hello, sorry for the probably silly question, but I've tried for days now. I'm a beginner. I have a list of genomic coordinates in GenomicRanges format, e.g.

> coords_GR
GRanges object with 3320 ranges and 0 metadata columns:
         seqnames            ranges strand
            <Rle>         <IRanges>  <Rle>
     [1]        1           1950293      *
     [2]        1 18336405-20142656      *
     [3]        1          29684252      *
     [4]        1 71218722-73861224      *
     [5]        1 80888506-83526065      *
     ...      ...               ...    ...
  [3316]       20    716956-2488855      *
  [3317]       21 41901419-44757190      *
  [3318]       22 29255810-31043932      *
  [3319]       22 35134992-38911889      *
  [3320]        X         154717327      *
  -------
  seqinfo: 24 sequences from an unspecified genome; no seqlengths

And I would like to annotate all SNVs in the range, obtaining information such as RSid, MAF, and, if exonic, if it's synonimous/non-synonimous. Also it would be good to be able to annotate non-SNP variations such as indels, CNVs, deletions etc for potential protein coding/clinical impact. I found a number of tools that can do that for genotyping data in VCF (such as ensembl VEP), but no tools that can do it based on genomic coordinates alone.

Many thanks for your help

Emanuele

Hi,

As a first step you could use a SNPlocs package to map your genomic ranges to rs ids:

library(SNPlocs.Hsapiens.dbSNP144.GRCh38)
snps <- SNPlocs.Hsapiens.dbSNP144.GRCh38
my_regions <- GRanges(c("1:1950293", "1:18336405-20142656", "1:29684252",
                        "1:71218722-73861224", "22:35134992-38911889"))
my_snps <- snpsByOverlaps(snps, my_regions)
my_snps
# UnstitchedGPos object with 411257 positions and 2 metadata columns:
#            seqnames       pos strand |   RefSNP_id alleles_as_ambig
#               <Rle> <integer>  <Rle> | <character>      <character>
#        [1]        1   1950293      * | rs145840463                R
#        [2]        1  18336405      * | rs147541038                R
#        [3]        1  18336413      * | rs542452980                R
#        [4]        1  18336416      * | rs575136982                R
#        [5]        1  18336432      * | rs182869403                K
#        ...      ...       ...    ... .         ...              ...
#   [411253]       22  38911800      * | rs535959691                R
#   [411254]       22  38911824      * | rs555183244                R
#   [411255]       22  38911852      * | rs538293072                K
#   [411256]       22  38911861      * | rs189740977                Y
#   [411257]       22  38911889      * | rs181225505                Y
#   -------
#   seqinfo: 25 sequences (1 circular) from GRCh38.p2 genome

Then use your preferred tool for retrieving details about the rs ids.

Notes:

1. Unfortunately the regions of origin of the SNPs are not reported but you can easily retrieve them with:

   hits <- findOverlaps(my_snps, my_regions)
   region_idx <- as(hits, "IntegerList")
   mcols(my_snps)$region_idx <- region_idx
   mcols(my_snps)$region <- extractList(my_regions, region_idx)
   my_snps
   # UnstitchedGPos object with 411257 positions and 4 metadata columns:
   #            seqnames       pos strand |   RefSNP_id alleles_as_ambig
   #               <Rle> <integer>  <Rle> | <character>      <character>
   #        [1]        1   1950293      * | rs145840463                R
   #        [2]        1  18336405      * | rs147541038                R
   #        [3]        1  18336413      * | rs542452980                R
   #        [4]        1  18336416      * | rs575136982                R
   #        [5]        1  18336432      * | rs182869403                K
   #        ...      ...       ...    ... .         ...              ...
   #   [411253]       22  38911800      * | rs535959691                R
   #   [411254]       22  38911824      * | rs555183244                R
   #   [411255]       22  38911852      * | rs538293072                K
   #   [411256]       22  38911861      * | rs189740977                Y
   #   [411257]       22  38911889      * | rs181225505                Y
   #               region_idx               region
   #            <IntegerList>        <GRangesList>
   #        [1]             1            1:1950293
   #        [2]             2  1:18336405-20142656
   #        [3]             2  1:18336405-20142656
   #        [4]             2  1:18336405-20142656
   #        [5]             2  1:18336405-20142656
   #        ...           ...                  ...
   #   [411253]             5 22:35134992-38911889
   #   [411254]             5 22:35134992-38911889
   #   [411255]             5 22:35134992-38911889
   #   [411256]             5 22:35134992-38911889
   #   [411257]             5 22:35134992-38911889
   #   -------
   #   seqinfo: 25 sequences (1 circular) from GRCh38.p2 genome
   

2. If you want SNVs in addition to SNPs, do the same as the above with the XtraSNPlocs.Hsapiens.dbSNP144.GRCh38 package.

See ?SNPlocs and ?XtraSNPlocs in the BSgenome package for more information.

Best,

H.