京橋のバイオインフォマティシャンの日常

南国のビーチパラソルの下で、Rプログラムを打ってる日常を求めて、、Daily Life of Bioinformatician in Kyobashi of Osaka

【R言語とWebスクレイピングとアニメーショングラフ】仮想通貨や株価変動をバー・チャート・レースで観察してみた件【その1 日経225・Nikkei225 編】

はじめに

グラフをアニメーションにすることで、静的なグラフでは得られなった、気づきを与えてくれます。

さらには、ある時系列の動的な変化を眺めることで、データのダイナミクスをより正確に感じれます。

同ブログは、「アニメーショングラフ」推しです。

YouTubeを見てて、以前から近いものを作りたいと思ってて、 この記事では、「バー・チャート・レース」を取り上げてみます。

第一弾として、アニメーショングラフにするデータは、日経平均(日経225)の時系列データ(日足)を使います。

さぁ、アニメーショングラフを設定していきましょう。

実行環境

実行環境
macOS Big Sur (バージョン11.5.2)
MacBook Air (M1, 2020)
チップ Apple M1
メモリ 16GB

Rパッケージの準備

関連のRパッケージを準備します。結構あります。 ここをスキップして、以下のスクリプトを読み込むのでもOKです。

#パッケージ・インストール
pack <- c("rvest", "quantmod", "magrittr", "purrr", "tidyr", "ggplot2", "treemapify", "gganimate", "gapminder", "gifski")
install.packages(pack[!(pack %in% unique(rownames(installed.packages())))])

#ロード
library(rvest)
library(quantmod)
library(magrittr)
library(purrr)
library(tidyr)
library(ggplot2)
library(treemapify)
library(gganimate)
library(gapminder)
library(gifski)

以下のgetNikkei225.Rスクリプトを読み込んでも、Rパッケージがインストールできます。

同時に、今回使用する、getNikkei225_list、Nikkei225_ChartData、Nikkei225_ChartData_mod、Nikkei225_animationといった関数群を読み込みます。

#関連関数群の読み込み
source("https://gist.githubusercontent.com/kumeS/c9883239311546c5ab1753676fa86514/raw/4c7a3a6c2520afdf7dc0e56c1fc5a161ce74ee20/getNikkei225.R")

ではでは、早速、実行してみましょう。

日経225銘柄のデータ取得

日経平均株価 - ウィキペディア

まず、getNikkei225_list関数では、ウェブスクレイピングを介して、ウィキペディアの該当ページから日経225の銘柄リストを取得します。時々、ページが更新されて、表の構成が変わるので、厄介です・・・

そして、次に、Nikkei225_ChartData関数でデータを取得します。

#Wikipediaページから日経225銘柄リストを取得する
Nikkei225List <- getNikkei225_list()

#表示
head(Nikkei225List)
#  Ticker                        Company Sector
#1   2002           日清製粉グループ本社   食品
#2   2269           明治ホールディングス   食品
#3   2282                       日本ハム   食品
#4   2501       サッポロホールディングス   食品
#5   2502 アサヒグループホールディングス   食品
#6   2503         キリンホールディングス   食品

#銘柄のセクター数
table(Nikkei225List$Sector)

#日経225のデータ取得: 実行時間 5分くらい
ChartData <- Nikkei225_ChartData(Dat=Nikkei225List, term=c("2022-01-01", "2022-12-31"))

#いったん、RDS保存
saveRDS(ChartData, "Nikkei225_ChartData.Rds")
#ChartData0 <- readRDS("Nikkei225_ChartData.Rds")

#または、Nikkei225_ChartData.Rdsをダウンロードできるようにもしています。
#system("which svn")
#system("svn export  https://github.com/kumeS/Blog/trunk/R_BarChart/Nikkei225_ChartData.Rds")

#表示
head(ChartData0)
#                         Company Ticker Sector       date close dclose dclose2 ranking
#1           日清製粉グループ本社   2002   食品 2022-01-04   100      0 #F0F5F0       1
#2           明治ホールディングス   2269   食品 2022-01-04   100      0 #F0F5F0       2
#3                       日本ハム   2282   食品 2022-01-04   100      0 #F0F5F0       3
#4       サッポロホールディングス   2501   食品 2022-01-04   100      0 #F0F5F0       4
#5 アサヒグループホールディングス   2502   食品 2022-01-04   100      0 #F0F5F0       5
#6         キリンホールディングス   2503   食品 2022-01-04   100      0 #F0F5F0       6

このステップで、日経225のうちで、224銘柄のリストとデータが取得できました。

日経225銘柄でのバー・チャート・レース

日経225 上位20銘柄でのバー・チャート・レース作成

次に、日経225の上昇20銘柄に絞って、バー・チャート・レースを作成します。

Nikkei225_ChartData_mod関数の引数「top = 20」で、日経225の上昇20銘柄のデータフレームにします。 次に、Nikkei225_animation関数を実行して、チャート設定をします。また、日本語設定も別途します。

#データ補正: 上位20銘柄を取得する
ChartData1 <- Nikkei225_ChartData_mod(ChartData0, top = 20)

#表示
head(ChartData1)
#                  Company Ticker             Sector              CompanyTicker
#8353 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#8577 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#8801 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#9025 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#9249 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#9473 DOWAホールディングス   5714 非鉄金属・金属製品 DOWAホールディングス(5714)
#                                CompanySector       date close dclose dclose2 ranking
#8353 DOWAホールディングス(非鉄金属・金属製品) 2022-03-01 111.7   11.7 #F0F5F0      17
#8577 DOWAホールディングス(非鉄金属・金属製品) 2022-03-02 118.9   18.9 #D2E3D2      10
#8801 DOWAホールディングス(非鉄金属・金属製品) 2022-03-03 125.3   25.3 #D2E3D2       8
#9025 DOWAホールディングス(非鉄金属・金属製品) 2022-03-04 122.0   22.0 #D2E3D2       8
#9249 DOWAホールディングス(非鉄金属・金属製品) 2022-03-07 122.0   22.0 #D2E3D2      11
#9473 DOWAホールディングス(非鉄金属・金属製品) 2022-03-08 112.7   12.7 #F0F5F0      12

