Annotation of Genomic Regions
1
0
Entering edit mode
elcombe.c • 0
@elcombec-23949
Last seen 3.7 years ago

I have a list of locations on the genome, which I have converted to GRanges object, and now wish to annotate.

Sample data:

 head(CpG_Ranges)
GRanges object with 6 ranges and 10 metadata columns:
      seqnames              ranges strand | n_samples      pvalue
         <Rle>           <IRanges>  <Rle> | <integer>   <numeric>
  [1]     chr1 110821406-110821412      * |         2 0.000253685
  [2]     chr1 110820767-110820776      * |         2 0.000887079
  [3]    chr14     1862428-1862429      * |         2 0.001082515
  [4]    chr24   40883339-40883354      * |         2 0.001207324
  [5]    chr14   49361102-49361103      * |         2 0.001705802
  [6]    chr14   56132052-56132053      * |         2 0.001876952

I have created a txdb file, and found and fetched the appropriate annotation package using Annotation hub as below:

gffFile <- "G:\\oviAri4.ncbiRefSeq.gtf"
txdb <- makeTxDbFromGFF(file=gffFile, format=c("gtf"))

ah <- AnnotationHub()
query(ah, c("Ovis aries", "OrgDb"))
org.Oa.eg.db <- ah[["AH80684"]]

Then I'm using annotatePeak as such:

peakAnno <- annotatePeak(CpG_Ranges, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Oa.eg.db")

Which runs fine for binning the data into genomic regions (exon, intron, TSS, etc), but fails in the annotation step:

>> preparing features information...         2020-08-05 12:09:06 
>> identifying nearest features...       2020-08-05 12:09:06 
>> calculating distance from peak to TSS...  2020-08-05 12:09:07 
>> assigning genomic annotation...       2020-08-05 12:09:07 
>> adding gene annotation...             2020-08-05 12:09:09 
>> assigning chromosome lengths          2020-08-05 12:09:09 
>> done...                   2020-08-05 12:09:09 
Warning messages:
1: In .Seqinfo.mergexy(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': chrUn_NW_014639147v1, chrUn_NW_014639227v1, chrUn_NW_014639425v1, chrUn_NW_014639454v1, chrUn_NW_014639514v1, chrUn_NW_014639584v1, chrUn_NW_014639757v1, chrUn_NW_014639941v1, chrUn_NW_014640025v1, chrUn_NW_014640441v1, chrUn_NW_014640601v1, chrUn_NW_014640620v1, chrUn_NW_014640625v1, chrUn_NW_014640781v1, chrUn_NW_014640783v1, chrUn_NW_014640831v1, chrUn_NW_014640847v1, chrUn_NW_014641029v1, chrUn_NW_014641066v1, chrUn_NW_014641250v1, chrUn_NW_014641439v1, chrUn_NW_014641483v1, chrUn_NW_014641578v1, chrUn_NW_014641740v1, chrUn_NW_014641910v1, chrUn_NW_014642173v1, chrUn_NW_014642238v1, chrUn_NW_014642275v1
  - in 'y': chrM, chrUn_NW_014639041v1, chrUn_NW_014639042v1, chrUn_NW_014639057v1, chrUn_NW_014639070v1, chrUn_NW_014639071v1, chrUn_NW_014639073v1, chrUn_NW_014639074v1, chrUn_NW_014639084v1, chrUn_NW_014639092v1, chrUn_NW_014639127v1, chrUn_NW_014639132v1, chrUn_NW_014639140v1, chrUn_NW_014639151v1,  [... truncated]
2: In .Seqinfo.mergexy(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': chrUn_NW_014639147v1, chrUn_NW_014639227v1, chrUn_NW_014639425v1, chrUn_NW_014639454v1, chrUn_NW_014639514v1, chrUn_NW_014639584v1, chrUn_NW_014639757v1, chrUn_NW_014639941v1, chrUn_NW_014640025v1, chrUn_NW_014640441v1, chrUn_NW_014640601v1, chrUn_NW_014640620v1, chrUn_NW_014640625v1, chrUn_NW_014640781v1, chrUn_NW_014640783v1, chrUn_NW_014640831v1, chrUn_NW_014640847v1, chrUn_NW_014641029v1, chrUn_NW_014641066v1, chrUn_NW_014641250v1, chrUn_NW_014641439v1, chrUn_NW_014641483v1, chrUn_NW_014641578v1, chrUn_NW_014641740v1, chrUn_NW_014641910v1, chrUn_NW_014642173v1, chrUn_NW_014642238v1, chrUn_NW_014642275v1
  - in 'y': chrM, chrUn_NW_014639041v1, chrUn_NW_014639042v1, chrUn_NW_014639057v1, chrUn_NW_014639070v1, chrUn_NW_014639071v1, chrUn_NW_014639073v1, chrUn_NW_014639074v1, chrUn_NW_014639084v1, chrUn_NW_014639092v1, chrUn_NW_014639127v1, chrUn_NW_014639132v1, chrUn_NW_014639140v1, chrUn_NW_014639151v1,  [... truncated]
3: In .Seqinfo.mergexy(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': chrUn_NW_014639147v1, chrUn_NW_014639227v1, chrUn_NW_014639425v1, chrUn_NW_014639454v1, chrUn_NW_014639514v1, chrUn_NW_014639584v1, chrUn_NW_014639757v1, chrUn_NW_014639941v1, chrUn_NW_014640025v1, chrUn_NW_014640441v1, chrUn_NW_014640601v1, chrUn_NW_014640620v1, chrUn_NW_014640625v1, chrUn_NW_014640781v1, chrUn_NW_014640783v1, chrUn_NW_014640831v1, chrUn_NW_014640847v1, chrUn_NW_014641029v1, chrUn_NW_014641066v1, chrUn_NW_014641250v1, chrUn_NW_014641439v1, chrUn_NW_014641483v1, chrUn_NW_014641578v1, chrUn_NW_014641740v1, chrUn_NW_014641910v1, chrUn_NW_014642173v1, chrUn_NW_014642238v1, chrUn_NW_014642275v1
  - in 'y': chrM, chrUn_NW_014639041v1, chrUn_NW_014639042v1, chrUn_NW_014639057v1, chrUn_NW_014639070v1, chrUn_NW_014639071v1, chrUn_NW_014639073v1, chrUn_NW_014639074v1, chrUn_NW_014639084v1, chrUn_NW_014639092v1, chrUn_NW_014639127v1, chrUn_NW_014639132v1, chrUn_NW_014639140v1, chrUn_NW_014639151v1,  [... truncated]
4: In .Seqinfo.mergexy(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': chrUn_NW_014639147v1, chrUn_NW_014639227v1, chrUn_NW_014639425v1, chrUn_NW_014639454v1, chrUn_NW_014639514v1, chrUn_NW_014639584v1, chrUn_NW_014639757v1, chrUn_NW_014639941v1, chrUn_NW_014640025v1, chrUn_NW_014640441v1, chrUn_NW_014640601v1, chrUn_NW_014640620v1, chrUn_NW_014640625v1, chrUn_NW_014640781v1, chrUn_NW_014640783v1, chrUn_NW_014640831v1, chrUn_NW_014640847v1, chrUn_NW_014641029v1, chrUn_NW_014641066v1, chrUn_NW_014641250v1, chrUn_NW_014641439v1, chrUn_NW_014641483v1, chrUn_NW_014641578v1, chrUn_NW_014641740v1, chrUn_NW_014641910v1, chrUn_NW_014642173v1, chrUn_NW_014642238v1, chrUn_NW_014642275v1
  - in 'y': chrM, chrUn_NW_014639041v1, chrUn_NW_014639042v1, chrUn_NW_014639057v1, chrUn_NW_014639070v1, chrUn_NW_014639071v1, chrUn_NW_014639073v1, chrUn_NW_014639074v1, chrUn_NW_014639084v1, chrUn_NW_014639092v1, chrUn_NW_014639127v1, chrUn_NW_014639132v1, chrUn_NW_014639140v1, chrUn_NW_014639151v1,  [... truncated]
5: In .Seqinfo.mergexy(x, y) :
  Each of the 2 combined objects has sequence levels not in the other:
  - in 'x': chrUn_NW_014639147v1, chrUn_NW_014639227v1, chrUn_NW_014639425v1, chrUn_NW_014639454v1, chrUn_NW_014639514v1, chrUn_NW_014639584v1, chrUn_NW_014639757v1, chrUn_NW_014639941v1, chrUn_NW_014640025v1, chrUn_NW_014640441v1, chrUn_NW_014640601v1, chrUn_NW_014640620v1, chrUn_NW_014640625v1, chrUn_NW_014640781v1, chrUn_NW_014640783v1, chrUn_NW_014640831v1, chrUn_NW_014640847v1, chrUn_NW_014641029v1, chrUn_NW_014641066v1, chrUn_NW_014641250v1, chrUn_NW_014641439v1, chrUn_NW_014641483v1, chrUn_NW_014641578v1, chrUn_NW_014641740v1, chrUn_NW_014641910v1, chrUn_NW_014642173v1, chrUn_NW_014642238v1, chrUn_NW_014642275v1
  - in 'y': chrM, chrUn_NW_014639041v1, chrUn_NW_014639042v1, chrUn_NW_014639057v1, chrUn_NW_014639070v1, chrUn_NW_014639071v1, chrUn_NW_014639073v1, chrUn_NW_014639074v1, chrUn_NW_014639084v1, chrUn_NW_014639092v1, chrUn_NW_014639127v1, chrUn_NW_014639132v1, chrUn_NW_014639140v1, chrUn_NW_014639151v1,  [... truncated]
6: In annotatePeak(CpG_Ranges, tssRegion = c(-3000, 3000), TxDb = txdb,  :
  Unknown ID type, gene annotation will not be added...

I could manually find which genes these are located in/near, but really need to automate this.

Any help would be greatly appreciated!
annotation genomic regions R • 813 views
ADD COMMENT •
2
Entering edit mode
@james-w-macdonald-5106
Last seen 4 hours ago
United States

Your problem arises from the GTF file you used. It's not ideal because the gene id column contains gene symbols instead of, like, Gene IDs.

> z <- import("oviAri4.ncbiRefSeq.gtf.gz")
> z
GRanges object with 1364113 ranges and 9 metadata columns:
            seqnames              ranges strand |     source        type
               <Rle>           <IRanges>  <Rle> |   <factor>    <factor>
        [1]     chr1         49237-54549      + | ncbiRefSeq  transcript
        [2]     chr1         49237-49423      + | ncbiRefSeq  exon      
        [3]     chr1         52671-52819      + | ncbiRefSeq  exon      
        [4]     chr1         54488-54549      + | ncbiRefSeq  exon      
        [5]     chr1         72197-75389      - | ncbiRefSeq  transcript
        ...      ...                 ...    ... .        ...         ...
  [1364109]     chrX 134863497-134863637      - | ncbiRefSeq 5UTR       
  [1364110]     chrX 134864183-134864433      - | ncbiRefSeq exon       
  [1364111]     chrX 134864183-134864433      - | ncbiRefSeq 5UTR       
  [1364112]     chrX 134863494-134863496      - | ncbiRefSeq start_codon
  [1364113]     chrX 134862348-134862350      - | ncbiRefSeq stop_codon 
                score     phase      gene_id  transcript_id    gene_name
            <numeric> <integer>  <character>    <character>  <character>
        [1]        NA      <NA> LOC106990386 XR_001433683.1 LOC106990386
        [2]        NA      <NA> LOC106990386 XR_001433683.1 LOC106990386
        [3]        NA      <NA> LOC106990386 XR_001433683.1 LOC106990386
        [4]        NA      <NA> LOC106990386 XR_001433683.1 LOC106990386
        [5]        NA      <NA>         RTP5 XM_015091617.1         RTP5
        ...       ...       ...          ...            ...          ...
  [1364109]        NA      <NA>       PABPC5 XM_004022513.3       PABPC5
  [1364110]        NA      <NA>       PABPC5 XM_004022513.3       PABPC5
  [1364111]        NA      <NA>       PABPC5 XM_004022513.3       PABPC5
  [1364112]        NA         0       PABPC5 XM_004022513.3       PABPC5
  [1364113]        NA         0       PABPC5 XM_004022513.3       PABPC5

So when you make a TxDb, your gene IDs are gene symbols rather than either NCBI Gene IDs or Ensembl Gene IDs, and annotatePeak can't do anything with that. The fix is to replace the gene symbols with actual Gene IDs, which is relatively trivial to do.

> library(RMySQL)
> library(DBI)
> z <- import("oviAri4.ncbiRefSeq.gtf.gz")
> con <- dbConnect(MySQL(), user = "genome", host = "genome-mysql.soe.ucsc.edu", dbname = "oviAri4")
> ids <- paste(unique(mcols(z)$transcript_id), collapse = "','")
> llids <- dbGetQuery(con, paste0("select locusLinkId, ncbiRefSeq.name from ncbiRefSeq inner join ncbiRefSeqLink on ncbiRefSeq.name=ncbiRefSeqLink.id where ncbiRefSeq.name in ('", ids, "');"))
## NOTE HERE that the parentheses have extra single quotes! In other words, it goes parenthesis, single quote, double quote
> mcols(z)$gene_id <- llids[match(mcols(z)$transcript_id, llids$name),1]
## there are a couple missing IDs we have to fix
> badids <- mapIds(orgdb, mcols(z)[is.na(mcols(z)$gene_id),"gene_name"], "GID","SYMBOL")
> mcols(z)$gene_id[is.na(mcols(z)$gene_id)] <- badids
> txdb <- makeTxDbFromGRanges(z)
 Warning message:
 In .get_cds_IDX(mcols0$type, mcols0$phase) :
   The "phase" metadata column contains non-NA values for features of type
   stop_codon. This information was ignored.
> transcriptsBy(txdb)
 GRangesList object of length 24491:
 $`100034665`
 GRanges object with 7 ranges and 2 metadata columns:
       seqnames          ranges strand |     tx_id        tx_name
          <Rle>       <IRanges>  <Rle> | <integer>    <character>
   [1]    chr21 5510049-5691733      + |     41334 XM_015103062.1
   [2]    chr21 5510049-5691733      + |     41335 XM_015103065.1
   [3]    chr21 5510252-5691733      + |     41336 XM_015103061.1
   [4]    chr21 5510252-5691733      + |     41337 XM_015103064.1
   [5]    chr21 5512619-5691733      + |     41338 XM_012101598.1

And now that should work for you.
ADD COMMENT •
0
Entering edit mode

I forgot, you need a library(rtracklayer) before the import step.

ADD REPLY •
0
Entering edit mode

Amazing! Thank you, I would have never noticed gene_id did not contain Gene IDs!!!

This worked perfectly, and I easily converted the csAnno object to a dataframe and exported! Also, this has highlighted how much I still have to learn - your use of the phrase "relatively trivial" before casually hooking into and querying the ucsc database 😂!

Thank you again! Chris

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 604 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6