2015.07.27版

配布するUSBメモリ中のhogeフォルダを
デスクトップにコピーしておいてください。

NGSハンズオン
講習会：R基礎東京大学・大学院農学生命科学研究科
アグリバイオインフォマティクス教育研究プログラム
門田幸二（かどた こうじ）
kadota@iu.a.u-tokyo.ac.jp
http://www.iu.a.u-tokyo.ac.jp/~kadota/
Jul 29 2015

1

Contents（全体）















7月22日（水）：84→83名。Bio-Linux 8とRのインストール状況確認。基本自習（門田・寺田先生）
7月23日（木）：92→90名。Linux基礎。LinuxコマンドなどUNIXの基礎の理解（門田）
7月24日（金）：85→83名。スクリプト言語。シェルスクリプト（アメリエフ株式会社 服部恵美先生）
7月27日（月）：93→91名。スクリプト言語。Perl（アメリエフ 服部先生）
7月28日（火）：91→90名。スクリプト言語。Python（アメリエフ 服部先生）
7月29日（水）：94→88名。データ解析環境R（門田） R基礎1(初級)：R言語の基礎（インストールから利用まで）
 R基礎2(初級)：ファイルの読み込み、行列演算の基本
 R各種パッケージ(中級)：パッケージのインストール法と代表的なパッケージの利用法
7月30日（木）：96→91名。データ解析環境R（門田） Bioconductorの利用法1(中級)：データの型やバージョンの違い
 Bioconductorの利用法2(中級)：FASTA/FASTQファイルの各種解析
8月3日（月）：89→84名。NGS解析。基礎（アメリエフ 山口昌雄先生）
8月4日（火）：85→80名。NGS解析。ゲノムReseq、変異解析（アメリエフ 山口先生）
8月5日（水）：86 →81名。NGS解析。RNA-seq、統計解析（前半：山口先生、後半：門田）
8月6日（木）：104 →98名。NGS解析。ChIP-seq（理研 森岡勝樹先生）

Jul 29 20152

各種ソフトの場所④①②Excelは行列データファイルの確認用。門田はEmEditorというテキストエディ
タを使っています。③「受講生の心構え」
でも書いていますが、貸与PCのほとん
どはR ver. 3.1.2, 3.1.3, 3.2.0, 3.2.1のいず
れか（または複数）がインストールされています。基本的には最新版を利用。④エ
ディタはR付属のものを推奨。主目的は二重クォーテーション問題の回避。

②①

Jul 29 2015

3

Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 2015

4

おさらい

hogeフォルダ中のr_seq.htmlをダブルクリックしてローカルで利用するのがいいかもしれません。ここで示すようなクリック
して眺めるだけのネットサーフィン系の部
分は、手を動かさずに前のスライドを見ているだけのほうがいいかもしれません。
①

Jul 29 2015

5

おさらい

①「基本的な利用法」PDF中の「解析基礎2」をおさらいします①

Jul 29 20156

解析基礎2

目的：アノテーションファイル（annotation.txt）中の第1列目に対して、リストファイル（genelist1.txt）中の文字列と一致する行を抜き出して、hoge1.txtというファイル名で出力したい

入力1：アノテーションファイル（annotation.txt）

出力：hoge1.txt

入力2：リストファイル（genelist1.txt）

Jul 29 20157

解析基礎2

Jul 29 2015

目的：アノテーションファイル（annotation.txt）中の第1列目に対して、リストファイル（genelist1.txt）中の文字列と一致する行を抜き出して、hoge1.txtというファイル名で出力したい

8

解析基礎2

Jul 29 2015

作業ディレクトリは「デスクトップ – hoge」。
hogeフォルダ中にannotation.txtとgenelist1.txtが存在するという前提。貸与
PCは黒矢印部分が「kadota」ではなく「iu」。

9

基本はコピペ①一連のコマンド群をコピーして②R
作業ディレクトリは「デスクトッ
Console画
面上でペースト。ブラウザがInternet
プ – hoge」。hogeフォルダ中
Explorerの場合は、CTRLとALTキーを押しながらコードのにannotation.txtとgenelist1.txt
枠内で左クリックすると、全選択できます。
が存在するという前提。

①②

Jul 29 2015

10

実行結果

①出力ファイル名として指定したhoge1.txtが生成されているのがわかる。「list.files()で表示される結果」と「実行後のhogeフォルダの中身」は当然同じです
実行前のhogeフォルダ①
実行後のhogeフォルダ①

Jul 29 2015

11

実行結果

①outというオブジェクトの中身をwrite.tableという関数でファイルに出力しています。この場合、出力ファイルhoge1.txtの中身は、Rコンソール
画面中でoutと打ち込むことで見られる。