#日本語設定: 必須
ggplot2::update_geom_defaults("text", list(family = "HiraKakuPro-W3"))
ggplot2::update_geom_defaults("label", list(family = "HiraKakuPro-W3"))

##実行
p1 <- Nikkei225_animation(ChartData1, Label="Nikkei225 Top 20 銘柄")

#gif出力: 実行時間 30秒くらい
animate(p1, fps = 2, duration = 125, width = 600, height = 600,
        renderer = gifski_renderer("Nikkei225_BarChartRace_animation_p1.gif"))

日経225 上位20銘柄でのバー・チャート・レース

日経225 下位20銘柄でのバー・チャート・レース作成

次に、日経225の下位20銘柄に絞って、バー・チャート・レースを作成します。

ここで、Nikkei225_ChartData_mod関数の引数「top = -20」でマイナス設定にすると 日経225の下位20銘柄のデータフレームを得ます。

#データ補正: 下位20銘柄を取得する
ChartData2 <- Nikkei225_ChartData_mod(ChartData0, top = -20)

#日本語設定: 必須
ggplot2::update_geom_defaults("text", list(family = "HiraKakuPro-W3"))
ggplot2::update_geom_defaults("label", list(family = "HiraKakuPro-W3"))

##実行: marの`unit(c(top, right, bottom, left), "cm")`のc(top, right, bottom, left)に対応
p2 <- Nikkei225_animation(ChartData2, Label="Nikkei225 Lowest 20 銘柄",
                          mar=c(0.5, 1.5, 0.5, 0.75), Size1=4, Size2=5)

#gif出力: 実行時間 30秒くらい
animate(p2, fps = 2, duration = 125, width = 600, height = 600,
        renderer = gifski_renderer("Nikkei225_BarChartRace_animation_p2.gif"))

日経225 下位20銘柄でのバー・チャート・レース

まとめ

今回、日経225銘柄の2022年11月まで値動きのバー・チャート・レースを作成しました。 よく眺めてみると、日本株も結構上昇している銘柄が多い印象ですね。。

実際のチャート設定としては、レースへの参加者数、つまりはグループ数が多いと、レースのアニメーションはシンプルにした方がいいかも知れないと思うところです。

また、設定を間違い、レンダリングにかなり時間が取られ、今回は日経225のみでタイムオーバーでした。

参考資料

blog.hoxo-m.com

ill-identified.hatenablog.com

stackoverflow.com

www.anarchive-beta.com

towardsdatascience.com

www.r-bloggers.com

ggplot2.tidyverse.org

flourish.studio

stackoverflow.com

Rコード: getNikkei225.R

gist.github.com

シーケンスリードのクオリティコントロール解析: fastqc、MultiQC、quack、fastpの基本設定と使い方

はじめに

この記事では、シーケンスリードのクオリティコントロール解析、つまりは品質管理を行うツール群を紹介します。品質管理は、シーケンスデータを解析する前の重要なステップとなります。

データ解析は、クオリティコントロールに初まり、クオリティコントロールに終わる。まさに、柔の精神ですね。

京橋の仙人さん、クオリティコントロールって何ですか?

うむ、まずはデータを知り、柔らかに対応するということじゃ。

前置きはこれくらいで、この記事では、fastqc、MultiQC、Quack、fastpについての基本設定とそれらの使い方を概説します。

実行環境

macOS Big Sur (バージョン11.5.2)
MacBook Air (M1, 2020)
チップ Apple M1
メモリ 16GB

テストデータのダウンロード

GAGE (Genome Assembly Gold-standard Evaluation)から、 Staphylococcus aureus (黄色ブドウ球菌)のシーケンスデータをダウンロードします。

簡単なサンプル情報は下記に示します。

# テストデータ
# Staphylococcus aureus: Data download page
# Library 1: Fragment
# Avg Read length: 101bp
# Insert length: 180bp
# # of reads: 1,294,104

まずは、Read1とRead2のFastqファイル(frag_1.fastq.gz と frag_2.fastq.gz)をそれぞれwgetコマンドでダウンロードします。

#ダウンロード
wget http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/frag_1.fastq.gz http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/frag_2.fastq.gz

#ファイルの詳細表示
ls -lh | grep "fastq.gz"
#-rw-r--r--  1 sas  staff    61M  4 15  2011 frag_1.fastq.gz
#-rw-r--r--  1 sas  staff    65M  4 15  2011 frag_2.fastq.gz

ls -lh | grepでデータサイズを確認したところ、fastq.gz形式で約60Mバイトでした。

それでは、このデータを使って、各ツールでのQCを実施していきます。

関連記事

以下の過去の関連記事です。

skume.net

Javaの設定: スキップ可です

Mac環境にJavaをインストールします。

Javaのインストールは過去記事を参照ください。

skume.net

ここでは、パスとバージョンの確認をします。

#パス確認
which java
#/usr/bin/java

#バージョン確認
java --version
#openjdk 18.0.2 2022-07-19
#OpenJDK Runtime Environment Homebrew (build 18.0.2+0)
#OpenJDK 64-Bit Server VM Homebrew (build 18.0.2+0, mixed mode, sharing)

fastqcの使い方

Fastqcのインストール

Mac Homebrewでインストールするか、あるいは、 本家のBabraham bioinformatics からバイナリファイルをダウンロードすることができます。

さっき試したところ、本家からのダウンロードは通信が遅過ぎました。 無難に、Homebrewでインストールするのがお勧めです。

Homebrewでのインストールの場合には、 openjdk--19.0.1.arm64_big_sur.bottle.tar.gz が一緒にインストールされました。

