/BDS

Haruo Suzuki (haruo[at]g-language[dot]org)
Last Update: 2015-11-09

"Bioinformatics Data Skills by Vince Buffalo (O’Reilly). Copyright 2015 Vince Buffalo, 978-1-449-36737-4."

Print & Ebook

News & Reviews

Bioinformatics Data Skills: Reproducible and Robust Research With Open Source Tools

バイオインフォマティクス・データスキル：オープンソースのツールによる再現可能で頑強な研究

Table of Contents

I. Ideology: Data Skills for Robust and Reproducible Bioinformatics

II. Prerequisites: Essential Skills for Getting Started with a Bioinformatics Project

III. Practice: Bioinformatics Data Skills

Preface

データスキルは、試練を経たオープンソースのツールを利用するので、同じスキルで次世代のデータにも適応できる。

The Approach of This Book

本書は、ソフトウェアの実行方法は扱わない。
本書は、複雑で大規模なデータセットから意味を抽出し探索する技術を扱う。
本書を通して、頑強（robust）で再現可能（reproducible）な方法を強調する。

Why This Book Focuses on Sequencing Data

本書は、主に配列データを扱う。配列データは豊富にあり、配列データ解析に必要なテキスト処理技術は、他のデータにも適用できる。

Audience

生物学者と計算機科学者の両方を対象

The Difficulty Level of Bioinformatics Data Skills

ハード

Assumptions This Book Makes

前提知識は以下の通り。

Supplementary Material on GitHub

GitHubリポジトリの補足資料を取得する。

Computing Resources and Setup

例題は、Unix系のOS（Mac OS XやLinux）で動作する。 パッケージ管理システム（Ubuntu/Debianのapt-getやMac OS XのHomebrew）を使用する。

Organization of This Book

本書は3部構成：第I部はイデオロギーに関する1章、第II部は基礎編、第III部は実践編。

Part I. Ideology: Data Skills for Robust and Reproducible Bioinformatics

第I部. イデオロギー：頑強で再現可能なバイオインフォマティクスのためのデータスキル

Chapter 1. How to Learn Bioinformatics

第1章. バイオインフォマティクスの学び方

Why Bioinformatics? Biology’s Growing Data

Figure 1-1. DNA配列決定コストの減少 DNA Sequencing Costs

Figure 1-2. Sequence Read Archiveの指数的成長

ショードリード配列マッピング・ツール

Learning Data Skills to Learn Bioinformatics

New Challenges for Reproducible and Robust Research

Reproducible Research

再現可能な研究
データ、コード、ソフトウェアのバージョンとダウンロードした日時を記録する。

Robust Research and the Golden Rule of Bioinformatics

"garbage in, garbage out"「ゴミを入れればゴミが出てくる」

バイオインフォマティクスの黄金律:

NEVER EVER TRUST YOUR TOOLS (OR DATA)

ツールやデータを絶対に信用しない。

Adopting Robust and Reproducible Practices Will Make Your Life Easier, Too

Recommendations for Robust Research

頑強な研究のススメ

Pay Attention to Experimental Design

実験計画

Write Code for Humans, Write Data for Computers

Style guides for Google-originated open-source projects

Let Your Computer Do the Work For You

Make Assertions and Be Loud, in Code and in Your Methods

前提条件のチェックに表明（assertion）を使用する。Pythonのassert()やRのstopifnot()

Test Code, or Better Yet, Let Code Test Code

単体テスト(unit testing)
例えば、add()関数をテストするtest_add()関数をPythonで書く:

Use Existing Libraries Whenever Possible

なるべく既存のライブラリを使う
歴史が長く、閲覧者が多いので、バグが少ない。

Treat Data as Read-Only

データを読み取り専用として扱う
Excelでセルの値を変更して保存するのはダメ。プログラムが、データを読み取り、新しい別の結果ファイルを作成するのが良い。元のファイルを変更してしまうと、再試行・再現できなくなる。

Spend Time Developing Frequently Used Scripts into Tools

Let Data Prove That It’s High Quality

探索的データ解析 (Exploratory Data Analysis; EDA)を通してデータの質を証明する。Chapter 8でR言語を用いてEDAを学ぶ。

Recommendations for Reproducible Research

再現可能な研究のススメ

Release Your Code and Data

データとコードを公開する

Document Everything

全て記録する

プレーンテキスト形式のREADMEファイルに解析の各ステップを記録する。

Make Figures and Statistics the Results of Scripts

図表を出力するスクリプトを書く

Use Code as Documentation

ドキュメントとしてコードを使用する

Continually Improving Your Bioinformatics Data Skills

Part II. Prerequisites: Essential Skills for Getting Started with a Bioinformatics Project

第II部. 前提条件：バイオインフォマティクス・プロジェクト入門のための必須スキル

Chapter 2. Setting Up and Managing a Bioinformatics Project

第2章. バイオインフォマティクス・プロジェクトの作成と管理

Project Directories and Directory Structures

プロジェクト・ディレクトリの構造

プロジェクトの全ファイルを1つのディレクトリに格納する。

例えば、トウモロコシ（学名Zea mays）のSNP検出プロジェクトのディレクトリ（zmays-snps/）を作成する:

  • data/ディレクトリにデータを格納する。
  • scripts/ディレクトリにスクリプトを格納する。
  • analysis/ディレクトリに解析結果を格納する。
What’s in a Name?

ファイル（ディレクトリ）名には、 スペース（空白）を使わない、 英数字かアンダースコアかダッシュ（ A-z a-z 0-9 _ - ）を使う。 拡張子を付ける。（例 osativa-genes.fasta）

スクリプトがプロジェクト・ディレクトリ内の他のファイルを参照する場合には、絶対パス（例 /home/vinceb/projects/zmays-snps/data/stats/qual.txt）ではなく、相対パス（例 ../data/stats/qual.txt）を使う。

Project Documentation

プロジェクトの記録

記録する情報の例は以下の通り。

  • 方法とワークフロー。全コマンドラインをコピー＆ペースト。デフォルト値も
  • データの入手元（ウェブサイトのURL等）
  • データをダウンロードした日付
  • データのバーション（例 TAIR10、WS231）
  • データのダウンロード方法（例 MySQL、UCSC Genome Browser）
  • ソフトウェアのバーション（なければ、日付やURL）

以上の情報をプレーンテキスト形式のREADMEファイルに記録する。プレーンテキストはコマンドラインから簡単に読込・検索・編集できる。

READMEファイルはプロジェクトの主ディレクトリに格納する。

data/READMEファイルに、data/ディレクトリのデータファイルの説明（いつ・どこから・どのようにダウンロードしたのか）を記載する。touchコマンドでサイズが0の空ファイルを作成する:

Use Directories to Divide Up Your Project into Subprojects

プロジェクトをサブプロジェクトに分割するディレクトリを作成

Organizing Data to Automate File Processing Tasks

ファイル処理を自動化するために、データをサブディレクトリに編成し、明確で一貫性のあるファイル名を付ける。

Shell Expansion Tips シェルの展開

cd ~でホームディレクトリに移動。ワイルドカードのアスタリスク（*）は全ての文字列にマッチする。
Brace expansion ブレース展開の例:

zmays-snps/プロジェクト・ディレクトリを作成:

3つのサンプル（zmaysA, zmaysB, zmaysC）毎にペア（R1, R2）の空データファイルを作成する:

ワイルドカードのアスタリスク（*）を用いて、サンプル名zmaysBを持つ全てのファイルを表示する:

Wildcards and "Argument list too long"

偶然の一致を避けるために、ワイルドカードを可能な限り限定する。例えば、zmaysB*の代わりに、zmaysB*fastqまたはzmaysB_R?.fastqを用いる（?は任意の1文字）。