①

実行後のhogeフォルダ

Jul 29 201512

Rコード中の色の使い分けについて説明します

色の説明

Jul 29 2015

13

応用

Jul 29 2015

このサンプルコードは1列目でキーワード検索する場合。別のリストファイルを読み込んで4列目で検索したい場合のやり方を示します。

14

解答例
1.
2.

「メモ帳」など任意のエディタでリストファイル(list.txt)を作成

目的のキーワードリストを含むファイルを作成し（例：list.txt）
該当箇所を変更し、Rコンソール画面上でコピペ

Jul 29 2015

15

解答例
1.
2.

一連の作業手順を記述したスクリプトを1つのファイルとして保存することをお勧め目的のキーワードリストを含むファイルを作成し（例：list.txt）
該当箇所を変更し、Rコンソール画面上でコピペ

Jul 29 201516

ありがちなミス1

Jul 29 2015

作業ディレクトリの変更を忘れているため、in_f1で指定した最初のファイルの読み込み段階でエラーが出る。つまり、実際に行ったフォルダ中にはannotation.txtというファイルは存在しないということ。

17

ありがちなミス2

Jul 29 2015

必要な入力ファイルが作業ディレクトリ中に存在しない。この場合、in_f2で指定したgenelist1.txtが存在しないため、それの読み込み段階でエラーが出ている。それゆえ、その情報を用いているコマンド部分でエラーが出ている。

18

ありがちなミス3

Jul 29 2015

出力予定のファイル名と同じものをエクセ
ルなど別のプログラムで開いているため、最後のwrite.table関数のところでエラーが出る。対処法は、出力ファイル名を変更す
るか、開いている別のプログラムを閉じる。

19

ありがちなミス4

Jul 29 2015

実行スクリプトをコピーする際、最後の行のところで改行を含ませずにR Console画面上でペーストしたため、最後のコマンドが実行されない（出力ファイルが生成されない）。これも比較的ありがちなパターンです。コピペ後に無意識にリター
ンキーを押すことを心がけるだけでもよいでしょう。

20

警告メッセージ①

Jul 29 2015

list.txtファイル作成時に、membraneと打った後に改行を①入れた場合と②入れない場合の挙動の違いを把握し、後学のために警告メッセージの意味を理解しておくとよい。この場合は結果には影響していな
いことがわかる。Rは警告メッセージの記述内容が比較的分かりやすいのでよく読むべし。

②

21

Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 201522

コード内部の説明

Jul 29 2015

コードの中身を説明します。
黒枠部分を再度コピペ。

23

①in_f1で指定したファイルを読み込め②読み込むファイルの最初の行はヘッダー部分③ファイルの区切り文字はタブです④読み込んだ結果をdataという名前で取り扱う

読み込み

④

Jul 29 2015

①

②

③

24

行列data②

dataと打ってリターン。入力ファイルの中身を正しく読み込めていることがわかる。
②header=TRUEとしているので、③このように見えて列名として認識される。

③

Jul 29 201525

dimで行数と列数を表示②

Jul 29 2015

①オブジェクトdataの行数と列数は11と4。
②ウェブページ中の表記が灰色なのは、特にやらなくてもいいコマンドだから。

①

26

行列の要素へのアクセス

行列dataの要素へのアクセスは[行, 列]。
①humeiは、読み込み元ファイルのannotation.txt中では7行×4列目だが、②1行目をヘッダー行としているので③6行×4列目とする必要がある。利用例
は、ファイル読み込み時に「x行×y列目に不具合がある」のようなエラーが出た時のトラブルシューティングなど。

②③①

Jul 29 201527

Tips：上下左右の矢印キー

Jul 29 2015

上矢印キーを押すと、直前に打ったコマンドが表示される。最初から全部打ち直すのではなく、上下左右の矢印キーを有効に利用し最小限の労力で打つべし!

28

行列の要素へのアクセス

行列dataの要素へのアクセスは[行, 列]。
①2行目の情報のみ抽出。読み込み時にhead=TRUEとしていたので、ヘッダー行がついていることが分かる。

①②

Jul 29 201529

行列の要素へのアクセス

行列dataの要素へのアクセスは[行, 列]。①2列目の情報のみ抽出。

①

Jul 29 201530

行列dataの要素へのアクセスは[行, 列]。①param列目の情報のみ抽出。②paramには1という数値が代入されていたのでこうなる。

行列の要素へのアクセス

①

②

②

Jul 29 201531

参考 行列dataの最初の数行を表示したい場合
は、head関数を利用。①n=3というオプションを利用すると最初の3行分のみ表示。関
数ごとに様々なオプションを利用可能です。
このあたりは②Linuxとよく似ている。