バージョンは、どちらも同じ、FastQC v0.11.9です。

##FastQCインストール
#(1) Homebrewでインストール
brew install fastqc

#(2) Fastqcバイナリファイルを本家からダウンロード
wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.zip

#(3) Mac用のdmgファイルを本家からダウンロード
wget https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.dmg

今回は、私自身は(1)のHomebrewでインストールしました。

Fastqcの基本動作

Fastqcをインストールした後、zsh環境で実行しています。

#パス確認
which fastqc
#/opt/homebrew/bin/fastqc

#バージョン確認
fastqc -v
#FastQC v0.11.9

#ヘルプ表示
fastqc -h
#出力省略

Fastqcの実行

それでは、いよいよ、テストデータ(fastq.gzファイル)に対して、fastqcを実行してみます。

実行方法としては、個別にファイルを指定する方法と、ワイルドカードでファイルを与えてる方法があります。

#fastq.gzを個別に指定する
fastqc --nogroup -o "./" frag_1.fastq.gz frag_2.fastq.gz

#ワイルドカード表記も使える
fastqc --nogroup -o "./" *.fastq.gz

#htmlファイルの表示
ls | grep "fastqc.html"
#frag_1_fastqc.html
#frag_2_fastqc.html

--nogroupはクオリティスコアの丸め込みを防ぐオプションで、-oの引数では出力フォルダを指定します。 "./"はカレントディレクトリ出力なので省略可能です。

次に、ターミナル上からのhtmlファイルの開き方を説明します。

ブラウザは、Google Chromeか、Firefoxかで開くのがお勧めです(経験上、Safariは非推奨)。

#デフォルト・ブラウザで開く
open frag_1_fastqc.html frag_2_fastqc.html

#Google Chrome指定で開く
open frag_1_fastqc.html frag_2_fastqc.html -a Google\ Chrome

#Firefox.app指定で開く
open frag_1_fastqc.html frag_2_fastqc.html -a Firefox.app

FastQC Report

以下のURLで結果レポートをみえます。

frag_1_fastqc.html: https://kumes.github.io/Blog/sh_fastqc/frag_1_fastqc.html

frag_2_fastqc.html: https://kumes.github.io/Blog/sh_fastqc/frag_2_fastqc.html

fastqcレポートの見方

bi.biopapyrus.jpさんのページが分かり易いので、そちらをご参照のこと。

bi.biopapyrus.jp

fastqcの原著

Andrews, S. (2010). 
FastQC:  A Quality Control Tool for High Throughput Sequence Data [Online]. 

www.bioinformatics.babraham.ac.uk

MultiQCの使い方

次には、複数の分析結果を1つに集約して分析・レポート結果を複合化できる、MultiQCについて紹介します。

例えば、Fastqcコマンド実行後に、生成されたレポートhtmlのあるディレクトリ内でMultiQCを実行します。

MultiQCのインストール

MultiQCは、PyPI(The Python Package Index)のリポジトリから取得するので、 ターミナルを起動して、pipコマンドでインストールします。

#pip3のパス確認
which pip3
#/opt/homebrew/bin/pip3

#PyPIからインストールする
pip3 install multiqc

#multiqcパス確認
which multiqc
#/opt/homebrew/bin/multiqc

#multiqcバージョン確認
multiqc --version
#multiqc, version 1.13 (759ab2d)

#multiqcヘルプ表示
multiqc -h
#省略

MultiQCの実行

それでは、MultiQCを実行します。htmlレポートがあるディレクトリに移動して、MultiQCを実行することで、MultiQCのレポートが生成されます。

#fastqcファイルの確認
ls | grep "fastqc.html"                      
#frag_1_fastqc.html
#frag_2_fastqc.html

#MultiQC実行
multiqc .

#  /// MultiQC 🔍 | v1.13 (759ab2d)
#
#|           multiqc | Search path : /Users/sas/Desktop/fastqc
#|         searching | ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 15/15  
#|            fastqc | Found 2 reports
#|           multiqc | Compressing plot data
#|           multiqc | Report      : multiqc_report.html
#|           multiqc | Data        : multiqc_data
#|           multiqc | MultiQC complete

実行が完了すると、multiqc_report.htmlmultiqc_dataフォルダが出力されます。

それでは、レポート結果multiqc_report.htmlを開いていみます。

open multiqc_report.html

下図のレポート結果multiqc_report.htmlが表示されます。

https://kumes.github.io/Blog/sh_fastqc/multiqc_report.html

MultiQCの原著

MultiQC: Summarize analysis results for multiple tools and samples in a single report.
Philip Ewels, Måns Magnusson, Sverker Lundin and Max Käller
Bioinformatics (2016)
doi: 10.1093/bioinformatics/btw354
PMID: 27312411

MultiQC: summarize analysis results for multiple tools and samples in a single report | Bioinformatics | Oxford Academic

quackの使い方

quackは、高速のFASTQ品質評価ツールです。 出力結果は最小限に、実行速度にこだわった開発理念になっています。

quackのインストール

quackは、GitHubからダウンロードします。 binフォルダにquackのバイナリがあり、そのバイナリが使うか、、 もし使えなかったら、ソースからコンパイルしてと記載されていました。

gitコマンドでリポジトリをダウンロードして、quackのバイナリが実行できるかを確認してみます。

quackのバイナリは、Mac版である./quack/bin/MacOS_10.14.6/quackを指定します。

#gitコマンドのパス確認
which git
#/usr/bin/git

#リポジトリのダウンロード
git clone https://github.com/IGBB/quack.git

#実行権限の付与
chmod +x ./quack/bin/MacOS_10.14.6/quack

#/opt/homebrew/binに移動しておく(パス設定済み)
mv ./quack/bin/MacOS_10.14.6/quack /opt/homebrew/bin

