【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
.....