Tips：関数とオプション①②

Jul 29 2015

32

参考

Tips：タブ補完

列番号を指定する以外にも特定の列を表示するやり方がある。head=TRUEで入力
ファイルを読み込むと、列の名前を利用することができる。①subcellular_location列の情報を抽出したい場合は、②「data$su」
くらいまで打ち込んでからTabキーを押す。

①②

Jul 29 2015

33

参考

Tips：タブ補完

列名中のsuからはじまる文字列を補完して表示してくれる。「Tabキーを用いた補完機能」という意味で「タブ補完」という。このテクニックはLinuxでも利用可能。

①②

Jul 29 2015

34

参考

Tips：table関数

Jul 29 2015

table関数は、ベクトル中の要素ごとの出現回数を返す。「NGSデータ中の特定のリードの出現回数(後述)」や、「アノテーションファイル中の染色体ごとの遺伝子数」など、様々な局面で利用可能。

35

参考

Tips：ソート

sort関数と併用することで全体像を俯瞰可能。例えば①nuclearに局在する遺伝子数が最も多く4個であった、などが簡単にわかる。

①

Jul 29 201536

Tips：is.element関数

Jul 29 2015

hogeベクトルに対して、”nuclear”の文字が存在する場所をTRUE、存在しない場所を
FALSEとして返す。as.character関数は、文
字列ベクトルとして取り扱いたい場合に利用。

37

Tips：“二重クォーテーション”

Jul 29 2015

二重クォーテーションが自動で変更されるエディタは非推奨です。日本語の二重クォーテー
ションもだめです。Microsoft
WordやPDFファイル中のコードのコピペ時によくハマります。

38

目的をおさらいJul 29 2015

目的は、数万～数百万行からなるファ
イルを読み込んで特定のキーワードを含む行のみ取り出すテクニックを習得。

39

目的をおさらい

論理値ベクトルobjを用いてTRUEの要素に対応する行を抽出している

入力2：リストファイル（genelist1.txt）

Jul 29 201540

目的をおさらいJul 29 2015

コード作成当時はas.character関数を用いてデータの型
を文字列ベクトルに揃えていた。少なくとも現在(R ver.
3.1.3以降)は、この関数がなくても大丈夫なようだ。同じ関数でもバージョンによって挙動が異なるということ(バージョンの違いの一例)。

41

1と12は手順が異なるだ
けで実質的に同じですgenelist1.txt

Jul 29 201542

このコードはヘッダー行がある場合のものです

入力: annotation.txt
出力: hoge12.txt

Jul 29 201543

このコードはヘッダー行がない場合のものです

入力: annotation2.txt

出力: hoge13.txt

Jul 29 201544

ヘッダー行がある場合

ヘッダー行がない場合

Jul 29 201545

Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 201546

教科書p81-89

カウントデータ
教科書
目的サンプルのRNA-Seqデータ
mapping

リファレンス配列：ゲノム

遺伝子1

Jul 29 2015

カウントデータとは、「マップされたリード数」をカウントしたデータのこと。以下の例では1サンプルなので1列分のデータしかないが、一般には複数サンプルのデータを取得し、サンプル間比較が行われるので複数の列からなる。それゆえ、
数値ベクトルではなく数値行列。詳細は8/5のRNA-seq前半で。

遺伝子2

count

遺伝子3

遺伝子4

47

数値行列
①

実験の詳細には立ち入らないが、3生物種間比較を行った公共RNA-seqカウントデータ(Blekhman et al.,Genome Res., 2010)を用いて、Rの王道的な使い方である数値行列解析のテクニックを伝授。8/5の統計解
析のところでこのデータを利用予定です。

②

Jul 29 2015

48

xls形式ファイルもOK①

Jul 29 2015

①xls形式のエクセルファイルを読み込むことができる。（但し、このファイルは壊れているなどというメッセージが出ており、実はタブ区切りテキス
トファイルなのに、.xlsという拡張子が無理やりつ
けられているというオチかもしれない…）②それほど大きなサイズでなければ、ネットワーク経由
で直接読み込むこともできる。他に、read.csvやreadLines関数などを駆使してファイルを読み込
むことができる。
②

49

#以降は無視される③先頭に#がついているものは無視される(実行されない)。つまり②のコマンドは無効で、①のコマンドのみが実行される。①だけではこのファイルをどこから取得したのかわからないが、このようにコメントアウト(#をつけること)して完全なURL情報がわかるようにしている。このあたりはLinuxのシェルスクリプトと同じ。

③①

Jul 29 2015②50

