コンテナ環境にVirFinderのインストール

VirFinder というソフトウェアのインストールに若干つまづいたのでメモしておく。

VirFinder概要

論文：https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-017-0283-5

• ウイルス由来の配列を抽出するソフトウェアで、VirSorter(2015, 引用>200)に次いで使われている(2017, 引用>60)
• 後発なので論文でVirSorterとの比較もしていて、特にショートコンティグでVirSorterより精度が出るとのこと
• VirSorterはめちゃくちゃアドホックに複数の手法(既知ウイルスとの相同性、遺伝子の長さ、未知遺伝子の割合など)を組み合わせて独自の基準を算出し、ウイルスかどうかを判断するが、VirFinderはウイルスと非ウイルスのデータを集めて単純な機械学習の判別問題として解いているのがイケてる
• 手法としてはL1 正則化を加えたロジスティック回帰。SVMとかナイーブベイズも試したが結果単純なロジスティック回帰が一番よかったそうだ
• VirSorterはパイプラインの中身を理解しないと理解できない複数のアウトプットがあったりしてややこしく、それらの解釈がむずいのだがVirFinderは単純に確率出してくれるだけなのでわかりやすい
• Jed A. Fuhrman 著者に入ってる、Fuhrman先生は海洋微生物学の巨人なのですがめちゃくちゃ時代に適応してバリバリにバイオインフォがらみの仕事しててすごい

インストール

• Rのパッケージとして提供されている ので、のでコンテナ使うのが無難だろう
• rockerを元にメタゲノム解析用にRのいろんなバッケージを突っ込んで雑に使っているコンテナがあるのでそこに入れる(https://cloud.docker.com/u/ryujhbv/repository/docker/ryujhbv/r-for-metagenome)
• 基本的にはVirFinderのリポジトリをcloneして、 VirFinder_1.1.tar.gz をRからインストールすれば良い
• VirFinderのGitHubにはsource("*/biocLite.R") とか biocLite()とか書いてあるがこれは古いやり方。Bioconductor 3.8からは BiocManager を使う。
• 以下DockerFile (from bioconductor 以下の部分がVirFinderのinstall部分)
• CyVerse使うともっと楽にいける説あるがCyVerseよくわかんないので自力で入れた

FROM rocker/verse:3.6.0
MAINTAINER YOHEI_KUMAGAI

# R

RUN apt-get update \
  && apt-get install -y --no-install-recommends \
    lbzip2 \
    libfftw3-dev \
    libgdal-dev \
    libgeos-dev \
    libgsl0-dev \
    libgl1-mesa-dev \
    libglu1-mesa-dev \
    libhdf4-alt-dev \
    libhdf5-dev \
    libjq-dev \
    liblwgeom-dev \
    libpq-dev \
    libproj-dev \
    libprotobuf-dev \
    libnetcdf-dev \
    libsqlite3-dev \
    libssl-dev \
    libudunits2-dev \
    netcdf-bin \
    postgis \
    protobuf-compiler \
    sqlite3 \
    tk-dev \
    wget \
    unixodbc-dev \
  && install2.r --error \
    RColorBrewer \
    RandomFields \
    RNetCDF \
    classInt \
    deldir \
    gstat \
    hdf5r \
    lidR \
    mapdata \
    maptools \
    mapview \
    ncdf4 \
    proj4 \
    raster \
    rgdal \
    rgeos \
    rlas \
    sf \
    sp \
    spacetime \
    spatstat \
    spatialreg \
    spdep \
    geoR \
    geosphere \
    vegan \
    dichromat \
    ## from bioconductor
    && R -e "BiocManager::install(c('rhdf5','qvalue'))" \ # ここからVirSorterのインストール、ただしrhdf5はVirSorter関係ない
    && R -e "install.packages(\"glmnet\", dependencies=TRUE)" \
    && R -e "install.packages(\"Rcpp\", dependencies=TRUE)"
WORKDIR /root
RUN git clone https://github.com/jessieren/VirFinder.git
RUN R CMD INSTALL /root/VirFinder/linux/VirFinder_1.1.tar.gz

動作確認

• 遺伝研スパコンにSingularityに変換した上記コンテナを突っ込む。

(base) [kumay@at028 VirFinder]$ module load singularity
(base) [kumay@at028 VirFinder]$ singularity exec ~/test_container/r-for-metagenome-latest.img R

R version 3.6.0 (2019-04-26) -- "Planting of a Tree"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(VirFinder)
Loading required package: glmnet
Loading required package: Matrix
Loading required package: foreach
Loaded glmnet 2.0-18

Loading required package: qvalue
> predResult <- VF.pred("Input_contigs.fna")
[1] "Input_contigs.fna [1] line1-3 k141_2 flag=1 multi=2.0000 len=336 len 336 score 0.9184 pvalue 0.01"
[1] "Input_contigs.fna [2] line3-5 k141_3 flag=1 multi=1.0000 len=328 len 328 score 0.7846 pvalue 0.0494"
[1] "Input_contigs.fna [3] line5-7 k141_4 flag=1 multi=7.0000 len=303 len 303 score 0.6339 pvalue 0.1167"
...

動いてるっぽい。

こんな感じのRスクリプトを実行すると、ウイルスと予測される確率が高い順にソートして出力してくれる

library(VirFinder)
setwd("/home/kumay/virus_data/VirFinder")
predResult <- VF.pred("Input_contigs.fna")
predResult$qvalue <- VF.qvalue(predResult$pvalue)
write(predResult[order(predResult$qvalue),],file="predicted_virus_sorted_by_qvalues")