#ヘルプ表示
quack --help
#Usage: quack [OPTION...]
#quack -- A FASTQ quality assessment tool
#
#  -1, --forward file.1.fq.gz      Forward strand
#  -2, --reverse file.2.fq.gz      Reverse strand
#  -a, --adapters adapters.fa.gz   (Optional) Adapters file
#  -n, --name NAME                 (Optional) Display in output
#  -u, --unpaired unpaired.fq.gz   Data (only use with -u)
#  -?, --help                      Give this help list
#      --usage                     (use alone)
#  -V, --version                   Print program version (use alone)
#Report bugs to <thrash@igbb.msstate.edu>.
#Usage: quack [OPTION...]
#Try `quack --help' or `quack --usage' for more information.

#もしソースからコンパイルする場合
#cd quack/
#make && make test

quackの実行

Quackでは、入力データとして、gzipされたFASTQ形式のファイル(fastq.gz)を渡して、SVG形式の画像出力を標準出力としています。

片側リードか、ペアエンドリードかで、引数が少し変わってきます。

出力結果のSVGファイルはブラウザ(Google Chrome)で開きます。

#片側リードでの実行(name & adapters の設定なし)
quack \
    -u frag_1.fastq.gz \
    > test_u.svg

#ペアエンドリードでの実行(name & adapters の設定なし)
quack \
    -1 frag_1.fastq.gz \
    -2 frag_2.fastq.gz \
    > test_p.svg

#Google Chrome指定で開く
open test_u.svg test_p.svg -a Google\ Chrome

quackの実行結果

出力結果のSVGファイルには、 片側リードの場合、1つの結果パネルが、 ペアエンドリードの場合、2つの結果パネルが表示されます。

quackのレポートの見方

quackのレポート結果の見方(from オフィシャル)

  1. アレイの各カラムに含まれる各塩基の割合を示す塩基量分布。

  2. 各カラムの配列品質分布を示すヒートマップと、アレイ全体の品質スコアの平均値を示すライン

  3. あるスコアに一致する塩基の割合を示すスコア分布グラフで、グラフの左側が100%、右側が0%。スコアの高いデータがグラフの上部に表示されます。

  4. 長さ分布グラフ:ある長さのリードが占める割合を示すグラフ

  5. Adapter content distributionグラフ:アレイ中のAdapter contentの分布を示すグラフ

quackの原著

A. Thrash, M. Arick, and D. G. Peterson, 
“Quack: A quality assurance tool for high throughput sequence data,” 
Analytical Biochemistry, vol. 548, pp. 38–43, 2018. 

doi.org

fastpの使い方

fastpは、FastQファイルの前処理をオールインワンで高速に行うためのツール群を提供しています。また、C++で開発され、マルチスレッドに対応した高性能ツールでもあります。

fastpでは、フィルタリングやトリミングなども処理実施できるパイプラインも提供されていますが、 この記事では、クオリティコントロールの箇所のみを紹介します。

minicondaの諸設定

まずは、condaを設定するために、minicondaをダウンロード・インストールします。

#Python3.9版 minicondaのダウンロード
wget https://repo.anaconda.com/miniconda/Miniconda3-py39_4.12.0-MacOSX-arm64.sh

#minicondaのインストール
bash ./Miniconda3-py39_4.12.0-MacOSX-arm64.sh -b -p $HOME/miniconda

#パス設定を追記
echo 'export PATH="/Users/{Your Name}/miniconda/bin:$PATH"' >> .zshrc 
cat .zshrc

#condaパス確認: HOMEに作成
which conda
#/Users/{Your Name}/miniconda/bin/conda

#zshの場合、実行??
#source $HOME/miniconda/bin/activate
#conda init zsh

#設定: 必須ですね
conda config --add channels bioconda
conda config --env --set subdir osx-64

ややハマったので、このあたりの掲示板を参照しました。arm/M1 Mac版としては配布してないようですね。

stackoverflow.com

stackoverflow.com

fastpのインストール

Minicondaとbiocondaの設定が終わったら、ようやくcondaコマンドで、fastpをインストールします。

#fastpのインストール
conda install -c bioconda fastp

#パス確認
which fastp
#/Users/[Your Home]/miniconda/bin/fastp

#ヘルプ表示
fastp --help
#省略

fastpの実行

fastpには色々なパイプラインがありますが、クオリティコントロールの方法を紹介します。

トリミングなどを実行しない場合には、入力ファイルとしてfastqあるいはfastq.gzのみを指定する。 フィルタリングはされてQC結果が出力されます。

#fastqファイルの確認
ls | grep "fastq.gz"
#frag_1.fastq.gz
#frag_2.fastq.gz

#片側リードの場合
fastp -i frag_1.fastq.gz -j out1_fastp.json -h out1_fastp.html

#ペアエンドリードの場合
fastp -i frag_1.fastq.gz -I frag_2.fastq.gz -j out2_fastp.json -h out2_fastp.html

-jの引数でjsonのファイル名、-hの引数でhtmlのファイル名を指定します。

実行が終わると、html形式のレポートが出力されます。

出力ファイル(-o, -O)をしていないと、処理結果は出力されません。 今回はQCだけなので、指定していません。

下図は、片側リードの結果です。出力結果htmlは各URLから確認できます。

https://kumes.github.io/Blog/sh_fastqc/out1_fastp.html

下図は、ペアエンドリードの結果です。

: htmlhttps://kumes.github.io/Blog/sh_fastqc/out2_fastp.html

出力結果では、クオリティフィルタ後の結果であったり、インサートサイズの分布を計算されたりと、いい感じの作りになってますね。また、グラフィックが、plotlyで作られているのもオツですね。

やや、htmlレポートが長くて、比較が難しいのが厄介かもですが。

ペアエンドリードでのクオリテイィトリミングの実行例

fastpでは、クオリティフィルタリングの前後でQCレポートを生成します(そのため、htmlレポートがやや長い)。