文字列[AB]や文字の範囲[A-B]にマッチするワイルドカードを用いて、サンプルCを排除する:

ワイルドカードは存在するファイルを展開するのに対して、brace expansion（例 snps_{10..13}.txt）はファイルやディレクトリが存在するか否かに関係なく展開する。

Table 2-1. Unixのワイルドカード

Leading Zeros and Sorting

ファイル名の先頭に0を付ける（例 file-21.txtではなくfile-0021.txtにする）と、lsで辞書順にファイルがソートされる。
touch gene-{1..14}.txt
printf "gene-%03d.txt " {1..14} | xargs touch

Markdown for Project Notebooks

プレーンテキスト形式で書かれたプロジェクト・ノートは、コマンドラインやネットワーク経由で読込・検索・編集できる。

Markdown（マークダウン）

Markdown Formatting Basics

John Gruberのホームページ（Daring Fireball: Markdown Syntax Documentation）

Figure 2-1. MarkdownノートのHTML表示

Table 2-2. Markdown記法

見出し（Header）、リスト、コードの書き方

見出しのレベル（1～6）は、#の個数で表す:

リストは、行頭にダッシュ（-）、アスタリスク（*）、プラス（+）か、番号ピリオド（1.）:

コードは、行頭に「半角スペースを4つ」か「タブを1つ」を追加:

リストの項目内にコードを配置する場合、「半角スペースを8つ」か「タブを2つ」にする。

MultiMarkdown
GitHub Flavored Markdown - User Documentation

Using Pandoc to Render Markdown to HTML

Pandocのインストール

Chapter 3. Remedial Unix Shell

第3章. Unixシェル

本章では、ストリーム、リダイレクト、パイプ、プロセス、コマンド置換（command substitution）を扱う。

Why Do We Use Unix in Bioinformatics? Modularity and the Unix Philosophy

UNIX哲学

The Many Unix Shells

bashを使う。
echo $SHELL (echo $0) で現在のシェルを確認
chshでログインシェルを変更

Working with Streams and Redirection

ストリームとリダイレクト

Redirecting Standard Out to a File

標準出力をファイルにリダイレクト

タンパク質のアミノ酸配列データを記述したFASTA形式ファイル tb1-protein.fastaと tga1-protein.fastaを見る。

catコマンドで tb1-protein.fasta ファイルを標準出力する:

複数のファイルを標準出力する:

記号>（上書き）や>>（追記）で標準出力をファイルにリダイレクトする:

Figure 3-1.

ls -lrtで更新日時の逆順にソートする（詳細はman lsを参照）。

Redirecting Standard Error

標準エラー出力をリダイレクト

ls -l tb1.fasta leafy1.fastaを実行すると、存在するファイル（tb1.fasta）は標準出力に、存在しないファイル（leafy1.fasta）は標準エラー出力に送られる。
記号>と2>を用いて、標準出力（standard output）と標準エラー出力（standard error）を別のファイルにリダイレクトする:

記号2>は上書き、2>>は追記。

File Descriptors

ファイル記述子 2>

擬似デバイス（pseudodevice）の /dev/null は、あらゆる入力を受け付けて捨てる。

Using tail -f to Monitor Redirected Standard Error

tail -fでリダイレクトされた標準エラー出力を監視する。Control-Cで動作中のプロセスを停止。

Using Standard Input Redirection

標準入力リダイレクト演算子<よりも、Unixパイプ（例 cat inputfile | program > outputfile）を使う方が一般的。

The Almighty Unix Pipe: Speed and Beauty in One

Figure 3-2.

Pipes in Action: Creating Simple Programs with Grep and Pipes

パイプとgrepを用いて、FASTAファイルに含まれるATGC以外の文字を探す:

ハイライトされたYはpYrimidine塩基[CT]を示す（Nucleic acid notation）。 正規表現はクオーテーションで囲む（例 ">"）。grep -v > tb1.fastaとした場合、シェルは>をリダイレクト演算子と解釈し、ファイルを上書きしてしまう！

Combining Pipes and Redirection

2>&1演算子は標準エラー出力を標準出力にリダイレクトする:

Even More Redirection: A tee in Your Pipe

tee

Managing and Interacting with Processes

プロセス操作の基本：バックグラウンドでプロセスを実行・管理、プロセスを強制終了、プロセスの終了ステータスを確認

Background Processes

コマンドの末尾にアンパサンド（&）を追加して、プログラムをバックグラウンドで実行する:

[1]はジョブ番号、26577はプロセスID（PID）。 jobsでバックグランド・ジョブを表示する:

fgでバックグラウンド・プロセスをフォアグラウンド（foreground）へ戻す。fg %ジョブ番号で指定。バックグラウンド・プロセスが1つしかない場合には、fgとfg %1は同じ:

Background Processes and Hangup Signals

バックグラウンドのプロセスがハングアップ・シグナル（SIGHUP） を受け付けないようにするには、nohupを使うか、Tmux内で走らせる。Chapter 4で詳細に述べる。

Control-z キーで中断させたジョブを bgコマンドを用いてバックグラウンド（background）で再開:

複数の実行中のプロセスがある場合、バックグラウンドで再開させるジョブをbg %ジョブ番号で指定する。

Killing Processes

top、 ps、 killコマンド。 GitHub上の本章のREADMEを参照されたい。

Exit Status: How to Programmatically Tell Whether Your Command Worked

終了ステータス（exit status） 慣習的に正常終了時はゼロ、異常終了時はゼロ以外を返すのが一般的である

Warning Exit Statuses

終了ステータスの値は、シェルの特殊変数$?に設定される。

終了ステータスを判定してコマンドを実行する。

&&は、コマンドが成功した場合のみ次のコマンドを実行する:

||は、コマンドが失敗した場合のみ次のコマンドを実行する:

&&と||をテストするのに、正常終了trueまたは異常終了falseを返すUnixコマンドを用いる:

Command Substitution

コマンド置換 - `command`ではなく$(command)を使う。なぜか？ネストできるから。

date +%Fコマンドを用いて日付ディレクトリを作成する:

このディレクトリ名は年代順にソートされる:

Storing Your Unix Tricks

add aliasを用いて ~/.bashrc（Mac OS Xでは ~/.profile）ファイルに追加する。例えば、常に同じディレクトリ構造のプロジェクト・ディレクトリを作成する:

Chapter 4. Working with Remote Machines

第4章. リモートマシンの操作

Connecting to Remote Machines with SSH

secure shell (SSH)

IPアドレスを用いてマシンに繋げるには ssh 192.169.237.42

ポートとユーザー名を指定する:

ホストに接続できない場合、ssh -vでデバッグする。verboseの-vは-vvや-vvvでより冗長に。詳細は、man ssh。

  • [管理者必見！ ネットワーク・コマンド集 - sshコマンド：ITpro](管理者必見！ ネットワーク・コマンド集 - sshコマンド：ITpro)｜-v デバッグ・モードを有効にする。このオプションは最大三つまで重ねて指定できる。重ねていけば，より詳細な情報が出力される。
  • 入門OpenSSH / 第7章 うまくいかない時は｜ssh をデバッグモードで実行する
Storing Your Frequent SSH Hosts

~/.ssh/configファイルを作る:

Host bio_serv
     HostName 192.168.237.42
     User cdarwin
     Port 50453

リモートホストのデフォルトと同じならPortとUserを指定する必要は無い。ssh -p 50453 cdarwin@192.169.237.42とタイプする代わりに、エイリアスssh bio_servを用いて192.168.236.42にSSH接続できる。

hostnameコマンドは、ホスト名を表示する。
whoamiコマンドは、ユーザー名を表示する。

Quick Authentication with SSH Keys

ssh-keygenコマンドを用いて、秘密鍵（~/.ssh/id_rsa）と公開鍵（~/.ssh/id_rsa.pub）を生成する。