list.files, file.info①getwd()で作業フォルダの確認。②list.files()で解析したいファイルの存在確認。”supp”を含むファイル名のもののみ出力させるテクニック。③file.info()で④ファイルサイズ(約4.5MB)などの詳細情報がわかる。

①
②
③
④

Jul 29 201551

(ほぼ)対応するLinuxコマンド

参考Linuxの場合
①
②
③
④

①
②
③
④

Jul 29 201552

読み込み確認①黒枠部分をコピペして、読み込めていることを確認。②コピー時に、③灰色部分は「反転」しないのでコピーできているか不安かもしれない
が、ちゃんとコピーできているので気にしない。
①②③

Jul 29 201553

右クリックで①ペースト

読み込み確認①

Jul 29 201554

読み込み確認

read.table関数を用いてsuppTable1.xlsを読み込
む際、①ヘッダー行あり(header=T)として、また②(行名として用いるため)1列目を行名(row.names=1)としている。このため、残りのデータは①20,689行×55列となる。

①

Jul 29 201555

suppTable1.xls

確かに入力ファイル(suppTable1.xls)は、①の幅的にも55列くらいありそうだと納得できる。
また、②2列目以降からすぐにカウントデータになっているわけではないこともわかる。

②

①

①

Jul 29 201556

行列の一部を抽出して表示。行列hogeの①1-7
行目、および②1-6列目を抽出して表示。こんな感じでうまく読み込めていることを確認する。

suppTable1.xls①②

Jul 29 201557

head関数①head関数を用いて、最初の1行分のみ表示。55列分もあるので、1行だけ表示させるのでも結構な画面サイズを要する。

①

Jul 29 201558

suppTable1.xls①別の表現方法。黒枠の列以降が目的のカウント情報であることが読み取れる。これはIlluminaのRNA-seqカウント
データ。Illuminaは実験単位をラン(Run; R)で表現する。また1つのRun中に複数のレーン(Lane; L)があるので複数サンプルを流せる。それゆえR1L1.HSM1は、Run1のLane1に流したHSM1というサンプルのカウントデータと読み取る。

①

Jul 29 201559

suppTable1.xls


ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



このデータは、3種類の生物種間比較。ヒト(Homosapiens; HS)、チンパンジー(Pan troglodytes; PT)、アカゲザル(Rhesus macaque; RM)。生物種ごとにオス3匹、メス3匹。雄雌を考慮しなければbiologicalreplicates (生物学的な反復)は6。①黒枠はヒトのオスで、個体識別番号が3のデータ(HSM3)と解釈する。

チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)①

Jul 29 201560

suppTable1.xls


ヒト(HS)

よく見ると、①Run3のLane6で流したHSM3 (i.e., R3L6.HSM3)以外にも、②Run4のLane1で流した同じHSM3のデータ(i.e.,
R4L1.HSM3)が存在する。これらは、同一個体由来データである。つまり、technical replicates (技術的な反復)は2である。

 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)①②

Jul 29 201561

suppTable1.xls


個体ごとにサンプルを分割して得たデータが全個体について存在する。technical replicates (技術的な反復)は2。

ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 201562

colnames関数


行列hogeの①列名を表示。目的のカウン
トデータが20列目以降にあることが分かる。

ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)①

 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015②63lengthは要素数


ヒト(HS)
 オス3匹(M1, M2, M3)①行列hogeの列名の20-55番目の要素のみを表示。(55– 20 +1)=36個の要素数と手計算できるが、length関数を用いて、②オリジナルが55個の要素、③サブセットの要素数が36個という結果を得ることもできる。lengthは、要素数分だけループを回したりする際にも用いられる。

 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)①②③
Jul 29 201564

列の並びがイマイチヒト(HS)
 オス3匹(M1, M2, M3)①行列hoge中の20-55番目の列を抽出した結果をdataに格納。これがsubsettingの基本形。②行列
dataの最初の1行目のみ表示。うまく抽出できていることがわかる。③しかしよく見ると、生物種ごとのようなきれいな並びになっていないので、イマイチ。

 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

①

アカゲザル(RM)②
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 201565

①全角はアリエマセン

嘘のようなホントの話

ヒト(HS)←

 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)①
Jul 29 201566

列名で並び替えheader=TRUEとしたおかげで、data$列名、hoge$列
名などとして特定の列のみ取り扱える。

ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 201567

任意のベクトル同士を列(column)方向で結合(bind)するのがcbind関数。①列を単純に結合することができる。

cbind関数

ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM)①

 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015

68

cbind関数


ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



こんな書き方もできる。①の拡張版として「HSM1, HSM2,
HSM3, HSF1, …, RMF2, RMF3」のような並びにしておけば、後の発現変動解析のときにいろいろと便利。②同一個体の反復データを足す場合。これはtechnical replicatesデータをマージ(合併)させることに相当する。一般的な発現変動解析は、
technical replicatesデータをマージして、biological replicatesのみからなるデータにしたものを入力として行う。

チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM)①

 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)②

Jul 29 201569

元データを整形

ここまでの説明で、例題41の下記コードの中身がかなり理解できるはずです①②

Jul 29 201570

元データを整形

Jul 29 2015

例題41をコピペ。Internet Explorerのヒト
はCTRLとALTキーを押しながらコードの枠内で左クリックすると全選択できます。

71

正常終了時の状態。①出力ファイル(sample_blekhman_36.txt)の中身は、ヘッ
ダー行や行名部分を除くと、②20,689行
×36列からなるカウントデータ行列。

元データを整形①②

Jul 29 201572

データ概観Jul 29 2015

最初の2行分を表示。ヒト(HS)、チンパンジー(PT)、アカゲザル(RM)で意図通りに並び替えできていることがわかる。

73

データ概観

全体的に、四角で囲ったtechnical replicates(同一個体の反復)間の類似度が、biological replicates(同一生物種の別個体)間の類似度よりも高そうであることがわかる。

ヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015

74

EXCELで概観

出力ファイル(sample_blekhman_36.txt)をEXCELで眺めるとこんな感じ。①はENSG00000000971という遺伝子領域上に2,262リードマップされたことを表す。②はENSG00000001460の遺伝子領域上に3リードマップされたことを表す。もしこの2つの配列長が同じなら、マップされたリード数が多い前者①の発現レベルが高いという理解でよい。

①②

Jul 29 2015

75

EXCELで概観

Jul 29 2015

サンプル(列)ごとにマップされた総リード数を計算した結果。サンプル間比較の場合には、この総リード数を揃えるのが基本戦略。

76

EXCELで概観

もし揃えずに、例えば①と②のサンプル間比較(発現変
動遺伝子(DEG)検出)を行うと、①のほうが②に比べて全体的に(1,801,009 / 1,346,515 = 1.34)倍高発現な状態であることを意味するので、①で高発現となるDEGが多く検出されるだろう。もちろんそれは間違い。

①

Jul 29 2015②

77

colSums関数は、列ごとの総リード数を調べるときに便利。③総リード数の最小と最大を調べる場合は、range関数を利用する。

colSumsとrange
①
②

③

①

Jul 29 2015②

78

apply, min, max①

Jul 29 2015②

行列演算といえばapply関数。行列dataを入力として、①列ごと(MARGIN=2)に、②sum関数を実行せよ、という意味。総リード数の最小と最大は、range関数でなくてもminとmax関数を用いて別々に計算してもよい。
様々な関数を紹介しているが、自分が使う際はどれか一つでよい。一度でも見ておけば、少しでも記憶に残るだろうという思想のもと、羅列的に紹介している。

79

colMeans, rowMeans①①列ごとにマップされたリード数の平均を算出、②colMeans関数も同じ機能。③行ごとにマップされたリード数の平均を算出、④rowMeans
関数も同じ機能。⑤行ごとにマップされたリード数を算出。rowSums関数は、低発現遺伝子のフィルタリング時にも利用される。

②

③

④

⑤

Jul 29 2015

80

EXCEL上での見た目とも一致してますね

EXCELと比較

Jul 29 2015

81

参考

summary関数

サンプルごとの要約統計量を概観する場合によく用いる。
ここでは、最初の6サンプル分(HS群のメス)に絞って表示。
私の最初の着眼点は黒枠のあたり。特に、①1st Qu. (第一四分位数)が全6サンプルで0であることから、20,689遺伝子中の少なくとも25%はゼロカウントであることがわかる。
①①

Jul 29 2015

82

参考

summary関数

次に見るのは②Medianの値。これは2nd Qu. (第二四分位数)と同じである。サンプル全体にわたって、ここを概観する。そして、低発現遺伝子のフィルタリングの際に、（ここでは最初の6サンプル分しか示していないが）マップされたリード数が5以下のものを除く処理を行うと、半分以上が落とされるだろう、などの見込みをつける。
②②

Jul 29 2015

83

参考

summary関数

ちなみに私は、③Mean (平均値)をほとんど見ません。
一応見ますが重要視していません。黒枠内の数値の関係(Mean > 3rd Qu.)から、ごく一部の異常に高発現(リード数の多い)の遺伝子の影響がカナリ大きそうだから、この種の外れ値の効果を排除できないMeanのような要約統計量は使わないほうがよいと判断します。
③③

Jul 29 201584

実用上は…


総リード数補正(RPM補正)

