NCBI2R 패키지를 이용해 게놈 주석 불러오기

2012-11-29
Annotation Genome NCBI NCBI2R R

관련분석(association study)을 하다 보면 필연적으로 SNP와 게놈정보를 많이 이용하게 됩니다. 보통 클라이언트에게 제출하는 분석결과 보고서에는 각 SNP에 대한 위치정보와 대응하는 유전자 이름을 같이 넣게 되는데 보통 DNA chip 제조회사가 제공하는 주석(annotation) 정보를 이용하게 됩니다.

저는 관련분석 시 주로 plink의 분석 결과 파일을 R로 불러와 그래프 작성 및 보고용 파일을 만드는데, plink의 결과 파일에는 유전자에 대한 정보가 포함되어 있지 않기 때문에 별도로 DNA chip 제조회사가 제공하는 주석 파일에서 필요한 부분을 추출하여 보고용 파일을 만들곤 합니다. 대략 아래와 같이

클라이언트에게 돌려 주는 결과

클라이언트에게 돌려 주는 결과

또한, 클라이언트에 따라 특정 유전자영역에 대해서만 분석을 의뢰하는 예도 있어서 특정 유전자 영역 내에 있는 SNP의 목록이 필요할 때도 있습니다. 이때도 보통 리눅스의 셸 스크립트를 이용하거나 R에서 작업합니다.

그리고 간혹 NCBI의 dbSNP에 접속해 SNP 정보를 확인하는데 관련 지식이 미천한지라 제가 필요로 하는 그 이상의 정보가 존재하기 때문에 필요한 부분만 가지고 올 수 없을까 하는 생각을 했었습니다.

dbSNP 검색결과 화면

dbSNP 검색결과 화면

물론 R에는 온라인 데이터베이스를 탐색할 수 있는 많은 패키지가 존재하며 대표적인 것이 생물정보학을 위한 Bioconductor라는 프로젝트입니다.

그러던 중 NCBI2R 라는 패키지를 이용해 NCBI의 필요한 정보만을 R로 가지고 오는 방법을 알게 되어 소개하고자 합니다.

  • NCBI2R 패키지 인스톨
> devtools::install_github("cran/NCBI2R")
  • NCBI dbDNP에서 SNP 정보 가지고 오기
> library(NCBI2R)
> snplist <- c("rs12345", "rs333", "rs624662")
> mySNPs <- GetSNPInfo(snplist)
#> NCBI2R beta version. Your version: 1.4.7 Latest Release 1.4.6
> mySNPs
#>     marker        genesymbol        locusID chr   chrpos
#> 1  rs12345          CRYBB2P1           1416  22 25459492
#> 2    rs333 LOC102724297,CCR5 102724297,1234   3 46373456
#> 3 rs624662      LOC101927120      101927120  11 56853544
#>                                     fxn_class      species dupl_loc
#> 1                       nc-transcript-variant Homo sapiens         
#> 2 frameshift-variant,intron-variant,reference Homo sapiens         
#> 3                              intron-variant Homo sapiens         
#>   current.rsid flag
#> 1      rs12345    0
#> 2        rs333    0
#> 3     rs624662    0
  • 특정 유전자 안에 있는 SNP 목록 가지고 오기
> locusID <- GetIDs("GNA12[SYM] human")
> mySNPs <- GetSNPsInGenes(locusID)
#> Error in getListFromXML(webget, MaxRet = MaxRet, sme = sme, smt = smt): NCBI2R error: parsing XML problem
> head(mySNPs)
#>     marker        genesymbol        locusID chr   chrpos
#> 1  rs12345          CRYBB2P1           1416  22 25459492
#> 2    rs333 LOC102724297,CCR5 102724297,1234   3 46373456
#> 3 rs624662      LOC101927120      101927120  11 56853544
#>                                     fxn_class      species dupl_loc
#> 1                       nc-transcript-variant Homo sapiens         
#> 2 frameshift-variant,intron-variant,reference Homo sapiens         
#> 3                              intron-variant Homo sapiens         
#>   current.rsid flag
#> 1      rs12345    0
#> 2        rs333    0
#> 3     rs624662    0
  • NCBI에서 유전자 정보 가지고 오기
> GetGeneInfo(locusID)
#>   locusID org_ref_taxname org_ref_commonname   OMIM     synonyms
#> 1    2768    Homo sapiens              human 604394 RMP gep NNX3
#>   genesummary                   genename phenotypes pathways GeneLowPoint
#> 1             G protein subunit alpha 12                          2728105
#>   GeneHighPoint ori chr genesymbol                               build
#> 1       2870656   -   7      GNA12 Homo sapiens Annotation Release 108
#>           cyto approx
#> 1 7p22.3-p22.2      0
  • HapMap 프로젝트 데이터를 이용해 LD 정보 가지고 오기
> LD <- GetLDInfo("4", 123456789, 123556789)
> head(LD)

물론 많은 양의 SNP에 대한 정보를 가지고 올 때는 시간이 오래 걸릴 수 있습니다. 그러나 관련분석 결과 통계적으로 유의한 SNP는 보통 많아야 수십 개기 때문에 NCBI2R 패키지를 이용해 작업하는 편이 훨씬 편리하리라 생각합니다. 다만 DNA chip 제조회사가 제공하는 주석 파일의 위치정보와 NCBI의 위치정보가 다를 수도 있으므로 주의해야 합니다. NCBI2R 패키지의 모든 함수목록과 사용법은 이곳에서 확인하실 수 있습니다.

참고사이트: Genome annotation with NCBI2R

꼬리말: 이 글은 RStudio와 knitr 패키지를 이용해 html 파일을 만들고 워드프레스에서 편집기에서 약간의 수정을 하였습니다. 처음부터 워드프레스 에서 글 작성하는 것에 비하면 완전 딴세상이네요. 강추!!!

comments powered by Disqus