リモートホストにログインして、~/.ssh/authorized_keysに公開鍵ファイル（id_rsa.pub）の内容を追加する。ssh-copy-idコマンドを用いて自動登録も可能。

Maintaining Long-Running Jobs with nohup and tmux

ハングアップシグナル（SIGHUP）
nohupとTmux

nohup

nohupは、プロセスID（PID）を返す。

Working with Remote Machines Through Tmux

terminal multiplexer：TmuxとGNU Screen

Installing and Configuring Tmux

TmuxをMac OS XにHomebrewでインストールする:

設定ファイル（.tmux.conf）をホームディレクトリに置く。シェルが設定を ~/.profile や ~/.bashrc から読み込むように、Tmuxは設定を ~/.tmux.conf から読み込む。

Creating, Detaching, and Attaching Tmux Sessions

Tmuxの新しいセッションを開始する:

Tmuxのプレフィックスキー（キーバインド）は、デフォルトではControl-bを用いるが、設定ファイル（.tmux.conf）でControl-aに変更した。

セッションのデタッチは Control-a d

セッションの確認: tmux list-sessions

セッションのアタッチ: tmux a

Working with Tmux Windows

Tmuxのマニュアルページを見る: man tmux

Table 4-1. Tmuxの基本的キー操作

Table 4-2. Tmuxの基本的コマンド

Emacsで、文字通り Control-a を入力するには、Control-a a とする。

Chapter 5. Git for Scientists

第5章. 科学者のためのGit

Why Git Is Necessary in Bioinformatics Projects

Git Allows You to Keep Snapshots of Your Project

Git Helps You Keep Track of Important Changes to Code

Git Helps Keep Software Organized and Available After People Leave

Installing Git

Gitのインストールは、Mac OS XでHomebrew (brew install git)を、Linuxでapt-get (apt-get install git)を使う。Git website

Basic Git: Creating Repositories, Tracking Files, and Staging and Committing Changes

Git Setup: Telling Git Who You Are

git init and git clone: Creating Repositories

Seqtk (SEQuence ToolKit)

Tracking Files in Git: git add and git status Part I

Staging Files in Git: git add and git status Part II

Figure 5-1.

git commit: Taking a Snapshot of Your Project

テキストエディタを変更

Some Advice on Commit Messages

git commit -a -m "your commit message"

Seeing File Differences: git diff

例えば、README.mdファイルに一行追加して、git diffを実行:

ファイルをステージすると、git diffは何も出力しない。

直近のコミット 比較

Seeing Your Commit History: git log

git log and Your Terminal Pager

クローンしたseqtkリポジトリでgit log

Moving and Removing Files: git mv and git rm

Telling Git What to Ignore: .gitignore

無視させたいファイルを記載した .gitignore ファイルを作成:

.gitignore ファイルをステージし、コミット:

バイオインフォマティクス・プロジェクトで無視させたいファイルの例:

  • 巨大なファイル
  • 中間ファイル（SAMやBAMファイル）
  • テキストエディタ（EmacsやVim）の一時ファイル（例 textfile.txt~ や #textfile.txt#）。.gitignore ではワイルドカード（*~や\#*\#）が使える。
  • 一時コードファイル（Pythonのoverlap.pyc）

Mac OS Xで作成される隠しファイル .DS_Store

グローバルな*.gitignoreファイルを~/.gitignore_global*に作成し、これを使用するようにGitを設定する:

Undoing a Stage: git reset

Collaborating with Git: Git Remotes, git push, and git pull

Figure 5-3. git push (a); git clone (b)
Figure 5-4. git push (a); git pull (b)

Creating a Shared Central Repository with GitHub

the Create a New Repository page
zmays-snps

Authenticating with Git Remotes

Connecting with Git Remotes: git remote

git remote rm <repository-name>

Pushing Commits to a Remote Repository with git push

Pulling Commits from a Remote Repository with git pull

両リポジトリは同じコミットを持つ git logで確認

オリジナルのzmays-snps/ローカル・リポジトリのファイルを修正し、コミットし、pushする:

Barbaraのリポジトリ (zmays-snps-barbara) で、セントラル・リポジトリからpullする。

Barbaraのリポジトリが最新のコミットを含むことをgit logで確認

Working with Your Collaborators: Pushing and Pulling

Merge Conflicts

Figure 5-5

More GitHub Workflows: Forking and Pull Requests

pull request

Using Git to Make Life Easier: Working with Past Commits

Getting Files from the Past: git checkout

Stashing Your Changes: git stash

More git diff: Comparing Commits and Files

SPECIFYING REVISIONS RELATIVE TO HEAD

Undoing and Editing Commits: git commit --amend

Working with Branches

Creating and Working with Branches: git branch and git checkout

隣のアスタリスクは、現在いるブランチを示す。

README.mdを編集:

この変更をreadme-changesブランチにコミットする。masterブランチに戻り、このコミットが存在しないことを確認する:

masterブランチに戻り、adapters.faファイルを追加し、この変更をコミットする:

今、両方のブランチに新しいコミットがある。Figure 5-6.

Merging Branches: git merge

テキスト・エディタが開く。

Branches and Remotes

Continuing Your Git Education

git checkout 変更したファイルを戻す
git stash 修正をいったん退避する
git branch ブランチ操作

Scott ChaconとBen StraubのPro Git book

Chapter 6. Bioinformatics Data

第6章. バイオインフォマティクス・データ

大規模で複雑なゲノム・データの課題:

  • Retrieving data データの取得
  • Ensuring data integrity データ完全性の確保
  • Compression 圧縮

Retrieving Bioinformatics Data

Downloading Data with wget and curl

wgetとcurlは、データをウェブからダウンロードするコマンドラインのプログラム。パッケージ管理システム（Homebrewやapt-get）でインストールできる。

wget

wgetを用いて、GRCh37（hg19）ヒト22番染色体をダウンロードする:

HTTP or FTP の認証は wgetの--user=と--ask-passwordオプションを用いる。

man wgetでオプション一覧を見る。

Table 6-1. wgetのオプション

Curl

curlは、デフォルトでは標準出力するので、wgetと同じようにするには、リダイレクトする（または-O <filename>を使う）:

curlはwgetよりも多くのプロトコル [SFTP (secure FTP) や SCP (secure copy) を含む] を用いてファイルを転送できる。 -L/--locationオプション。 RCurlとpycurlはCurlの機能を提供する（ラッパー）。

Rsync and Secure Copy (scp)

-avzオプションで、-aはwrsyncのアーカイブモード（-rlptgoDと同じ）、-zはデータを圧縮、-vは経過を表示。SSHでリモートホストに繋ぐので、-e sshを用いる。ディレクトリをコピーするコマンドは以下の通り:

末尾のスラッシュを追加するか否か (例えば、data/ と data) で挙動が変わる。

Secure copy (scp)

例えば、GTFファイルを192.168.237.42:/home/deborah/zea_mays/data/に転送する:

Data Integrity

データ完全性

チェックサムで転送データの整合性を検証。

SHA and MD5 Checksums

MD5とSHA-1チェックサム

SHA-1チェックサム。shasum（一部のシステムではsha1sum）プログラムに任意の文字列を渡す:

md5sum（Mac OS Xではmd5）プログラムはMD5ハッシュ値を計算する。

Looking at Differences Between Data

データの違いを見る

diffコマンドで gene-1.bedと gene-2.bed ファイルの差分を出力する:

Compressing Data and Working with Compressed Data

データの圧縮

gzip

gzip

gzipを用いて、trimmer（架空のプログラム）の出力を、ディスクに書き込む前に、圧縮する:

gzipコマンドで圧縮:

gunzipコマンドで解凍:

-cオプションを用いて圧縮・解凍の結果を標準出力に書き出す:

解凍しないで圧縮ファイルに結合する:

Working with Gzipped Compressed Files

圧縮ファイルを直接操作できるコマンドは、zgrep、zcat（Mac OS Xではgzcat）、zdiff、zless

Case Study: Reproducibly Downloading Data

Ensemblのウェブサイトの Mouse の "Download DNA sequence (FASTA)" ftp://ftp.ensembl.org/pub/release-82/fasta/mus_musculus/dna/ を開く。 Genome Reference Consortium GRCm38マウス参照ゲノムをwgetでダウンロードする:

zgrepコマンドを用いて正規表現"^>"で圧縮ファイルのFASTAヘッダを確認:

EnsemblのCHECKSUMSファイルと比較:

SHA-1サムを計算:

GTFとCHECKSUMSをEnsemblからダウンロード:

CHECKSUMSファイルと比較し、SHA-1サムを計算:

Markdownノート（README.md）の例:

Part III. Practice: Bioinformatics Data Skills

第III部. 実践：バイオインフォマティクス・データスキル

Chapter 7. Unix Data Tools

第7章. Unixのデータ・ツール

Unix Data Tools and the Unix One-Liner Approach: Lessons from Programming Pearls

Unixコマンドをパイプで繋ぐことにより、データをパースし操作し集計する1行プログラム（ワンライナー）を構築する。

When to Use the Unix Pipeline Approach and How to Use It Safely

Inspecting and Manipulating Text Data with Unix Tools

タブ区切り

Tabular Plain-Text Data Formats

本章では、3列のBED形式とGTF形式のファイルを用いる。

Inspecting Data with Head and Tail

headでファイルの先頭部分を表示する:

tailでファイルの末尾部分を表示する:

tailでファイルのヘッダを削除する:

ファイルの始まりと終わりの両方を見る:

設定ファイル（~/.bashrcや*~/.profile*）にショートカットを作る:

ショートカットのi（inspect）コマンドを実行:

パイプでgrepの標準出力をheadに渡す:

  • シグナル (ソフトウェア) シグナルSIGPIPEは、読み手のいないパイプへの書き込み
    シグナルSIGINTは、割り込み端末から割り込みキー（通常 CTRL + C）を押下したときに発生

less

lessでcontaminated.fastqファイルを見る:

lessを終了するには、qを押す。hでヘルプを表示する。

Table 7-1. lessの操作方法

lessは、テキスト検索（マッチした部分をハイライト）、パイプラインのデバッグや構築に利用できる。

lessでcontaminated.fastqを開いて、 / を押して、AGATCGGを入力。結果は Figure 7-1

Plain-Text Data Summary Information with wc, ls, and awk

wc（word count）で行数、単語数、文字数を表示:

ls -lでファイルのサイズを確認:

Data Formats and Assumptions

空白（スペース、タブ、改行）を除くには、grepを使う:

awkでファイルの列（フィールド）数を表示:

Mus_musculus.GRCm38.75_chr1.gtfファイルのヘッダを除いてから、列（フィールド）数をawkで表示:

grepを用いて、"#"で始まる行を除く:

Working with Column Data with cut and Columns

cutでタブ区切りファイルの2列目を抽出:

grepでメタデータ行を削除し、cutで1,4,5列（染色体、開始位置、終了位置）を抽出:

cut -dで区切り文字を指定する。CSVファイル:

Formatting Tabular Data with column

タブ区切りファイルの出力は（要素が何列目に属するのか）みにくい:

column -tで整形:

column -sで区切り文字を指定:

The All-Powerful Grep

Figure 7-2
トウモロコシ・ゲノムで文字列"AGATGCATG"を検索した実行時間を、4手法（grep、sed、awk、Pythonスクリプト）間で比較した結果、grepが最速。

ヒト1番染色体の全タンパク質のEnsembl遺伝子識別子と遺伝子名が含まれている Mus_musculus.GRCm38.75_chr1_genes.txtファイル内の遺伝子"Olfr418-ps1"をgrepで見つける:

grep --color=autoオプションでマッチング部分を色付けする。

GNU, BSD, and the Flavors of Grep

GNU coreutilsをMac OS Xにインストール:

grep -vでマッチしない行を返す。"Olfr1413"以外の"Olfr"を含む遺伝子群を出力:

grep -wで（空白で囲まれた）単語全体にマッチする:

パターンにマッチした行の前（-B）、後（-A）、前後（-C）を出力するオプション:

grepは、正規表現 POSIX Basic Regular Expressions (BRE)をサポート。 Ensembl遺伝子識別子"Olfr1413"と"Olfr1411"を見つける:

grep -Eオプションで、POSIX Extended Regular Expressions (ERE) を用いる。 "Olfr218"または"Olfr1416"にマッチ:

grep -cオプションで、パターンにマッチした行数を表示:

grep -oオプションで、マッチした部分だけを抽出:

Example 7-1. ユニークな（重複のない）ソートされた遺伝子名のリストを出力

Decoding Plain-Text Data: hexdump

ASCII。man ascii

何かが正常に動作しない場合、ファイルの文字コードを疑い、file, hexdump, grepコマンドで確認する。

Sorting Plain-Text Data with Sort

Sortで行を並べ替える。

Using Different Delimiters with sort

sort -t（例えば、CSVファイルは-t","）で列の区切り文字を指定する。

sortのオプション。-kで列の範囲（start,end）を指定してソート、-nで数値としてソート。 1列目（染色体 chromosome）でソートし（-k1,1）、1列目が同じもの（例、"chr1"や"chr3"）は2列目で数値としてソートする（-k2,2n）:

Sorting Stability

sort -s

-cオプションで、既にソートされているかチェックする:

-rオプションで逆順（降順）にソートする:

一部のオプションは、GNU sortでのみ利用可能（Mac OS XのBSD版ではない）。例えば、 -Vオプションは英数字を（文字列の中の数字を理解して）ソートする:

merge sortでメモリを超えるファイルをソートできる:

-Sオプションで、K（キロバイト）、M（メガバイト）、G（ギガバイト）、%（総メモリの割合）を指定する。

--parallelオプションで、4つのコアを指定してソートする:

Finding Unique Values in Uniq

uniqは、連続する重複行を削除して出力する:

-iオプションで、大文字と小文字の区別をつけない。
-cオプションで、重複行の数も表示:

Unixコマンド（grep, cut, sort, uniq）を組み合わせて、表形式データの列を要約:

-dオプションで、重複行のみを表示:

Join

example.bedとexample_lengths.txtファイルを使う:

両方のファイルがソートされていない限り、joinは期待通り動作しない:

基本的な構文は、join -1 <file_1_field> -2 <file_2_field> <file_1> <file_2>:

example_lengths.txtにchr3が無い場合:

-aオプションで、ペアにならない行も生成:

Text Processing with Awk

Awk

Gawk versus Awk

Mac OS XはBSD Awk。GNU Awk（Gawk）のマニュアルはman gawk

Awkは、入力データをレコード（行）に分割し、各レコードをフィールドに分割する。レコード全体は変数$0に格納され、フィールドは$1, $2, $3, ...と分割される。

Awkは算術演算子（+, -, *, /, %, ^）をサポートしている。featureの長さ（終了位置 - 開始位置）が18を超える行だけを出力:

Table 7-2. Awkの比較・論理演算子

論理演算子 && (AND), || (OR), ! (NOT) でパターンを繋ぐ。1番染色体上で長さ>10の行を出力:

2番染色体と3番染色体について、長さ（終了位置 - 開始位置）の列を追加する:

BEGINとEND

NR（Number of Record）は現在の行番号

Setting Field, Output Field, and Record Separators

awk -F","でCSVファイルの区切り文字を指定する。