サンプル間比較の場合は、Rの発現変動解析
用Rパッケージをそのまま利用すればよい（うまくデータの正規化を行ってくれる）。8/5の統計
解析のところで発現変動解析を行う予定です。

 Mortazavi et al., Nat. Methods, 2008
 総リード数を100万など一定の値に揃えるベーシックな補正。外れ値に影響されやすい



TMM補正(edgeRパッケージ)
 Robinson and Oshlack, Genome Biol., 2010
 高発現側と低発現側で一定数をトリムして外れ値の影響を排除

TbT補正(TCCパッケージ)
 Kadota et al., Algorithms Mol. Biol., 2012
 TMMを含むedgeR (やDESeq)を内部的に利用して、高発現側と低発現側の外れ値に相当する発現変動遺伝子(DEG)をより正確に排除することで頑健な正規化を達成。
DEG-elimination strategy (DEGES)の基本形を提唱した論文


DEGES補正(TCCパッケージ)
 Sun et al., BMC Bioinformatics, 2013。TCC原著論文
 DEGESを一般化して、より高速かつ頑健な正規化を達成。edgeRやDESeq (後にDESeq2)の通常の手順を内部的に繰り返し実行して頑健な結果を得る枠組みを提供。
 Multi-group comparisonでもTCCの枠組みが有効であることを示した論文が近々…。
Jul 29 201585

クラスタリング

入力ファイルは20,689遺伝子×36サンプルのカウント
データファイル。ヒト(HS)、チンパンジー(PT)、アカゲザル(RM)の3生物種のデータ。各12サンプル。TCCパッケージを用いて、これのサンプル間クラスタリングを行います。

①②

Jul 29 201586

①出力は、hoge7.pngという名前のPNGファイル。②サイズは、700×400ピクセル。これは論文の図としても使えるレベル（実際我々の論文中でも使っている）

クラスタリング①②
hoge7.png

Jul 29 2015

ヒト(HS)

チンパンジー(PT)

アカゲザル(RM)

87

全個体について、同一個体を分割したtechnicalreplicatesのデータで末端のクラスターを形成していることが分かる。これはtechnical replicatesのデータ同士の類似度が非常に高いことを示している。妥当ですよね。

クラスタリングヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)

hoge7.png

 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015

ヒト(HS)

チンパンジー(PT)

アカゲザル(RM)88クラスタリングヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)

統計的手法で2群間比較(例えばMales vs. Females)をする目
的は、同一群内の別個体(biological replicates)のばらつきの程度を見積もっておき(モデル構築)、比較する2群間で発現に変動がないという前提(帰無仮説)からどれだけ離れているのかをp値で評価することである。ｐ値が低ければ低いほど「発現変動していない(帰無仮説に従う)」とは考えにくく、帰無仮説を棄却して「発現変動している(DEGである)」と判定することになる。

hoge7.png

 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015

ヒト(HS)

チンパンジー(PT)

アカゲザル(RM)89

サブセット抽出と整形①

統計的手法の多くは、biological replicatesのデータを前提としている。technical replicatesのデータをマージ(merge; collapseともいうらしい)したものを作成。③出力ファイルはsample_blekhman_18.txt。サンプル名部分は必要最小限の情報のみにしている。

②③

Jul 29 201590

クラスタリング

20,689遺伝子×18サンプルの
biological replicatesのみからなる
カウントデータでクラスタリング。

①②

Jul 29 201591

36サンプルのときの結果と同様、全体的なトポロジーは同じ。このクラスタリング結果を眺めるだけで、DEG検出
結果のイメージは大体つかめる。例1：「HS vs. RMで得られるDEG数」のほうが「HS vs. PTで得られるDEG数」よりも多そう。例2：ヒトは「オス vs. メス」でのDEG数は0に近いだろう。例3：RMM2が外れサンプルっぽいので、これを除去すれば生物種間比較時にDEG数が増えるだろう。
8/5の統計解析のところで発現変動解析を行う予定。

クラスタリングヒト(HS)
 オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)



チンパンジー(PT) オス3匹(M1, M2, M3)

hoge8.png

 メス3匹(F1, F2, F3)



アカゲザル(RM) オス3匹(M1, M2, M3)
 メス3匹(F1, F2, F3)

Jul 29 2015

ヒト(HS)

チンパンジー(PT)

アカゲザル(RM)

92

Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 201593

https://ja.wikipedia.org/wiki/FASTA Rでmulti-FASTAファイルを読み込んで自在に解析できます。ゲノム配列解析≒FASTA形式ファイルの解析。ここでは全体像を完全に把握すべくhoge4.fa
ファイルを仮想ゲノム配列ファイルとして取り扱う。

FASTA形式

Jul 29 201594

ゲノム配列Jul 29 2015