下記の条件は、 塩基レベルでの品質値30以上(-q 30)、3'末端のクオリテイィフィルタリング(-3)、ReadのN塩基数10未満(-n 10)、read1の末尾1塩基の削除(-t 1 )、read2の末尾1塩基の削除(-T 1)、長さフィルタリング25塩基以上( -l 25 )、3スレッド(-w 3)での実行例です。

#ペアエンドリードの場合
fastp -i frag_1.fastq.gz -I frag_2.fastq.gz -3 -q 30 -n 10 -t 1 -T 1 -l 25 -w 3

fastpのオプションの説明

オプション 概要
-i, --in1 read1の入力ファイル名(文字列[=])
-o, --out1 read1の出力ファイル名(文字列[=])
-I, --in2 read2の入力ファイル名(文字列[=])
-O, --out2 read2の出力ファイル名(文字列[=])
-A, --disable_adapter_trimming アダプタのトリミングはデフォルトで有効。このオプションが指定された場合、アダプタのトリミングは無効となる。
-a, --adapter_sequence read1用のアダプタ。SE data: 指定しない場合、アダプタは自動検出される。PE data: R1/R2が重複していないことが判明した場合に使用される。(文字列 [=auto])
--adapter_fasta FASTAファイルを指定し、read1, read2 (PE)をこのFASTAファイル中の全配列でトリミングする (文字列 [=])
-f, --trim_front1 read1の先頭何塩基をトリミングするか(デフォルト: 0) (int [=0])
-t, --trim_tail1 read1の末尾何塩基をトリミングするか(デフォルト: 0) (int [=0])
-F, --trim_front2 read2の先頭何塩基をトリミングするか(デフォルト: 0) (int [=0])
-T, --trim_tail2 read2の末尾何塩基をトリミングするか(デフォルト: 0) (int [=0])
-5, --cut_front 5'末端(先頭)から後方へウィンドウを移動し、その平均品質 < threshold ならウィンドウ内の塩基を削除し、そうでないなら停止する。
-3, --cut_tail 3'末端(末尾)から前方へウィンドウを移動し、その平均品質 < threshold ならウィンドウ内の塩基を削除し、そうでないなら停止する。
-Q, --disable_quality_filtering クオリティフィルタリングはデフォルトで有効。このオプションが指定された場合、品質フィルタリングは無効となる。
-q, --qualified_quality_phred 塩基レベルでの品質値を指定する。デフォルト15で、phred quality >=Q15 でトリミングすることを意味する。(int [=15])
-n, --n_base_limit あるReadのN塩基数がn_base_limitを超えた場合、そのread/pairは破棄される。デフォルトは5 (int [=5])
-e, --average_qual あるReadが、平均品質スコア < avg_qual ならば、そのread/pairは破棄される。デフォルト=0は要求しないことを意味する (int [=0])
-L, --disable_length_filtering 長さフィルタリングはデフォルトで有効。このオプションが指定された場合、長さフィルタリングは無効となる。
-l, --length_required length_requiredより短い、Readは捨てられる。 デフォルトは15。
-w, --thread ワーカースレッド数、デフォルトは3 (int [=3])

fastpの原著

Shifu Chen, Yanqing Zhou, Yaru Chen, Jia Gu; 
fastp: an ultra-fast all-in-one FASTQ preprocessor
Bioinformatics, Volume 34, Issue 17, 1 September 2018, Pages i884–i890, 

doi.org

まとめ

4つのクオリティチェックのツールをインストールから基本的な使い方を概説しました。

最近は、htmlレポートが標準出力なので、見易いですね。 下流解析から戻って再確認するにも、human-Readableって良いですね。

参考資料

formulae.brew.sh

kazumaxneo.hatenablog.com

github.com

multiqc.info

【R言語と時系列解析】Prophetパッケージを利用した、時系列データの未来予測: ブログ利用者数の推移予測とか、ビットコインや日経225の価格トレンド推移予測とか

はじめに

Prophetパッケージは、Meta社(当時、Facebook社)が2017年頃に開発した時系列予測のオープンソースのライブラリです。 R実装版もあり、今回そのパッケージを使用しました。

実際、私も実務で経験しましたが、、 時系列データの予測アルゴリズムの構築は途轍もなく困難を伴います。 その行為はまさしく、途方もなく、雲を掴む感じですね。

Prophetパッケージは、そのような時系列モデル構築の難易度を下げて、ユーザーをアシストしてくれます。 また、多くの場合、モデルの細かいチューニングをしなくても、デフォルト設定でも案外良好なアウトプットが得られるのも魅力的なところです。

この記事では基本的な関数紹介が中心ではなくて、より実践的な解析内容を中心に紹介します。 実データとしては、このブログ「京橋のバイオインフォマティシャンの日常」の利用ユーザー数のデータをもとにして、将来のユーザー数の推移予測をやってみました。2022年11月現在、ブログのユーザー数は月で7000-8000人あたりにいます。

また、それに加えて、日経225やビットコインの推移予測の実施事例も紹介したいと思います。

実行環境

macOS Big Sur (バージョン11.5.2)
MacBook Air (M1, 2020)
チップ Apple M1
メモリ 16GB

Prophetのデータフォーマット

まずは、Prophetのデータ形式についてです。

Prophetのインプット・データには、ds列とy列という、2列のデータを持つ、データフレームを準備します。

ds列 (=datestamp列) は時系列の列で、日付ならYYYY-MM-DD、タイムスタンプなら YYYY-MM-DD HH:MM:SS というスタイルになります。y列は時系列のそれぞれに対応する数値データで、予測したい測定値の過去データとなります。

データフォーマットも、非常にシンプルですね。

関連パッケージの準備

それでは、prophetパッケージをインストールします。

開発者版のprophet(ver 1.1.1)をインストールする方がバグが改善されています。

そのため、CRAN版パッケージ・インストールはスキップしてください。

#環境変数の全削除
rm(list=ls())