NRを用いて、行の範囲を抽出:

GTFファイル（Mus_musculus.GRCm38.75_chr1.gtf）をBEDファイル（3列）に変換:

Awkの連想配列（associative array）は、Pythonの辞書、Perlのハッシュのように振る舞う。

Table 7-3. Awkの組み込み関数

  • [The GNU Awk User's Guide - 組み込み関数](The GNU Awk User's Guide - 組み込み関数)

Unixコマンド（grep, cut, sort, uniq -c）を用いて、特定の遺伝子の特徴をカウントする:

Bioawk: An Awk for Biological Formats

Bioawkのソースからダウンロード、コンパイル、インストール、またはMac OS Xのパッケージ管理システムHomebrewでインストール:

bioawk -cで入力形式を指定する。サポートされている入力形式（bed, sam, vcf, gff, fastx）と設定変数を見る:

全タンパク質コード遺伝子の長さ（終了位置 - 開始位置）の列を追加:

FASTQをFASTAファイルに変換:

FASTQ/FASTAエントリ数をカウント:

配列の相補鎖を求める:

オプション-c hdr

Stream Editing with Sed

GNU Sed versus BSD Sed

chroms.txtファイルの染色体名を変換（"chrom1" → "chr1"）:

sedの文字列置換の構文: s/pattern/replacement/
gフラグで全ての文字列を置換する: s/pattern/replacement/g
iフラグで大文字と小文字の区別をつけない: s/pattern/ replacement/i

"chr1:28427874-28425431" (染色体名:開始位置-終了位置) を3列で出力:

（head -n 10と同様）ファイルの先頭10行を出力:

20〜50行まで出力:

Advanced Shell Tricks

Subshells

サブシェルでコマンドをグループ化する例:

メタデータのヘッダがあるGTFファイルをソートする:

ストリームを（gzipで圧縮してから）ファイルにリダイレクトする:

Named Pipes and Process Substitution

名前付きパイプをmkfifoコマンドで作成:

プロセス置換

Figure 7-3. プロセス置換

The Unix Philosophy Revisited

Chapter 8. A Rapid Introduction to the R Language

第8章. R言語入門

探索的データ解析 Exploratory Data Analysis (EDA)

参考:

Getting Started with R and RStudio

THE COMPREHENSIVE R ARCHIVE NETWORK (CRAN)

ggplot2パッケージをインストールする:

R言語とRStudio統合開発環境（IDE）をインストールする。

R Language Basics

Simple Calculations in R, Calling Functions, and Getting Help in R

簡単な計算

関数

Table 8-1. 数学関数

SIGNIFICANT DIGITS, PRINT(), AND OPTIONS IN R

round()関数

GETTING HELP IN R

Variables and Assignment

変数の代入

RStudio Assignment Operator Shortcut

global environment

Vectors, Vectorization, and Indexing

ベクトル

vectorization

recycle

Rの演算子（+, -, *, /）や数学関数（例 sqrt(), round(), log()）はベクトルに対応:

indexing

要素ラベル

subsetting

Out-of-Range Indexing

z[負整数ベクトル]は、対応する要素番号の要素を取り除く:

Negative Indexes and the Colon Operator

要素の並べ替え

比較演算子（Table 8-2 例 ==, !=, <, <=, >, >=）を用いて、論理ベクトル（TRUE; FALSE）を作成する:

v[論理ベクトル]は、TRUEの要素に対応した要素を取り出す（Example 8-1）:

論理演算子 & (AND), | (OR), ! (NOT)

Table 8-2. Rの比較演算子と論理演算子

Vector types

Numeric (double), Integer, Character (strings), Logical

Table 8-3. Rのベクトル型

R’S SPECIAL VALUES

(NA, NULL, Inf/-Inf, and NaN)

ベクトルの要素は全て同じ型

Type Coercion in R

Rの型の強制変換規則

Factors and classes in R

因子とクラス

関数factor()で、カテゴリーを要素としたベクトルを作成する:

関数levels()で、グループを確認:

Rには3種類のオブジェクト指向システム（S3, S4, R5のクラス）がある。

Working with and Visualizing Data in R

Spencer et al. (2006) "The influence of recombination on human genetic diversity."のデータDataset_S1.txtは、集団遺伝学の統計値を含む。例えば、塩基多様度（列PiとTheta）、組換え（列Recombination）、ヒトとチンパンジーのゲノム配列一致率（列Divergence）。他の列は、ウィンドウの開始位置と終了位置（列startとend）、シークエンシング深度（列depth）、GC含量（列%GC）など。

Loading Data into R

作業ディレクトリ

Rに読み込む前に、コマンドラインからファイルを見る:

LARGE GENOMICS DATA INTO R: COLCLASSES, COMPRESSION, AND MORE

colClasses vector to "NULL"
nrow in read.delim()
data.tableパッケージのfread()関数を使う。
RパッケージRSQLite
Rの関数は直接gzip圧縮されたファイルを読み込むこともできる。

read.csv()関数でCSVファイルを読み込む:

read.delim()関数でタブ区切りファイルを読み込む:

デフォルトでは、Rの関数read.delim()とread.csv()は、文字列を文字列(character)ではなく因子(factor)に強制変換する。これを無効にするには、引数stringsAsFactors=FALSE（またはas.is）を使う。

Table 8-4. read.csv() と read.delim() の引数

GETTING DATA INTO SHAPE

Table 8-5. 組織毎の遺伝子発現の計数表（wide形式）
Table 8-6. 組織毎の遺伝子発現の計数表（long形式）
reshape2パッケージはデータを変換する関数を提供する。melt()はwideデータをlongデータに変換し、cast()はlongデータをwideデータに変換する。

Exploring and Transforming Dataframes

read.csv()関数はファイルをデータフレームとして読み込む。
データフレームの各列はベクトル。

Creating Dataframes from Scratch

ドル・マーク($)で列を指定できる。

d[row, col]で行と列を指定できる。

Fragile Code and Accessing Rows and Columns in Dataframes

データフレームから単一の列を取り出す際、結果のベクトル化を防ぐには、引数にdrop=FALSEを指定する:

20番染色体のセントロメアの位置は、25,800,000 から 29,700,000 である。セントロメア領域か否か（TRUE/FALSE）を示す列centを、データフレームdに追加する:

塩基多様度の列Piの大きさを変更した新しい列diversityを作成する:

Exploring Data Through Slicing and Dicing: Subsetting Dataframes

Exploring Data Visually with ggplot2 I: Scatterplots and Densities

ggplot2パッケージ、base graphics、latticeパッケージ

ggplot2のオンライン・ドキュメント

library()関数でggplot2パッケージをロードする:

染色体の位置毎の塩基多様度の散布図（Figure 8-1）を作成する。まず、データフレームdに列position（各ウィンドウの中間点）を追加する。次に、データフレームdの列positionと列diversityの散布図をggplot2で作成する:

Figure 8-1. ggplot2の散布図：ヒト20番染色体の位置毎の塩基多様度

例えば、geom_point(), geom_line(), geom_bar(), geom_density(), geom_boxplot()など。aes()関数

Axis Labels, Plot Titles, and Scales

xlab(), ylab(), ggtitle()関数
scale_x_continuous(limits=c(start, end))、関数scale_x_log10()とscale_y_log10()

Example 8-2は、aes()をggplot()に含み、Figure 8-1と全く同じ散布図を作成する:

Example 8-3は、セントロメア領域か否か（列centのTRUE/FALSE）で色分けして、Figure 8-2を作図:

（同じ位置にプロットが重なっている）overplottingを回避するために、透明度（alpha）を調整して、Figure 8-3を作図:

geom_density()を用いて、多様度の密度を見る（Figure 8-4）:

多様度の密度を、セントロメア領域か否か（列centのTRUE/FALSE）で分けて、透明度（alpha）を半分にして、Figure 8-5を作図:

Exploring Data Visually with ggplot2 II: Smoothing

散布図と平滑化曲線を用いて、シークエンシング深度（列depth）とウィンドウのSNP合計数（列total.SNPs）の関係を見る（Figure 8-6）:

デフォルトでは、ggplot2は一般化加法モデル（generalized additive models; GAM）を用いて、平滑化曲線に合わせる。 help(stat_smooth)、geom_smooth(se=FALSE)

散布図と平滑曲線を用いて、シークエンシング深度に及ぼすGC含量の影響を見る（Figure 8-7）:

Binning Data with cut() and Bar Plots with ggplot2

Binning（離散化）

cut()関数でデータをbinに分割する（Example 8-4は、GC含量）:

ggplot2のgeom_bar()

Figure 8-8.
geom_bar()のxが、 因子（例えば、d$binned.GC）の場合には、ggplot2は計数値の棒グラフ（Figure 8-8の左）を作成する。 連続の数値（例えば、d$percent.GC）の場合には、自動的にデータを離散化してヒストグラム（Figure 8-8の右）を作図する。

Figure 8-9. GC含量でグループ分けされたシークエンシング深度（depth）の密度プロット

Finding the Right Bin Width

binwidthの値を 0.05, 0.5, 1, 5, 10 に変化させる。

Merging and Combining Data: Matching Vectors and Merging Dataframes

Rの%in%演算子

RepeatMaskerで発見されたヒトX染色体上のリピート領域のデータ（chrX_rmsk.txt）を読み込む:

列repClassは因子（factor）であることを確認:

複数のリピートのクラス（DNA、LTR、LINE、SINE、Simple_repeat）の行を選択するために、common_repclassベクトルを作成し、%in%を用いる:

上位5つの最も多いリピートのクラスを計算する:

[%in%]演算子は別の関数match()の簡易版

第1のデータセット（motif_recombrates.txt）は、各モチーフの40kb内の全ウィンドウについて、組換え確率（recombination rate）推定値を含む。第2のデータセット（motif_repeats.txt）は、各モチーフが出現するリピートを含む。2つのデータセットをマージして、特定のリピートに及ぼす各モチーフの組換えの局所的な効果を見る。

Creating These Example Datasets

データ生成コードはmotif-example/

rptsデータフレームの列nameをmtfsの列にマージすることが目的。 paste()関数を用いて、2つの列（chrとmotif_start）を1つの文字列に結合する:

列posは、2つのデータセット間の共通の鍵として機能する。
マージする前に、mtfsのモチーフとrptsの対応するエントリの個数を確認する:

マージする前に、match()を用いて、mtfs$posとrpts$posのインデックス・ベクトル（i）を作成する:

match()の結果を（iへの割り当てをスキップして）直接用いる:

merge()でデータを結合（マージ）する:

Using ggplot2 Facets

Figure 8-10. 配列モチーフとの距離と組換え確率

ggplot2のfacet_wrap()を用いて、これらのモチーフを分割する（Figure 8-11）:

Figure 8-11.

ggplot2のfacet_wrap()とfacet_grid()

Figure 8-12

facet_grid()とfacet_wrap()の何れも引数scalesを指定できる。例えば（Figure 8-13）:

More R Data Structures: Lists

データフレームはリスト; is.list(mtfs)

[ ]はリストを取り出す。 [[ ]]はリスト内の要素（ベクトル）を取り出す。

Peeking into R’s Structures with str()

Writing and Applying Functions to Lists with lapply() and sapply()

Using lapply()

lapply() in Parallel

引数を渡す:

Writing functions

無名関数(anonymous function)

meanRemoveNAVerbose()関数の定義:

Function Scope

スコープ

Digression: Debugging R Code

デバッグ
関数browser()

nで次の行を確認、cでコードの実行を継続、Qで終了

options(error=recover)も有用

More list apply functions: sapply() and mapply()

Other Apply Functions for Other R Data Structures

配列と行列 apply()とsweep()

Working with the Split-Apply-Combine Pattern

データをグループ化し、グループ毎に関数を適用し、結果を組み合わせる（"The Split-Apply-Combine Strategy for Data Analysis"）。最初にRの標準関数を用いて、次にdplyrパッケージを用いる。

split-apply-combineの簡単な例：GC含量でグループ化した3群の平均深度（Example 8-4, Figure 8-9）

Understanding do.call()

関数tapply()とaggregate()でグループ毎に要約する:

Exploring Dataframes with dplyr

dplyrは非常に高速。 dplyrでデータフレームを操作する関数は、select(), filter(), arrange(), mutate(), summarize()

tbl_df()関数を用いて、dデータフレームをtbl_dfオブジェクトに変換する:

selectは、d[, c("start", "end", "Pi", "Recombination", "depth")]に対応:

filter()は、d[d$Pi > 16 & d$percent.GC > 80, ]に対応:

arrange()は、d[order(d$percent.GC), ]に対応:

mutate()関数を用いて、データフレームに新しい列を追加できる。

パイプ（help('%>%')）を用いて、複数の処理を連結する:

mtfsデータフレームを用いる

summarize()関数を用いて要約を作成する:

新たに作成した要約の列max_recomでソートする:

Working with Strings

Rの文字列処理機能

nchar()で文字ベクトルの各要素の文字数を取得する:

grep()かregexpr()で文字ベクトル中のパターンを検索する。関数grep(pattern, x)は、patternにマッチするベクトルxの全要素の番号を返す:

match()と同様、grep()で一貫性のない染色体名のベクトルから6番染色体のエントリを抽出する:

Rの正規表現については help(regex)

The Double Backslash

grep()と異なり、regexpr(pattern, text)は、ベクトルtextの各要素でpatternにマッチした位置を返し、マッチしない場合には-1を返す:

返り値 5 は該当文字列の第5文字目以降にマッチしたことを示す。属性（attributes）"match.length"の値 2 は2文字分マッチしたことを示す。

substr(x, start, stop) は文字列xのstartとstopの間の文字を返す。

sub(pattern, replacement, x)は文字ベクトルxの各要素で最初に出現したpatternをreplacementで置換する:

いくつかの簡単な例:

Friendly Functions for Loud Code

paste()関数は、文字列を結合する:

Extracting Multiple Values from a String

sub()とstrsplit()を組み合わせる:

strsplit(x, split)は、文字列xをsplitで分割する:

Developing Workflows with R Scripts

Control Flow: if, for, and while

Iterating over Vectors

1:length(x) の代わりに seq_along(x) を使うと良いってごみ箱が言ってた - My Life as a Mock Quant

ifelse関数:

Working with R Scripts

Rスクリプトを用いた作業

スクリプトでは、作業ディレクトリを設定するsetwd()を使用しない。（同じディレクトリ構造を持っていない）他のシステムへの移植性が無くなるので。同じ理由で、データを読み込む時には、絶対パス（例 /Users/jlebowski/data/achievers.txt）ではなく、相対パス（例 data/achievers.txt）を使う。また、ユーザー向けに（コメントやREADMEファイルで）作業ディレクトリを指定するのが良い。

関数source()を用いて、Rスクリプトを実行する:

コマンドラインからバッチモードでRスクリプトを実行する:

コマンドライン引数を出力するRスクリプト:

を実行する:

Reproducibility with Knitr and Rmarkdown

Rmarkdown:

Reproducibility and sessionInfo()

Rのバージョンとパッケージを確認: sessionInfo()

Workflows for Loading and Combining Multiple Files

複数のファイルを読み込み、統合する。

hotspots/ディレクトリ:

Rの関数list.files()を用いて、ファイル名を取得する:

lapply()関数で各ファイルを読み込む:

ファイル名（フルパスなし）を用いて、各リストの要素に名前を付ける:

do.call()とrbindを用いて、このデータをマージする:

我々のloadFile関数を変更して、データに適用する:

basename()関数でファイルパスからファイル名を取得:

染色体ファイル毎にホットスポットの数と平均長を求める関数を作成して、実行する:

このデータを1つのデータフレームにマージする:

lapply()をmclapply()で置き換えることにより、データ処理を並列化できる。

Exporting Data

データフレームmtfsを hot‐spot_motifs.txt という名前のタブ区切りファイルに書く:

gzip圧縮をかけて出力する:

serialization
Rオブジェクトを保存・復元する関数はsave()とload()

save, save.image(), savehistory()

Further R Directions and Resources

Chapter 9. Working with Range Data

第9章. 範囲データの操作

Chapter 10. Working with Sequence Data

第10章. 配列データの操作

The FASTA Format

FASTA
FASTA形式ファイルは、">"で始まるコメント行と、配列データが記述される。

一般的な命名規則は、コメント行をスペースで2つの部分（IDと説明）に分割する:

The FASTQ Format

Fastq
NGS Surfer's Wiki | FASTQ

FASTQファイルは、以下のようなものである:

FASTQファイル内では、1本の配列は4行で記述される。
1行目：文字「@」で始まり、配列のID（オプションとして説明）を記述する。
2行目：塩基配列
3行目：文字「+」
4行目：クオリティ値 quality score（2行目の配列と同じ文字数でなければならない。ASCIIコードで表現される）

The Ins and Outs of Counting FASTA/FASTQ Entries

FASTA/FASTQエントリの計数

@は、クオリティ値として行頭にくることもある。

配列数を計数する頑強な方法はbioawkを用いる:

Nucleotide Codes

核酸コード
A、T、C、Gは、ヌクレオチドのアデニン、チミン、シトシン、グアニンを表す。

Table 10-1. Nucleotide base codes (IUPAC)

Base Qualities

塩基のクオリティ値 quality scoreは、ASCIIコードで表現される。man ascii

ASCIIコードを変換する関数（整数から1文字へ、1文字から整数へ）。Pythonの関数ord()とchr()。例えば、chr(97)は文字列 'a' を返し、ord('a')は整数 97 を返す。2. 組み込み関数 — Python 2.7ja1 documentation

Example: Inspecting and Trimming Low-Quality Bases

Rにqrqcをインストール:

sickleとseqtkをHomebrewを用いてMac OS Xにインストール:

FASTQファイル untreated1_chr4.fq をトリムする:

これらの結果をRで比較する。qrqcを用いて位置毎のクオリティの分布を収集し、ggplot2を用いて可視化する。lapplyで自動化する。

このRスクリプトは2つのプロット（Figures 10-1と10-2）を作成する。

A FASTA/FASTQ Parsing Example: Counting Nucleotides

FASTA/FASTQのパースには既存のライブラリを使うのがベスト。
readfqのPythonスクリプト（readfq.py）をダウンロードする。
from readfq import readfq

Indexed FASTA Files

Samtoolsのfaidxサブコマンドを用いて、FASTAファイルのインデックスを作成する:

インデックス・ファイル（Mus_musculus.GRCm38.75.dna.chromosome.8.fa.fai）が作成される。

特定の領域の部分配列にアクセスするには、samtools faidx <in.fa> <region>を実行する。ここで、<in.fa>は（インデックスを作成した）FASTAファイル、<region>はchromosome:start-endの形式。例えば:

Chapter 11. Working with Alignment Data

第11章. アライメントデータの操作

Getting to Know Alignment Formats: SAM and BAM

The SAM Header

Read Groups

メタデータ

samtools

samtools view -HでSAM/BAMヘッダ全体を見る :

引数なしのsamtools viewは、ヘッダなしでアラインメント全体を返す:

The SAM Alignment Section

SAMファイルのアラインメント部分は11フィールド以上から成る。

Bitwise Flags

CIGAR Strings

Mapping Qualities

Command-Line Tools for Working with Alignments in the SAM Format

Using samtools view to Convert between SAM and BAM

Samtools Sort and Index

Visualizing Alignments with samtools tview and the Integrated Genomics Viewer

Pileups with samtools pileup, Variant Calling, and Base Alignment Quality

Chapter 12. Bioinformatics Shell Scripting, Writing Pipelines, and Parallelizing Tasks

第12章. バイオインフォマティクスのシェルスクリプト、パイプライン、タスクの並列化

頑強で再現可能なパイプラインを構築する

Basic Bash Scripting

Writing and Running Robust Bash Scripts

A robust Bash header

template.sh 頑強なBashスクリプトのヘッダ:

  • 1行目：shebangは、スクリプトを実行するインタプリタを指定する。
  • 2行目：set -eは、異常終了時（非ゼロ終了ステータス）スクリプトを終了させる。
  • 3行目：set -uは、変数の値が設定されていない場合にスクリプトを終了させる。echo "rm $NOTSET/blah"
  • 4行目：set -o pipefailは、パイプで繋いだコマンドの何れかが非ゼロ終了ステータスを返したら、スクリプトを終了させる。

短縮してset -euo pipefailと書ける。

Running Bash scripts

Bashスクリプトを実行する方法:

  1. bashプログラムを用いる: bash script.sh
  2. プログラムとしてスクリプトを呼び出す: ./script.sh
    chmodコマンドでファイルの所有者（u）に実行権限を追加する（+x）: chmod u+x script.sh

Variables and Command Arguments

変数に値を割り当てる（=の前後にスペースを入れない）:

変数の値にアクセスするためには、変数名の前にドル記号を付ける（例 $results_dir）:

中括弧{}で変数名を囲む:

ダブルクォーテーション""で変数を囲む:

Command-line arguments

コマンドライン引数は、$1, $2, $3, ...に割り当てられる。変数$0はスクリプト名を格納する。

このファイルを実行すると、割り当てられた引数（$0, $1, $2, $3）を出力する:

変数$#にはコマンドライン引数の個数を割り当てる（スクリプト名$0は引数としてカウントしない）:

3未満の引数を指定して、このスクリプトを実行するとエラーになる（非ゼロ終了ステータス）:

Pythonのargparseモジュールは、Unixツールgetoptよりも簡単に使える。

Reproducibility and Environmental Variables

exportコマンドで環境変数

なお、some_var=3で変数を作成するスクリプトを実行しても、現在のシェルにsome_varは作成されない。

Conditionals in a Bash Script: if Statements

bashの条件文：if文

Bashでは、コマンドの終了ステータスが 真/成功 true/success (0) と 偽/失敗 false/failure (1) を与える。

if文の基本構文:

例えば、特定の文字列がファイルに含まれる（grepで"pattern"に一致する）場合にのみコマンドを実行する:

if文で&&（論理積AND）や||（論理和OR）演算子を用いる:

!で終了ステータスを否定する:

testコマンドの実行例（echo "$?"で終了ステータスを出力）:

Table 12-1. 文字列と整数の比較演算子

Table 12-2. ファイルとディレクトリのテスト式

if test -f some_file.txtをif [ -f some_file.txt ] で代用できる。角括弧[ ]の前後に半角スペースが必要。
この構文では、-a（論理積AND）、-o（論理和OR）、!（否定）を使える。&&と||演算子はtestでは使えない。 例として、スクリプトが十分な引数を持ち、入力ファイルが読み込み可能であることを確認する:

短絡評価（short-circuit evaluation）

Processing Files with Bash Using for Loops and Globbing

for文とglobでファイル処理

bash 配列、要素、添字

コマンド置換 (command substitution)

The Internal Field Separator, Word Splitting, and Filenames