実際のゲノム配列はここからも取得可能。Rで染色体ごとの配列長やGC含量の計算ができる。

95

入力と出力の関係

multi-FASTAファイルを読み込んで、トータルの配列長、染色体数（コンティ
グ数）、配列長の平均、中央値、最大値、最小値、N50、GC含量を計算した結果を返すコードを実行してみよう。

入力: hoge4.fa
出力: hoge1.txt

Jul 29 201596

①ここです。コードの最初のほうに入力ファ
イルと出力ファイルを記述するので、コピペで実行した結果としてどういう名前のファイルが出力されるべきかわかる。hoge4.faはhogeフォルダ中にもありますが、②ここからも右クリックでダウンロードできます。

基本情報取得①②

Jul 29 201597

getwdとlist.files①例題の入力ファイル（hoge4.fa)をダウンロード。
②R上で作業ディレクトリの確認。③作業ディレク
トリに解析したい入力ファイルがあることを確認。

①②③

Jul 29 2015

98

コピペ
①

①一連のコマンド群をコピーして②R Console画面上でペースト。
ブラウザがInternet Explorerの場
合は、CTRLとALTキーを押しな
がらコードの枠内で左クリックすると、全選択できます。

②

Jul 29 201599

コピペ後に①list.files()で、②出力ファイル名として指定したhoge1.txtが作成されていることを確認。

実行結果
②

①

Jul 29 2015100

出力ファイルをテキストエディタやExcelで眺めてもよいが、オブジェクトtmpの中身を出力しているだけなので①R上で眺めている。

実行結果

①

Jul 29 2015

101

実行結果

contig_1の配列が最短、contig_2の配列が最長であることがわかる。入力と出力の関係を確認。

入力: hoge4.fa
出力: hoge1.txt

Jul 29 2015102

averageだと外れ値の影響を受けやすく、medianだと短いコンティグが多くを占める場合に不都合らしい。

N50


アセンブル結果の評価基準の一つ

出力: hoge1.txt

 長いコンティグから足していってTotal_lengthの50%に達したときのコンティグの長さ
 一般に数値が大きいほどよい

contig_2 (103 bp)
contig_3 (65 bp)

contig_4 (49 bp)
contig_1 (24 bp)

Total_length ×0.5 (120.5 bp)
Total_length (241 bp)

Jul 29 2015103

コード内部の説明

Jul 29 2015

コードの中身を説明します。
黒枠部分を再度コピペ。

104

コード内部の説明

Jul 29 2015

入力ファイル情報を格納したものがfastaオブジェクト。widthの位置にあるのがコンティグごとの配列長情報。

105

コード内部の説明

Jul 29 2015

width(fasta)にsum関数を適用すれば、トータルの配列長（配列長の総和）になる。

106

コード内部の説明

Jul 29 2015

length関数は要素数を返す。この場合、fastaオ
ブジェクトの要素数（つまりコンティグ数）を返す。

107

Tips：条件判定

Jul 29 2015

50 bp以上のコンティグからなるサブセットの抽出ができそうだ!

108

Tips：条件判定

Jul 29 2015

コードの中身が分かると応用範囲が飛躍的に増大。一定以上のスキルをもつバイオインフォマティシャンは「例題を探す」よりも「自分で作る」ヒトのほうが多いかも…。

109

Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 2015110

任意の領域の切り出し

subseq関数を用いて、任意の領域の配列を切り出すことができます。

入力：sample1.fasta>kadota
AGTGACGGTCTT

出力：hoge1.fasta>kadota
TGACGGTJul 29 2015111

Tips：関数のオプション

subseq関数実行時に、①数
値を直接指定してもいいし、②オプション名を明記してもよい。

①
入力：sample1.fasta>kadota
AGTGACGGTCTT

Jul 29 2015

出力：hoge1.fasta②>kadota
TGACGGT

112

Tips：関数のオプション①原因既知状態でエラーを出す。②「3番目の位置から5塩基分抽出」という他のオプション(endではなくwidth)を利用。

①②

入力：sample1.fasta>kadota
AGTGACGGTCTT

Jul 29 2015

出力：hoge1.fasta>kadota
TGACGGT

113

Tips：関数の使用法①①「?関数名」で使用法を記したウェブページ
が開く。ページの下のほうに、大抵の場合使用例が掲載されている。使用法既知の関数のマニュアルをいくつか読んで、慣れておこう。

…

Jul 29 2015

114

任意の領域の切り出し

入力がmulti-FASTAファイル（hoge4.fa）
で、リストファイル（list_sub2.txt）で指定した複数領域を切り出したい場合①②

Jul 29 2015115

こんな感じの結果が得られます