#関連パッケージ・インストール
pack <- c("devtools", "ggplot2", "dygraphs")
install.packages(pack[!(pack %in% unique(rownames(installed.packages())))])

#CRAN版パッケージ・インストール
#ver 1.0; 2021-03-30
#pack <- c("prophet")
#install.packages(pack[!(pack %in% unique(rownames(installed.packages())))])

#開発者版パッケージ・インストール
#ver 1.1.1 インストール from GitHub (2022.09.08)
devtools::install_github('facebookincubator/prophet', subdir='R')

#ライブラリ準備
library(prophet)
library(ggplot2)
library(dygraphs)

#バージョン確認
packageVersion("prophet")
#[1] ‘1.1.1’

#Calc_residuals関数の読み込み
source("https://gist.githubusercontent.com/kumeS/486d049718bb981f980db651d550b59b/raw/3480fdfcea0d168952b62dc682558a537a24cf8f/Calc_residuals.R")

また、入力値と構築モデルの残差を計算する、Calc_residuals関数も読み込んでください。

ブログ利用ユーザー数の推移を未来予測する

実行には、prophetパッケージのprophet関数を使います。

prophet関数の最初の引数dfには、過去データのデータフレームを指定します。 このデータフレームは、過去データの日付を含むds列、及び過去データの数値データを含むy列を持つ必要があります。

それでは、実際のこのブログの利用ユーザー数のデータを使って、予測を実行してみます。 この時系列データは「2020/4/1」か「ら2022/11/12」まであります。

また、過去1日間の訪問ユーザー数(y_1day列)、過去7日間の訪問ユーザー数(y_7day列)、過去14日間の訪問ユーザー数(y_14day列)、過去28日間の訪問ユーザー数(y_28day列)の列があります。

過去14日間の訪問ユーザー数(y_14day列)を使って、ブログ利用ユーザー数の推移を予測する

まずは、過去14日間の訪問ユーザー数(y_14day列)のデータに対して、prophetを実行しました。

ここで、増加し続ける周期性を示すデータの予測を行う場合には、 「seasonality.mode = "multiplicative"」を指定します。

また、周期性の誤差を考慮するには、mcmc.samplesというパラメータを設定することで、ベイズ推論を行えます。結構計算時間がかかるので、今回は、デフォルトのMAP推定にしています。

それでは、実行していきます。

#利用者データのURL
url <- "https://gist.githubusercontent.com/kumeS/c677d0b3ceb56f74cd109aae6acfd1a8/raw/6cf600582b8ba5bc3812d2a600de959c01548eca/Blog_Users_20200401-20221112.csv"

#データの読み込み
df <- read.csv(url, header = T)

#表示
head(df)
#        ds y_1day y_7day y_14day y_28day
#1 2020/4/1      0      0       0       0
#2 2020/4/2      0      0       0       0
#3 2020/4/3      0      0       0       0
#4 2020/4/4      0      0       0       0
#5 2020/4/5      0      0       0       0
#6 2020/4/6      0      0       0       0

#14dayデータを使用する
df_14day <- data.frame(ds=as.Date(df$ds),
                       y=as.numeric(sub(",", "", df$y_14day)),
                       row.names=1:nrow(df))

#表示
head(df_14day)
#          ds y
#1 2020-04-01 0
#2 2020-04-02 0
#3 2020-04-03 0
#4 2020-04-04 0
#5 2020-04-05 0
#6 2020-04-06 0

#Prophet実行
m <- prophet(df=df_14day,
             growth = "linear",
             n.changepoints = 100,
             changepoint.range = 1,
             changepoint.prior.scale=0.05,
             seasonality.mode = "multiplicative",
             interval.width = 0.95)

#クラス表示
class(m)
#[1] "prophet" "list"   

prophet関数を実行して、mというモデル変数を出力します。

次に、make_future_dataframe関数を使って、予測する期間のデータフレームを生成します。 freqの引数としては、dayweekmonthquarteryearなどを選択できます。

さらには、モデルと予測期間のデータフレームをpredict関数に渡して、 plot関数などを使って予測結果を可視化します。

#365日分のデータフレーム作成
future <- make_future_dataframe(m, periods = 365, freq = "day")

#モデル予測
forecast <- predict(m, future)

#末尾表示
tail(forecast[c('ds', 'yhat', 'yhat_lower', 'yhat_upper')])
#             ds     yhat yhat_lower yhat_upper
#1316 2023-11-07 6524.384   5941.722   7090.318
#1317 2023-11-08 6533.519   5921.415   7090.767
#1318 2023-11-09 6536.082   5960.249   7101.341
#1319 2023-11-10 6554.383   5971.351   7139.919
#1320 2023-11-11 6564.437   5974.211   7160.411
#1321 2023-11-12 6576.142   5978.119   7156.932

#予測結果のプロット
plot(m, forecast)

#文字化け対策
Sys.setlocale("LC_TIME", "en_US")

#トレンド結果
prophet_plot_components(m, forecast)

#実測と予測された結果: 残差の計算
Calc_residuals(df_14day$y,
               forecast$yhat[1:nrow(df_14day)])
#       mse     mae    rmse    MAPE     R2
#1 1666.788 28.3356 40.8263 12.6803 0.9987

予測結果のプロット

青線の予測結果から、2週間ベースで、2022年末に4500人で、2023年明けに5000人を超えることが期待されますね。 また、2023年6月ごろには6000人を目指せそう結果になってますね。

prophet_plot_componentsのプロット結果

トレンドは、2021年の後半から上昇トレンドにありますね。 水曜日と金曜日に多い傾向にありますね。

他のプロット方法について

add_changepoints_to_plot関数で、潜在的変化点を重ね合わせできます。 また、plot_forecast_componentでは、特定の成分をプロットします。

dyplot.prophetでは、dygraphsパッケージによるインタラクティブなプロットは可能です。

