【R】mygeneパッケージで遺伝子の情報を取得する

概要

R用クライアントmygeneパッケージを用いて, 遺伝子情報取得ウェブサービスMyGene.infoで遺伝子に関する要約文(GeneRIF)を取得する

使用したRのバージョン

• R 3.4.3

パッケージ

• mygene (Bioconductor)

mygeneパッケージのインストール

source("https://bioconductor.org/biocLite.R")
biocLite("mygene")

• available fields に取得できる情報の一覧がある

• 例えばfield = "generif"とするとNCBI gene GeneRIFのテキストとpubmedのNo.が取得できる

Rのコード

遺伝子名 (symbol) でGeneRIFを古いものから5個取得する

x <- read.table("genelist.txt")
x

##      V1
## 1 BIRC6
## 2  RND3
## 3  NSD1
## 4 CNNM1
## 5 MCMBP

コード

• in_fに遺伝子名リスト (genelist.txt) を指定

• 拡張子をのぞいた文字列をprefixとして, prefix_info.tsvというタブ区切りのテキストができる

• この例ではqueryとしてsymbolを用いているが, available fields にあるidであればqueryMany(scopes = "ensembl.transcript")のようにscopesオプションで指定できる

get_generif.r

#パッケージ読み込み
library(mygene)

#遺伝子名リストのファイルを指定
in_f <- "genelist.txt" 

#出力ファイル名を指定
prefix <- gsub("\\..+$","",in_f)
out_f <- paste(prefix,"info.tsv",sep = "_")

data <- read.table(in_f,stringsAsFactors = FALSE)
d_list <- c(data$V1)

#出力ファイルのヘッダー
header <- c("query","gene_id","pubmed","text")

#symbolをgene idに変換
geneid_list <- queryMany(d_list,scopes = "symbol",
                         fields="entrezgene", species="human")

#はじめにヘッダーを書き込む
write(header,file = out_f,ncolumns = 4,sep = "\t")

#GeneRIFのtextを古いものから5個まで抽出して書き込む
for (i in 1:nrow(geneid_list)){
  if (is.na(geneid_list$entrezgene[i])){
    na_text <- paste(geneid_list$query[i],"NA","NA","NA",
                     sep = "\t")
    write(na_text,file = out_f,append = TRUE)
  }
  else{
    gene <- getGene(geneid_list$entrezgene[i],fields = "generif")


    if (length(gene[[1]]$generif) < 5){
      gr_text <- rep(NA,length(gene[[1]]$generif))
      for (j in 1:length(gene[[1]]$generif)){
        gr_text[j] <- paste(geneid_list$query[i],
                            geneid_list$entrezgene[i],
                        gene[[1]]$generif[[j]]$pubmed,
                        gene[[1]]$generif[[j]]$text,
                            sep = "\t")
      }
      write(gr_text,file = out_f,append = TRUE)
    }

    else{
      gr_text <- rep(NA,5)
      for (j in 1:5) {
        gr_text[j] <- paste(geneid_list$query[i],
                            geneid_list$entrezgene[i],
                        gene[[1]]$generif[[j]]$pubmed,
                        gene[[1]]$generif[[j]]$text,
                        sep = "\t")
      }
      write(gr_text,file = out_f,append = TRUE)
    }
  }
}

結果

symbol  gene_id pubmed  text
BIRC6   57448   14765125    Overexpression of Nrdp1 promotes ubiquitination and proteasomal degradation of BRUCE.
BIRC6   57448   15300255    Apollon binds to, ubiquitinates and facilitates proteasomal degradation of SMAC and caspase-9, results suggesting that Apollon has an essential function in preventing SMAC-induced apoptosis
.....