Question: locateVariants missing some mutation site
0
gravatar for noeD
9 weeks ago by
noeD70
noeD70 wrote:

Hello!

I am using Variant Annotation package from Bioconductor to annotate my variant but it missed some variants in the output. For example, in var I have the following mutation site :

CHROM      POS REF ALT     
9751     chr6 41184897   G  A

I used GRanges for the genomic locations and their associated annotations:

gr <-  GRanges(seqnames = var$CHROM, ranges = IRanges(as.numeric(var$POS), as.numeric(var$POS)), strand="*")

The output was:

    GRanges object with 1 range and 0 metadata columns:
      seqnames               ranges strand
         <Rle>            <IRanges>  <Rle>
  [1]     chr6 [41184897, 41184897]      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

After that, I used locateVariants to annotate the sites:

gr_var <- locateVariants(gr, txdb, AllVariants(), ignore.strand=T)

where txdb is txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene (library of R).

But I obtained this output:

GRanges object with 0 ranges and 9 metadata columns:
   seqnames    ranges strand | LOCATION  LOCSTART    LOCEND   QUERYID      TXID         CDSID      GENEID       PRECEDEID        FOLLOWID
      <Rle> <IRanges>  <Rle> | <factor> <integer> <integer> <integer> <integer> <IntegerList> <character> <CharacterList> <CharacterList>
  -------
  seqinfo: no sequences

I searched the mutation site in genome Browser https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg19&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr6%3A41184397-41185396&hgsid=612547501_yo5UeqjaPqdVW1X73qUXN4riW5iA and I found its genes: TREML3 (Entrez Gene ID: 340206)

Therefore, I looked into txdb if it was reported the 340206 gene ID:

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exon_columns <- columns(txdb)
exon_columns <- exon_columns[which(exon_columns %in% c("EXONID", "EXONNAME", "GENEID"))]
all_exons <- exons(txdb, columns=exon_columns)
temp <- as.data.frame(all_exons)
temp[temp$GENEID == "340206", ]

And I found that the Gene_ID is included in txdb:

seqnames    start      end width strand EXONID EXONNAME GENEID
    95166     chr6 41176292 41176892   601      -  95166     <NA> 340206
    95167     chr6 41177337 41177468   132      -  95167     <NA> 340206
    95168     chr6 41183500 41183557    58      -  95168     <NA> 340206
    95169     chr6 41184865 41184915    51      -  95169     <NA> 340206
    95170     chr6 41185391 41185685   295      -  95170     <NA> 340206

Do you have any suggestions as to how to solve?

Thank you in advance

Best

ADD COMMENTlink written 9 weeks ago by noeD70
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 658 users visited in the last hour