#重要な変化点の重ね合わせプロット
plot(m, forecast) + add_changepoints_to_plot(m)

##特定の成分をプロットする
#yhat
plot_forecast_component(m, forecast, name="yhat")
#trend
plot_forecast_component(m, forecast, name="trend")

#dyplotによるインタラクティブ・プロット
dyplot.prophet(m, forecast)

重要な変化点の重ね合わせプロット

特定の成分 yhat をプロットする

特定の成分 trend をプロットする

https://kumes.github.io/Blog/R_prophet/Fig_03_4.html

dyplot.prophet関数は、インタラクティブにプロットできて、使い勝手も良いですね。

過去28日間の訪問ユーザー数(y_28day列)を使って、ブログ利用ユーザー数の推移を未来予測する

次に、過去28日間の訪問ユーザー数(y_28day列)を使って、prophetを実行しました。

ここでは、少し実行内容を省略しています。

#利用者データのURL
url <- "https://gist.githubusercontent.com/kumeS/c677d0b3ceb56f74cd109aae6acfd1a8/raw/6cf600582b8ba5bc3812d2a600de959c01548eca/Blog_Users_20200401-20221112.csv"

#データの読み込み
df <- read.csv(url, header = T)

#28dayデータを使用する
df_28day <- data.frame(ds=as.Date(df$ds),
                       y=as.numeric(sub(",", "", df$y_28day)),
                       row.names = 1:nrow(df))

#Prophet実行
m <- prophet(df=df_28day,
             growth = "linear",
             n.changepoints = 100,
             changepoint.range = 1,
             changepoint.prior.scale=0.05,
             interval.width = 0.95,
             seasonality.mode = "multiplicative")

#365日分のデータフレーム作成
future <- make_future_dataframe(m, periods = 365, freq = "day")

#モデル予測
forecast <- predict(m, future)

#結果のプロット
plot(m, forecast) + add_changepoints_to_plot(m)

#dyplotによるインタラクティブ・プロット
dyplot.prophet(m, forecast)

#実測と予測された結果: 残差の計算
Calc_residuals(df_28day$y,
               forecast$yhat[1:nrow(df_28day)])
#       mse     mae    rmse   MAPE     R2
#1 998.6104 20.8671 31.6008 6.7593 0.9998

重要な変化点の重ね合わせプロット

28日の利用者数の予測結果は、14日の利用者数の予測結果の2倍の値を示します。

来年1月初めに、利用者数が激減しているけども、何か対策はあるのかな。。

https://kumes.github.io/Blog/R_prophet/Fig_04_2.html

何方にせよ、2023年もアップ・トレンドを期待したいです。

日経225インデックスの将来推移を予測してみたところ

https://www.macrotrends.net/2593/nikkei-225-index-historical-chart-data

最近、www.macrotrends.netで、日経225インデックスの67年間のチャートデータを見つけました。

ちょうどいいデータなので、今回使ってみることにしました。

まずは、Gistに成形した日経225インデックスのチャートデータをアップしていますので、それをダウンロードします。

#ダウンロードデータURLの定義
url <- "https://gist.githubusercontent.com/kumeS/0b901711f87e5a6854a6b354559b11ff/raw/d018a8700d150d43a528a347ee0ffb73cee322c8/nikkei-225-index-67-year-historical-chart-data.csv"

#日経225データの読み込み
df <- read.table(url, sep=",", header = TRUE)

#dyplotによるインタラクティブ・プロット
df0 <- df
rownames(df0) <- df0$ds
dygraph(df0) %>% dyRangeSelector()

https://kumes.github.io/Blog/R_prophet/Fig_06.html

2012年以降のデータで、prophetで株価トレンドを予測する

日経225の2012年以降のデータを使用って、prophetで株価トレンドを予測してみました。

changepoint.prior.scaleは、潜在的変化点の自動選択の柔軟性を調整するパラメータです。 大きな値を設定すると、多くの変化点が選択され、小さな値を設定すると、変化点は選択され難くなります。今回、そのスケール値を0.85に変えています。

また、make_future_dataframe関数では、24週間分のデータフレームを作成しています。

#2012年以降のデータを使用する
df2012 <- df[df$ds > "2012-01-01",]

#表示
head(df2012)
#              ds       y
#15605 2012-01-04 8560.11
#15606 2012-01-05 8488.71
#15607 2012-01-06 8390.35
#15608 2012-01-10 8422.26
#15609 2012-01-11 8447.88
#15610 2012-01-12 8385.59

#Prophet実行
m <- prophet(df=df2012,
             growth = "linear",
             n.changepoints = 50,
             changepoint.range = 1,
             changepoint.prior.scale=0.85,
             interval.width = 0.8,
             seasonality.mode = "multiplicative")

#24週分のデータフレーム作成
future <- make_future_dataframe(m, periods = 24, freq="week")

#モデル予測
forecast <- predict(m, future)

#dyplotによるインタラクティブ・プロット
dyplot.prophet(m, forecast)

https://kumes.github.io/Blog/R_prophet/Fig_07.html

まぁ、フィットが悪いんだけど、年明けから日経225が下がったら面白いね。 ふと思うと、アップトレンドをではないので、seasonality.mode = "additive"のほうがフィットは良いのかも・・・

ビットコイン(BTC)の価格推移を予測してみたところ

ここでは、getBTC関数を使って、 2010年7月18日から今日(現在、2022年11月20日)までのBTCの日足データを取得します。

それで、ビットコインの取得データは、2017年以降のデータを使用します。

#BTCのデータ取得
source("https://gist.githubusercontent.com/kumeS/c466e97d42f6facbb8ba5a22ed361bf4/raw/07728e5e4810e53aa8bdcc411ac51603943761df/getBTC.R")

#BTCのデータ取得
BTC <- getBTC()
BTC0 <- BTC[,c("Date", "Close")]
colnames(BTC0) <- c("ds", "y")