sample_files=($(cut -f 3 samples.txt)) Internal Field Separator (IFS)の値を調べる: printf %q "$IFS"

Unixプログラムbasenameは、ファイル名からパスや拡張子を削除する:

処理スクリプト:

最後に、ループで簡単に:

Automating File-Processing with find and xargs

Using find and xargs

Finding Files with find

lsとは異なり、findは再帰的に検索する。

findでプロジェクト・ディレクトリの構造を見る:

findの基本構文は、find path expression

Example 12-1. ファイル名の一致で検索

find’s Expressions

（名前が一致するディレクトリではなく）FASTQファイルのみを返すようにしたいので、-typeオプションで結果を制限する（fはファイル、dディレクトリ、lはリンク）:

findコマンドの論理演算子

AND（-and）:

OR（-or）:

Table 12-3. findの判別式と演算子

NOT（-notまたは!）:

seqs/zmaysB_R1-temp.fastqファイルを作成:

find’s -exec: Running Commands on find’s Results

findでマッチしたファイルを-exec command ;で処理する。rm -iでファイルの削除時に問い合わせを行う。

Deleting Files with find -exec:

-deleteオプション

xargs: A Unix Powertool

xargs

Playing It Safe with find and xargs

ファイル名には、英数字かアンダースコアかダッシュを使い、スペースや他の特殊文字を使わない。

xargs -nで、コマンドライン1つにつき使用する引数の個数を指定:

rmを実行する前に、findが返すファイルを確認する:

xargsを用いてコマンドを書いたスクリプトを作成し、確認してから実行する:

Using xargs with Replacement Strings to Apply Commands to Files

あるプログラムはオプションで引数を取る（例 program --in file.txt --out-file out.txt）、別のプログラムは位置指定引数を取る（例 program arg1 arg2）。xargsの-I{}オプション

BSD and GNU xargs

Mac OS Xでは、Homebrewを用いて、GNU Coreutilsをインストール

xargs and Parallelization

xargsのオプション-P <num>で、プロセスまで同時に実行する

xargs, Pipes, and Redirects

このスクリプトを実行する:
find . -name "*.fastq" | xargs -n 1 -P 4 bash script.sh
ここで、-n 1はコマンドライン1つにつき最大1個の引数を使用することを意味し、-Pオプションで同時に実行するプロセスの個数を指定できる（0ならできるだけ多くのプロセスを同時に実行しようとする）。

GNU Parallel

Make and Makefiles: Another Option for Pipelines

Make

宣言型プログラミング（英: Declarative programming）
各ルールは3つのパート（ターゲット、必要条件、レシピ）を持つ。

Makefile
をmake allで実行する

Chapter 13. Out-of-Memory Approaches: Tabix and SQLite

第13章. TabixとSQLite

Fast Access to Indexed Tab-Delimited Files with BGZF and Tabix

BGZF (Blocked GNU Zip Format)
bgzipとtabixプログラムは、Samtools (例 samtools-1.2/htslib-1.2.1) に含まれている。

Compressing Files for Tabix with Bgzip

ファイルの先頭にメタデータのヘッダがある。

サブシェルを使う。gzipの代わりにbgzipで:

Indexing Files with Tabix

tabixの-p引数を用いて、bgzipで圧縮されたGTFファイルにインデックスを張る:
インデックスファイルの拡張子は*.tbi*:

Using Tabix

例えば、Mus_musculus.GRCm38.75.gtf.bgz の16番染色体上の23,146,536（開始位置）から23,158,028（終了位置）までのフィーチャーにアクセスする:

このコマンドの標準出力をパイプでawkに渡し、列featureが「exon」の行を抽出する:

Introducing Relational Databases Through SQLite

flat file フラットファイル形式。UnixツールのjoinやR言語のmatch()やmerge()関数を用いてテーブルを結合できる。relational databases 関係データベース。

SQL (Structured Query Language)。

relational database management system (RDBMS) 関係データベース管理システムのSQLite。

When to Use Relational Databases in Bioinformatics

アプリケーションプログラミングインタフェース（API）

Installing SQLite

Homebrew（brew install sqlite）でMac OS Xにインストール、またはapt-get install sqlite3でUbuntuマシンにインストール。

Exploring SQLite Databases with the Command-Line Interface

SQLiteのコマンドラインツールsqlite3を用いて、データベースgwascat.dbファイルに接続すると、対話的（インタラクティブ）なSQLiteプロンプトが出る:

ドット・コマンドは.で始まる（空白で始めることはできない）。

Table 13-1. SQLite3のドットコマンド

SQLiteでは、列は型を持たないが、データ値は型を持つ。データ値は5タイプの何れか: text、integer、real、NULL、BLOB（binary large object）

Orderly Columns

Querying Out Data: The Almighty SELECT Command

SELECT文を用いて、テーブルの全ての列から全ての行を取得する:

Working with the SQLite Command-Line Tool

Control-u で、入力をクリアする。

sqlite3のコマンドラインツールは、（対話的なSQLiteのシェル代わりに）直接コマンドラインから問い合わせ可能。例えば、gwascatテーブル内の全てのデータを取得する:

列の分離形式（オプション-separatorに、CSVは","を、タブは"\t"を指定）や、列ヘッダの表示（-header）・省略（-noheader）を変更:

Limiting results with LIMIT

Selecting columns with SELECT

カンマ区切りで一部の列（trait, chrom, position, strongest_risk_snp, pvalue）を指定:

SQLiteの設定を変更:

Ordering rows with ORDER BY

SELECTが返す行はソートされていない。例えば、著者の姓（author）でソートして、研究に関連する列（author、trait、journal）を表示する:

ORDER BY句にDESCを追加すると、降順で結果を返す:

SQLiteの列は厳格な型を持たない。データ値の型が混在する列に、ORDER BYを使うと、次の順に従う。NULL値 > 整数と実数の値（数値でソート） > テキスト値 > BLOB値。ソートすると、NULL値が常に最初にくる。p値で昇順ソートして確認:

第8章で議論したように、データを並べ替えて異常値を探すことは、データの問題を探す有効な方法である。例えば、p値（0〜1の範囲を示す確率）を降順にソートして確認:

2つの誤ったp値はデータ入力ミス。p値を入力するときにマイナス記号を忘れる（例えば、9e-7 を 9e7 と記載）。

Filtering which rows with WHERE

WHERE句で特定の行を除外。例えば、strongest risk SNP = "rs429358"の行を表示:

大文字と小文字は区別される。比較の前にlower()関数で値を小文字に変換:

Table 13-2. WHERE句で使用される演算子

論理演算子ANDとORを用いて条件式を組み合わせる。例えば、例えば、第22染色体上でp-valueが10e-15未満の全ての行を取得:

WHERE pvalue IS NOT NULLでNULL値を排除してから並べ替える:

ORで何れかの条件を満たす行を選択する:

ORと=の代わりに、IN（またはNOT IN）を用いる:

The Human Genome Browser at UCSC

SQLite Functions

既存の列から新しい列を作成する。

||は文字列を連結する演算子。

全てのNULL値をNAに置き換えるには、ifnull()関数を用いる:

Table 13-3. SQLiteの関数

SQLite Aggregate Functions

  • count関数 - SQLite入門
    引数にはカラム名または「*」を指定。カラム名を指定した場合にはカラムに含まれる値の中でNULLのカラムを除いた行数を返す。「*」を指定した場合には行数を返す。

Table 13-4. SQLiteの集計関数

YYYY-MM-DD。xkcd: ISO 8601

Grouping rows with GROUP BY

Subqueries

Organizing Relational Databases and Joins

Organizing relational databases
Inner joins
Left outer joins
Other Types of Outer Joins

Writing to Databases

Chapter 14. Conclusion

Glossary

Bibliography

Index

オンライン教材