任意の領域の切り出し

Jul 29 2015

入力2: list_sub2.txt

116

妥当ですよね

任意の領域の切り出し
入力1: hoge4.fa

Jul 29 2015

入力2: list_sub2.txt

117FastQCと同じ結果を得る①100万リード、②107bpからなる③乳酸菌RNA-seqデータのFastQC解析結果のうち、例えば④のOverrepresented
sequencesと同じ結果をsubseqとtable関数を使って得ることができます。

③
①
②

④

Jul 29 2015

118

FastQCと同じ結果を得る①頻出する配列をリストアップ。②トッ
プは「CCCCGGTATA…」という50塩基の配列で14,383回出現。
Percentageは1.4383%。全部で100万リードなので妥当。オリジナル107 bpのうち最初の50 bpで解析している。

②①

Jul 29 2015

result_without_nogroup.html

119Overrepresented seq.

subseq関数を使っています。やってみましょう。

①②

Jul 29 2015120

Overrepresented seq.

Jul 29 2015

完璧に同じ結果を得られていることが分かります

121Contents

R基礎(初級) おさらい コード内部の説明（ファイルの読み込み、行列演算の基礎)
 リアルRNA-seqカウントデータ（数値行列データ）

R各種パッケージ(中級)：代表的なパッケージの利用法 (パッケージのインストール法)
 基本情報取得（コンティグ数、配列長、N50、GC含量）
 任意の領域の切り出し GC含量計算部分の説明

Jul 29 2015122

GC含量計算部分の説明

右のサイドバーを下に移動させるとGC含量計算部分を見られる。

①②

Jul 29 2015123

GC含量計算部分の説明

Jul 29 2015

fastaオブジェクトを出発点として、配列全体のGC含量（57.68%）を得るところの説明です。

124

GC含量計算部分の説明

黒枠部分を再度コピペしたのち、①fasta
オブジェクトの中身を表示させたところ。

①

Jul 29 2015

125

GC含量計算部分の説明

Jul 29 2015

alphabetFrequency関数は、塩基ごとの出現回数を返す。

126

GC含量計算部分の説明

DNA配列上のMは「A or C」、Rは「A orG」などというルールがあるようです。

http://en.wikipedia.org/wiki/FASTA_format
Jul 29 2015127

GC含量計算部分の説明①dim関数は行列の行数と列数を返す。alphabetFrequency関数出力結果は、4行×18列からなることが分かる。キーボードの上下キーを上手に利用して最小限の労力でキータイプ（あるいはコピペ）すべし!

①

Jul 29 2015128

GC含量計算部分の説明

Jul 29 2015

任意のサブセットを取得可能。
2:3やc(1,4)などをうまく利用。

129

GC含量計算部分の説明

Jul 29 2015

黒丸中の数値はcontig_1中のAの数が4個、赤丸中の数値は、contig_4中のTの数が10個であるということ。
rowSums関数は行ごとの和を返す。

130

GC含量計算部分の説明

Jul 29 2015

rowSums関数の入力として、
ACGTのみのカウント数を与えているが、その結果（返り値）は、配列
中にNなどを含まない場合は実質的にコンティグごとの配列長と同じ。

131

GC含量計算部分の説明

Jul 29 2015

オブジェクトCG中には、配列(コンティグ)
ごとのCとGのカウント数が格納されている。オブジェクトACGT中には、配列ごとのA, C, G, Tのカウント数が格納されている。例えば49塩基からなるcontig_4中に、ACGTの4種類の塩基が49個、CGの数は25個あることを意味する。sum関
数は、ベクトルの要素の和を返す。

132

GC含量計算部分の説明

ここでは①sum関数を用いて配列全体の総和でGC含量計算をしているが、②sum関数を用いずに「CG/ACGT」とやると、コンティグごとのGC含量を得られる。例えば、
contig_1はCGの数が16個で、ACGTの数が24個。それゆえGC
含量は16/24 = 0.6666667となる。

①②

Jul 29 2015133

配列ごとのGC含量計算

Jul 29 2015

sum関数を用いずに「CG/ACGT」とやって、コンティグごとのGC含量を得るための項目。記述内容がほぼ同じであることが分かる。

134配列ごとのGC含量計算

Jul 29 2015

出力ファイル(hoge1.txt)中の一番右側の列が配列ごとのGC含量です。

135

配列ごとのGC含量計算

Jul 29 2015

ACGT列は4種類の塩基のみの出現数、Length列は配列長情報を表す。配列長は、ACGT以外の全てを含むので、2つの数値の差分(Length - ACGT)がNなどのACGT以外の塩基のトータルの出現回数ということになる。Length ≧ ACGTという関係。

136