【R】Mygene.infoで遺伝子情報を取得（改）

##前回からの変更点
R3.5.3以上ではBioconductor経由のインストール方法がBiocManagerパッケージを使う方法に変わった
R 3.4.3でのインストール方法は→[mygeneパッケージで遺伝子の情報を取得する](https://qiita.com/motthy/items/4138e8d686dd2023b36f)

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

##使用したRのバージョン
- R 3.5.3

##パッケージ
- [mygene](https://www.bioconductor.org/packages/release/bioc/html/mygene.html) (Bioconductor) 

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

```{r}
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()
BiocManager::install("mygene")
```

- [available fields](http://docs.mygene.info/en/latest/doc/data.html#available-fields) に取得できる情報の一覧がある

- 例えば`field = "generif"`とすると[NCBI gene GeneRIF](https://www.ncbi.nlm.nih.gov/gene/about-generif)のテキストとpubmedのNo.が取得できる

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

```{r}
x <- read.table("genelist.txt")
x

##      V1
## 1 BIRC6
## 2  RND3
## 3  NSD1
## 4 CNNM1
## 5 MCMBP
```
<b>コード</b>

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

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

- この例ではqueryとしてsymbolを用いているが, [available fields](http://docs.mygene.info/en/latest/doc/data.html#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
.....
```