#2017年以降のデータを使用する
BTCdata <- BTC0[BTC0$ds > "2020-01-01",]

やっぱり、仮想通貨のように値動きの激しいデータは、パラメータ・チューニングに苦戦しますね。それでもかなりマシですけども。。。

アップトレンドで無いデータということで、seasonality.modeの引数はadditiveにしています。

#Prophet実行
m <- prophet(df=BTCdata,
             growth = "linear",
             n.changepoints = 250,
             changepoint.range = 1,
             changepoint.prior.scale=0.5,
             interval.width = 0.8,
             seasonality.mode = "additive")

#30日分のデータフレーム作成
future <- make_future_dataframe(m, periods = 30, freq="day")

#モデル予測
forecast <- predict(m, future)

#dyplotによるインタラクティブ・プロット
dyplot.prophet(m, forecast)

https://kumes.github.io/Blog/R_prophet/Fig_08.html

色々とパラメータ調整したからかも知れませんが、、 この結果は、興味深いですね。

ビットコインが10000ドルあたりで反転するのか、どうかですね。

まとめ

現在、ブログのユーザー数は月で7000-8000人を推移してますが、今後どう推移してくるのか?、Prophetでの予想通りに増えてくるのか?、1万人の大台を超えてくるのか、楽しみですね。今後、実際のデータとの比較もできるとさらに面白そうだ。。。

あとは、弱気相場も大詰めで、ビットコインの底値がどのあたりになるかも気になりますね。時系列予測がどこまで使えるか、気になりますよね。

補足

Calc_residuals関数

このCalc_residuals関数では、 MSE(平均二乗誤差)、MAE(平均絶対誤差)、RMSE(二乗平均平方根誤差)、MAPE(平均絶対パーセント誤差)、およびR2(決定係数)を計算します。

Calc_residuals <- function(original, predicted, digits=4){
#diff
d0 <- original - predicted
#MSE
d1 <- mean((d0)^2)
#MAE
d2 <- mean(abs(d0))
#RMSE
d3 <- sqrt(mean((d0)^2))
#MAPE
d4a <- abs((d0/original))
d4b <- mean(d4a[!is.infinite(d4a)])*100
#R2
d5 <- round(1-(sum((d0)^2)/sum((original-mean(original))^2)), digits)

return(data.frame(mse=round(d1, digits=digits),
                  mae=round(d2, digits=digits),
                  rmse=round(d3, digits=digits),
                  MAPE=round(d4b, digits=digits),
                  R2=round(d5, digits=digits),
                  row.names = 1))
}

prophet関数のオプション

引数 概要
df (オプション) 履歴を含むデータフレーム。ds (日付型) と y (時系列) の列を持たなければならない。成長がロジスティックである場合、dfは、各dsの容量を指定する列capも持っていなければならない。モデルをフィットさせるにはfit.prophet(m, df)を使ってください。
growth 文字列(String) 'linear', 'logistic', または 'flat' で、線形、ロジスティック、またはフラットトレンドを指定します。
changepoints 潜在的なチェンジポイントを含む日付のベクトル。指定しない場合、変化点候補は自動的に選択されます。
n.changepoints 含めるべき潜在的な変化点の数。入力 'changepoints' が与えられている場合は使用されません。changepoints' が与えられない場合、 df$ds の最初の 'changepoint.range' の割合から一様に n.changepoints 潜在的な変化点が選択されます。
changepoint.range トレンド変化点が推定される履歴の割合。デフォルトでは、最初の 80 'changepoints'が指定された場合、0.8 となります。
yearly.seasonality 年間の季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
weekly.seasonality 週ごとの季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
daily.seasonality 日次の季節性をフィットさせる。auto', TRUE, FALSE, あるいは生成するフーリエ項の数を指定することができます。
holidays holiday (文字) および ds (日付型) の列と、休日として含める日付を指定する lower_window および upper_window の列を持つデータフレームです。また、オプションでprior_scaleカラムを持ち、各祝日の事前スケールを指定することができます。
seasonality.mode 'additive' (default) or 'multiplicative'.
seasonality.prior.scale 季節性モデルの強さを調節するパラメータ。大きな値を設定すると、より大きな季節変動に対応できるようになり、小さな値を設定すると、季節性が弱くなる。個々の季節性に対して add_seasonality を用いて指定することができる。
holidays.prior.scale 休日入力で上書きされない限り、休日成分モデルの強さを調節するパラメータ。
changepoint.prior.scale チェンジポイントの自動選択の柔軟性を調整するパラメータ。大きな値を設定すると多くのチェンジポイントが選択され、小さな値を設定するとチェンジポイントはほとんど選択されません。
mcmc.samples 整数(Integer)、 0より大きい場合,指定された数のMCMCサンプルで完全なベイズ推定を行います.0の場合, MAP推定を行います.
interval.width 数値(Numeric)、予測に提供される不確実性区間の幅。mcmc.samples=0の場合、これは外挿された生成モデルのMAP推定値を使用したトレンドの不確実性のみとなる。mcmc.samples>0の場合、これは季節性の不確実性を含むすべてのモデルパラメータにわたって統合される。
uncertainty.samples 不確実性区間を推定するために使用されるシミュレーション描画の数。この値を0またはFalseに設定すると、不確実性推定が無効となり、計算が高速化されます。
fit true/false, FALSE の場合、モデルは初期化されるが、フィットしない。

make_future_dataframe関数のオプション

引数 概要
m Prophet モデルオブジェクト
periods 将来予測を行う期間数。
freq 'day', 'week', 'month', 'quarter', 'year', 1(1 sec), 60(1 minute) or 3600(1 hour).
include_history true/false; 予測用のデータフレームに過去の日付を含めます。

参考資料

github.com

www.slideshare.net

atmarkit.itmedia.co.jp

www.datatechnotes.com

rstudio.github.io

predictionoffindata.blog.jp