205.07.27 版 配布する USB メモリ中の hoge フォルダをデスクトップにコピーしておいてください NGS ハンズオン 講習会 :R 基礎 東京大学 大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp http://www.iu.a.u-tokyo.ac.jp/~kadota/
Contents( 全体 ) 7 月 22 日 ( 水 ):84 83 名 Bio-Linux 8 と R のインストール状況確認 基本自習 ( 門田 寺田先生 ) 7 月 23 日 ( 木 ):92 90 名 Linux 基礎 Linux コマンドなど UNIX の基礎の理解 ( 門田 ) 7 月 24 日 ( 金 ):85 83 名 スクリプト言語 シェルスクリプト ( アメリエフ株式会社服部恵美先生 ) 7 月 27 日 ( 月 ):93 9 名 スクリプト言語 Perl( アメリエフ服部先生 ) 7 月 28 日 ( 火 ):9 90 名 スクリプト言語 Python( アメリエフ服部先生 ) 7 月 29 日 ( 水 ):94 88 名 データ解析環境 R( 門田 ) R 基礎 ( 初級 ):R 言語の基礎 ( インストールから利用まで ) R 基礎 2( 初級 ): ファイルの読み込み 行列演算の基本 R 各種パッケージ ( 中級 ): パッケージのインストール法と代表的なパッケージの利用法 7 月 30 日 ( 木 ):96 9 名 データ解析環境 R( 門田 ) Bioconductor の利用法 ( 中級 ): データの型やバージョンの違い Bioconductor の利用法 2( 中級 ):FASTA/FASTQ ファイルの各種解析 8 月 3 日 ( 月 ):89 84 名 NGS 解析 基礎 ( アメリエフ山口昌雄先生 ) 8 月 4 日 ( 火 ):85 80 名 NGS 解析 ゲノム Reseq 変異解析 ( アメリエフ山口先生 ) 8 月 5 日 ( 水 ):86 8 名 NGS 解析 RNA-seq 統計解析 ( 前半 : 山口先生 後半 : 門田 ) 8 月 6 日 ( 木 ):04 98 名 NGS 解析 ChIP-seq( 理研森岡勝樹先生 ) 2
各種ソフトの場所 2 4 2Excel は行列データファイルの確認用 門田は EmEditor というテキストエディタを使っています 3 受講生の心構え でも書いていますが 貸与 PC のほとんどは R ver. 3..2, 3..3, 3.2.0, 3.2. のいずれか ( または複数 ) がインストールされています 基本的には最新版を利用 4 エディタは R 付属のものを推奨 主目的は二重クォーテーション問題の回避 3
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 4
おさらい hoge フォルダ中の r_seq.html をダブルクリックしてローカルで利用するのがいいかもしれません ここで示すようなクリックして眺めるだけのネットサーフィン系の部分は 手を動かさずに前のスライドを見ているだけのほうがいいかもしれません 5
おさらい 基本的な利用法 PDF 中の 解析基礎 2 をおさらいします 6
解析基礎 2 目的 : アノテーションファイル (annotation.txt) 中の第 列目に対して リストファイル (genelist.txt) 中の文字列と一致する行を抜き出して hoge.txt というファイル名で出力したい 入力 : アノテーションファイル (annotation.txt) 出力 :hoge.txt 入力 2: リストファイル (genelist.txt) 7
解析基礎 2 目的 : アノテーションファイル (annotation.txt) 中の第 列目に対して リストファイル (genelist.txt) 中の文字列と一致する行を抜き出して hoge.txt というファイル名で出力したい 8
解析基礎 2 作業ディレクトリは デスクトップ hoge hoge フォルダ中に annotation.txt と genelist.txt が存在するという前提 貸与 PC は黒矢印部分が kadota ではなく iu 9
基本はコピペ 一連のコマンド群をコピーして作業ディレクトリは デスクトッ 2R Console 画面上でペースト ブラウザがプ hoge Internet hogeフォルダ中 Explorerの場合は CTRLとALT にannotation.txt キーを押しながらコードのとgenelist.txt 枠内で左クリックすると 全選択できます が存在するという前提 2 0
実行結果 出力ファイル名として指定した hoge.txt が生成されているのがわかる list.files() で表示される結果 と 実行後の hoge フォルダの中身 は当然同じです 実行前の hoge フォルダ 実行後の hoge フォルダ
実行結果 out というオブジェクトの中身を write.table という関数でファイルに出力しています この場合 出力ファイル hoge.txt の中身は R コンソール画面中で out と打ち込むことで見られる 実行後の hoge フォルダ 2
色の説明 R コード中の色の使い分けについて説明します 3
応用 このサンプルコードは 列目でキーワード検索する場合 別のリストファイルを読み込んで 4 列目で検索したい場合のやり方を示します 4
解答例. 目的のキーワードリストを含むファイルを作成し ( 例 :list.txt) 2. 該当箇所を変更し Rコンソール画面上でコピペ メモ帳 など任意のエディタでリストファイル (list.txt) を作成 5
解答例. 目的のキーワードリストを含むファイルを作成し ( 例 :list.txt) 2. 該当箇所を変更し Rコンソール画面上でコピペ 一連の作業手順を記述したスクリプトを つのファイルとして保存することをお勧め 6
ありがちなミス 作業ディレクトリの変更を忘れているため in_f で指定した最初のファイルの読み込み段階でエラーが出る つまり 実際に行ったフォルダ中には annotation.txt というファイルは存在しないということ 7
ありがちなミス 2 必要な入力ファイルが作業ディレクトリ中に存在しない この場合 in_f2 で指定した genelist.txt が存在しないため それの読み込み段階でエラーが出ている それゆえ その情報を用いているコマンド部分でエラーが出ている 8
ありがちなミス 3 出力予定のファイル名と同じものをエクセルなど別のプログラムで開いているため 最後の write.table 関数のところでエラーが出る 対処法は 出力ファイル名を変更するか 開いている別のプログラムを閉じる 9
ありがちなミス 4 実行スクリプトをコピーする際 最後の行のところで改行を含ませずに R Console 画面上でペーストしたため 最後のコマンドが実行されない ( 出力ファイルが生成されない ) これも比較的ありがちなパターンです コピペ後に無意識にリターンキーを押すことを心がけるだけでもよいでしょう 20
警告メッセージ list.txt ファイル作成時に membrane と打った後に改行を 入れた場合と 2 入れない場合の挙動の違いを把握し 後学のために警告メッセージの意味を理解しておくとよい この場合は結果には影響していないことがわかる R は警告メッセージの記述内容が比較的分かりやすいのでよく読むべし 2 2
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 22
コード内部の説明 コードの中身を説明します 黒枠部分を再度コピペ 23
読み込み in_f で指定したファイルを読み込め 2 読み込むファイルの最初の行はヘッダー部分 3 ファイルの区切り文字はタブです 4 読み込んだ結果を data という名前で取り扱う 4 2 3 24
行列 data 2 data と打ってリターン 入力ファイルの中身を正しく読み込めていることがわかる 2header=TRUE としているので 3 このように見えて列名として認識される 3 25
dim で行数と列数を表示 オブジェクト data の行数と列数は と 4 2 ウェブページ中の表記が灰色なのは 特にやらなくてもいいコマンドだから 2 26
行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] humei は 読み込み元ファイルの annotation.txt 中では 7 行 4 列目だが 2 行目をヘッダー行としているので 3 6 行 4 列目とする必要がある 利用例は ファイル読み込み時に x 行 y 列目に不具合がある のようなエラーが出た時のトラブルシューティングなど 2 3 27
Tips: 上下左右の矢印キー 上矢印キーを押すと 直前に打ったコマンドが表示される 最初から全部打ち直すのではなく 上下左右の矢印キーを有効に利用し最小限の労力で打つべし! 28
行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] 2 行目の情報のみ抽出 読み込み時に head=true としていたので ヘッダー行がついていることが分かる 2 29
行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] 2 列目の情報のみ抽出 30
行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] param 列目の情報のみ抽出 2param には という数値が代入されていたのでこうなる 2 2 3
Tips: 関数とオプション 参考 行列 data の最初の数行を表示したい場合は head 関数を利用 n=3 というオプションを利用すると最初の 3 行分のみ表示 関数ごとに様々なオプションを利用可能です このあたりは 2Linux とよく似ている 2 32
Tips: タブ補完 参考 列番号を指定する以外にも特定の列を表示するやり方がある head=true で入力ファイルを読み込むと 列の名前を利用することができる subcellular_location 列の情報を抽出したい場合は 2 data$su くらいまで打ち込んでから Tab キーを押す 2 33
Tips: タブ補完 参考 列名中の su からはじまる文字列を補完して表示してくれる Tab キーを用いた補完機能 という意味で タブ補完 という このテクニックは Linux でも利用可能 2 34
Tips:table 関数 参考 table 関数は ベクトル中の要素ごとの出現回数を返す NGS データ中の特定のリードの出現回数 ( 後述 ) や アノテーションファイル中の染色体ごとの遺伝子数 など 様々な局面で利用可能 35
Tips: ソート 参考 sort 関数と併用することで全体像を俯瞰可能 例えば nuclear に局在する遺伝子数が最も多く 4 個であった などが簡単にわかる 36
Tips:is.element 関数 hoge ベクトルに対して nuclear の文字が存在する場所を TRUE 存在しない場所を FALSE として返す as.character 関数は 文字列ベクトルとして取り扱いたい場合に利用 37
Tips: 二重クォーテーション 二重クォーテーションが自動で変更されるエディタは非推奨です 日本語の二重クォーテーションもだめです Microsoft Word や PDF ファイル中のコードのコピペ時によくハマります 38
目的をおさらい 目的は 数万 ~ 数百万行からなるファイルを読み込んで特定のキーワードを含む行のみ取り出すテクニックを習得 39
目的をおさらい 論理値ベクトル obj を用いて TRUE の要素に対応する行を抽出している 入力 2: リストファイル (genelist.txt) 40
目的をおさらい コード作成当時は as.character 関数を用いてデータの型を文字列ベクトルに揃えていた 少なくとも現在 (R ver. 3..3 以降 ) は この関数がなくても大丈夫なようだ 同じ関数でもバージョンによって挙動が異なるということ ( バージョンの違いの一例 ) 4
と 2 は手順が異なるだけで実質的に同じです genelist.txt 42
このコードはヘッダー行がある場合のものです 入力 : annotation.txt 出力 : hoge2.txt 43
このコードはヘッダー行がない場合のものです 入力 : annotation2.txt 出力 : hoge3.txt 44
ヘッダー行がある場合 ヘッダー行がない場合 45
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 46
カウントデータ 目的サンプルの RNA-Seq データ mapping リファレンス配列 : ゲノム 教科書 p8-89 教科書 カウントデータとは マップされたリード数 をカウントしたデータのこと 以下の例では サンプルなので 列分のデータしかないが 一般には複数サンプルのデータを取得し サンプル間比較が行われるので複数の列からなる それゆえ 数値ベクトルではなく数値行列 詳細は 8/5 の RNA-seq 前半で count 遺伝子 遺伝子 2 遺伝子 3 遺伝子 4 47
数値行列 実験の詳細には立ち入らないが 3 生物種間比較を行った公共 RNA-seq カウントデータ (Blekhman et al., Genome Res., 200) を用いて R の王道的な使い方である数値行列解析のテクニックを伝授 8/5 の統計解析のところでこのデータを利用予定です 2 48
xls 形式ファイルも OK xls 形式のエクセルファイルを読み込むことができる ( 但し このファイルは壊れているなどというメッセージが出ており 実はタブ区切りテキストファイルなのに.xls という拡張子が無理やりつけられているというオチかもしれない )2 それほど大きなサイズでなければ ネットワーク経由で直接読み込むこともできる 他に read.csv や readlines 関数などを駆使してファイルを読み込むことができる 2 49
# 以降は無視される 3 先頭に # がついているものは無視される ( 実行されない ) つまり 2 のコマンドは無効で のコマンドのみが実行される だけではこのファイルをどこから取得したのかわからないが このようにコメントアウト (# をつけること ) して完全な URL 情報がわかるようにしている このあたりは Linux のシェルスクリプトと同じ 3 2 50
list.files, file.info getwd() で作業フォルダの確認 2list.files() で解析したいファイルの存在確認 supp を含むファイル名のもののみ出力させるテクニック 3file.info() で 4 ファイルサイズ ( 約 4.5MB) などの詳細情報がわかる 2 3 4 5
Linux の場合 2 3 参考 ( ほぼ ) 対応する Linux コマンド 4 2 3 4 52
読み込み確認 黒枠部分をコピペして 読み込めていることを確認 2 コピー時に 3 灰色部分は 反転 しないのでコピーできているか不安かもしれないが ちゃんとコピーできているので気にしない 2 3 53
読み込み確認 右クリックで ペースト 54
読み込み確認 read.table 関数を用いて supptable.xls を読み込む際 ヘッダー行あり (header=t) として また 2( 行名として用いるため ) 列目を行名 (row.names=) としている このため 残りのデータは 20,689 行 55 列となる 55
supptable.xls 確かに入力ファイル (supptable.xls) は の幅的にも 55 列くらいありそうだと納得できる また 22 列目以降からすぐにカウントデータになっているわけではないこともわかる 2 56
supptable.xls 行列の一部を抽出して表示 行列 hoge の -7 行目 および 2-6 列目を抽出して表示 こんな感じでうまく読み込めていることを確認する 2 57
head 関数 head 関数を用いて 最初の 行分のみ表示 55 列分もあるので 行だけ表示させるのでも結構な画面サイズを要する 58
supptable.xls 別の表現方法 黒枠の列以降が目的のカウント情報であることが読み取れる これは Illumina の RNA-seq カウントデータ Illumina は実験単位をラン (Run; R) で表現する また つの Run 中に複数のレーン (Lane; L) があるので複数サンプルを流せる それゆえ RL.HSM は Run の Lane に流した HSM というサンプルのカウントデータと読み取る 59
supptable.xls ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) このデータは 3 種類の生物種間比較 ヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) 生物種ごとにオス 3 匹 メス 3 匹 雄雌を考慮しなければ biological replicates ( 生物学的な反復 ) は 6 黒枠はヒトのオスで 個体識別番号が 3 のデータ (HSM3) と解釈する 60
supptable.xls ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) よく見ると Run3 の Lane6 で流した HSM3 (i.e., R3L6.HSM3) 以外にも 2Run4 の Lane で流した同じ HSM3 のデータ (i.e., R4L.HSM3) が存在する これらは 同一個体由来データである つまり technical replicates ( 技術的な反復 ) は 2 である 2 6
supptable.xls ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) 個体ごとにサンプルを分割して得たデータが全個体について存在する technical replicates ( 技術的な反復 ) は 2 62
colnames 関数 ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) 行列 hoge の 列名を表示 目的のカウントデータが 20 列目以降にあることが分かる 2 63
length は要素数 ヒト (HS) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) チンパンジー (PT) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) アカゲザル (RM) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) 行列 hoge の列名の 20-55 番目の要素のみを表示 (55 20 +)=36 個の要素数と手計算できるが length 関数を用いて 2 オリジナルが 55 個の要素 3 サブセットの要素数が 36 個という結果を得ることもできる length は 要素数分だけループを回したりする際にも用いられる 2 3 64
列の並びがイマイチ ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) 2 オス3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) 行列 hoge 中の 20-55 番目の列を抽出した結果を data に格納 これが subsetting の基本形 2 行列 data の最初の 行目のみ表示 うまく抽出できていることがわかる 3 しかしよく見ると 生物種ごとのようなきれいな並びになっていないので イマイチ 65
嘘のようなホントの話 ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) 全角はアリエマセン 66
列名で並び替え ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) header=true としたおかげで data$ 列名 hoge$ 列名などとして特定の列のみ取り扱える 67
cbind 関数 ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) 任意のベクトル同士を列 (column) 方向で結合 (bind) するのが cbind 関数 列を単純に結合することができる 68
cbind 関数 ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) こんな書き方もできる の拡張版として HSM, HSM2, HSM3, HSF,, RMF2, RMF3 のような並びにしておけば 後の発現変動解析のときにいろいろと便利 2 同一個体の反復データを足す場合 これは technical replicates データをマージ ( 合併 ) させることに相当する 一般的な発現変動解析は technical replicates データをマージして biological replicates のみからなるデータにしたものを入力として行う 2 69
元データを整形 ここまでの説明で 例題 4 の下記コードの中身がかなり理解できるはずです 2 70
元データを整形 例題 4 をコピペ Internet Explorer のヒトは CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できます 7
元データを整形 正常終了時の状態 出力ファイル (sample_blekhman_36.txt) の中身は ヘッダー行や行名部分を除くと 220,689 行 36 列からなるカウントデータ行列 2 72
データ概観 最初の 2 行分を表示 ヒト (HS) チンパンジー (PT) アカゲザル (RM) で意図通りに並び替えできていることがわかる 73
データ概観 ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) 全体的に 四角で囲った technical replicates( 同一個体の反復 ) 間の類似度が biological replicates( 同一生物種の別個体 ) 間の類似度よりも高そうであることがわかる 74
EXCEL で概観 出力ファイル (sample_blekhman_36.txt) を EXCEL で眺めるとこんな感じ は ENSG0000000097 という遺伝子領域上に 2,262 リードマップされたことを表す 2 は ENSG0000000460 の遺伝子領域上に 3 リードマップされたことを表す もしこの 2 つの配列長が同じなら マップされたリード数が多い前者 の発現レベルが高いという理解でよい 2 75
EXCEL で概観 サンプル ( 列 ) ごとにマップされた総リード数を計算した結果 サンプル間比較の場合には この総リード数を揃えるのが基本戦略 76
EXCEL で概観 もし揃えずに 例えば と 2 のサンプル間比較 ( 発現変動遺伝子 (DEG) 検出 ) を行うと のほうが 2 に比べて全体的に (,80,009 /,346,55 =.34) 倍高発現な状態であることを意味するので で高発現となる DEG が多く検出されるだろう もちろんそれは間違い 2 77
colsums と range colsums 関数は 列ごとの総リード数を調べるときに便利 3 総リード数の最小と最大を調べる場合は range 関数を利用する 2 3 2 78
apply, min, max 2 行列演算といえば apply 関数 行列 data を入力として 列ごと (MARGIN=2) に 2sum 関数を実行せよ という意味 総リード数の最小と最大は range 関数でなくても min と max 関数を用いて別々に計算してもよい 様々な関数を紹介しているが 自分が使う際はどれか一つでよい 一度でも見ておけば 少しでも記憶に残るだろうという思想のもと 羅列的に紹介している 79
colmeans, rowmeans 2 列ごとにマップされたリード数の平均を算出 2colMeans 関数も同じ機能 3 行ごとにマップされたリード数の平均を算出 4rowMeans 関数も同じ機能 5 行ごとにマップされたリード数を算出 rowsums 関数は 低発現遺伝子のフィルタリング時にも利用される 3 4 5 80
EXCEL と比較 EXCEL 上での見た目とも一致してますね 8
summary 関数 参考 サンプルごとの要約統計量を概観する場合によく用いる ここでは 最初の 6 サンプル分 (HS 群のメス ) に絞って表示 私の最初の着眼点は黒枠のあたり 特に st Qu. ( 第一四分位数 ) が全 6 サンプルで 0 であることから 20,689 遺伝子中の少なくとも 25% はゼロカウントであることがわかる 82
summary 関数 参考 次に見るのは 2Median の値 これは 2nd Qu. ( 第二四分位数 ) と同じである サンプル全体にわたって ここを概観する そして 低発現遺伝子のフィルタリングの際に ( ここでは最初の 6 サンプル分しか示していないが ) マップされたリード数が 5 以下のものを除く処理を行うと 半分以上が落とされるだろう などの見込みをつける 2 2 83
summary 関数 参考 ちなみに私は 3Mean ( 平均値 ) をほとんど見ません 一応見ますが重要視していません 黒枠内の数値の関係 (Mean > 3rd Qu.) から ごく一部の異常に高発現 ( リード数の多い ) の遺伝子の影響がカナリ大きそうだから この種の外れ値の効果を排除できない Mean のような要約統計量は使わないほうがよいと判断します 3 3 84
実用上は 総リード数補正 (RPM 補正 ) Mortazavi et al., Nat. Methods, 2008 総リード数を 00 万など一定の値に揃えるベーシックな補正 外れ値に影響されやすい TMM 補正 (edger パッケージ ) Robinson and Oshlack, Genome Biol., 200 高発現側と低発現側で一定数をトリムして外れ値の影響を排除 TbT 補正 (TCC パッケージ ) Kadota et al., Algorithms Mol. Biol., 202 TMM を含む edger ( や DESeq) を内部的に利用して 高発現側と低発現側の外れ値に相当する発現変動遺伝子 (DEG) をより正確に排除することで頑健な正規化を達成 DEG-elimination strategy (DEGES) の基本形を提唱した論文 DEGES 補正 (TCC パッケージ ) Sun et al., BMC Bioinformatics, 203 TCC 原著論文 サンプル間比較の場合は R の発現変動解析用 R パッケージをそのまま利用すればよい ( うまくデータの正規化を行ってくれる ) 8/5 の統計解析のところで発現変動解析を行う予定です DEGES を一般化して より高速かつ頑健な正規化を達成 edger や DESeq ( 後に DESeq2) の通常の手順を内部的に繰り返し実行して頑健な結果を得る枠組みを提供 Multi-group comparison でも TCC の枠組みが有効であることを示した論文が近々 85
クラスタリング 入力ファイルは 20,689 遺伝子 36 サンプルのカウントデータファイル ヒト (HS) チンパンジー (PT) アカゲザル (RM) の 3 生物種のデータ 各 2 サンプル TCC パッケージを用いて これのサンプル間クラスタリングを行います 2 86
クラスタリング 出力は hoge7.png という名前の PNG ファイル 2 サイズは 700 400 ピクセル これは論文の図としても使えるレベル ( 実際我々の論文中でも使っている ) 2 hoge7.png ヒト (HS) チンパンジー (PT) アカゲザル (RM) 87
クラスタリング ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) hoge7.png 全個体について 同一個体を分割した technical replicates のデータで末端のクラスターを形成していることが分かる これは technical replicates のデータ同士の類似度が非常に高いことを示している 妥当ですよね ヒト (HS) チンパンジー (PT) アカゲザル (RM) 88
クラスタリング ヒト (HS) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) チンパンジー (PT) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) アカゲザル (RM) オス 3 匹 (M, M2, M3) メス 3 匹 (F, F2, F3) hoge7.png 統計的手法で2 群間比較 ( 例えばMales vs. Females) をする目的は 同一群内の別個体 (biological replicates) のばらつきの程度を見積もっておき ( モデル構築 ) 比較する2 群間で発現に変動がないという前提 ( 帰無仮説 ) からどれだけ離れているのかをp 値で評価することである p 値が低ければ低いほど 発現変動していない( 帰無仮説に従う ) とは考えにくく 帰無仮説を棄却して 発現変動している (DEGである) と判定することになる ヒト (HS) チンパンジー (PT) アカゲザル (RM) 89
サブセット抽出と整形 統計的手法の多くは biological replicates のデータを前提としている technical replicates のデータをマージ (merge; collapse ともいうらしい ) したものを作成 3 出力ファイルは sample_blekhman_8.txt サンプル名部分は必要最小限の情報のみにしている 2 3 90
クラスタリング 20,689 遺伝子 8 サンプルの biological replicates のみからなるカウントデータでクラスタリング 2 9
クラスタリング ヒト (HS) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) チンパンジー (PT) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) アカゲザル (RM) オス3 匹 (M, M2, M3) メス3 匹 (F, F2, F3) hoge8.png 36 サンプルのときの結果と同様 全体的なトポロジーは同じ このクラスタリング結果を眺めるだけで DEG 検出結果のイメージは大体つかめる 例 : HS vs. RM で得られる DEG 数 のほうが HS vs. PT で得られる DEG 数 よりも多そう 例 2: ヒトは オス vs. メス での DEG 数は 0 に近いだろう 例 3:RMM2 が外れサンプルっぽいので これを除去すれば生物種間比較時に DEG 数が増えるだろう 8/5 の統計解析のところで発現変動解析を行う予定 ヒト (HS) チンパンジー (PT) アカゲザル (RM) 92
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 93
https://ja.wikipedia.org/wiki/fasta FASTA 形式 R で multi-fasta ファイルを読み込んで自在に解析できます ゲノム配列解析 FASTA 形式ファイルの解析 ここでは全体像を完全に把握すべく hoge4.fa ファイルを仮想ゲノム配列ファイルとして取り扱う 94
ゲノム配列 実際のゲノム配列はここからも取得可能 R で染色体ごとの配列長や GC 含量の計算ができる 95
入力と出力の関係 入力 : hoge4.fa multi-fasta ファイルを読み込んで トータルの配列長 染色体数 ( コンティグ数 ) 配列長の平均 中央値 最大値 最小値 N50 GC 含量を計算した結果を返すコードを実行してみよう 出力 : hoge.txt 96
基本情報取得 ここです コードの最初のほうに入力ファイルと出力ファイルを記述するので コピペで実行した結果としてどういう名前のファイルが出力されるべきかわかる hoge4.fa は hoge フォルダ中にもありますが 2 ここからも右クリックでダウンロードできます 2 97
getwd と list.files 例題の入力ファイル (hoge4.fa) をダウンロード 2R 上で作業ディレクトリの確認 3 作業ディレクトリに解析したい入力ファイルがあることを確認 2 3 98
コピペ 一連のコマンド群をコピーして 2R Console 画面上でペースト ブラウザが Internet Explorer の場合は CTRL と ALT キーを押しながらコードの枠内で左クリックすると 全選択できます 2 99
実行結果 コピペ後に list.files() で 2 出力ファイル名として指定した hoge.txt が作成されていることを確認 2 00
実行結果 出力ファイルをテキストエディタや Excel で眺めてもよいが オブジェクト tmp の中身を出力しているだけなので R 上で眺めている 0
実行結果 contig_ の配列が最短 contig_2 の配列が最長であることがわかる 入力と出力の関係を確認 入力 : hoge4.fa 出力 : hoge.txt 02
N50 アセンブル結果の評価基準の一つ 長いコンティグから足していって Total_length の 50% に達したときのコンティグの長さ 一般に数値が大きいほどよい average だと外れ値の影響を受けやすく median だと短いコンティグが多くを占める場合に不都合らしい 出力 : hoge.txt contig_2 (03 bp) Total_length 0.5 (20.5 bp) contig_3 (65 bp) contig_4 (49 bp) contig_ (24 bp) Total_length (24 bp) 03
コード内部の説明 コードの中身を説明します 黒枠部分を再度コピペ 04
コード内部の説明 入力ファイル情報を格納したものが fasta オブジェクト width の位置にあるのがコンティグごとの配列長情報 05
コード内部の説明 width(fasta) に sum 関数を適用すれば トータルの配列長 ( 配列長の総和 ) になる 06
コード内部の説明 length 関数は要素数を返す この場合 fasta オブジェクトの要素数 ( つまりコンティグ数 ) を返す 07
Tips: 条件判定 50 bp 以上のコンティグからなるサブセットの抽出ができそうだ! 08
Tips: 条件判定 コードの中身が分かると応用範囲が飛躍的に増大 一定以上のスキルをもつバイオインフォマティシャンは 例題を探す よりも 自分で作る ヒトのほうが多いかも 09
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 0
任意の領域の切り出し subseq 関数を用いて 任意の領域の配列を切り出すことができます 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT
Tips: 関数のオプション subseq 関数実行時に 数値を直接指定してもいいし 2 オプション名を明記してもよい 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta 2 >kadota TGACGGT 2
Tips: 関数のオプション 原因既知状態でエラーを出す 2 3 番目の位置から 5 塩基分抽出 という他のオプション (end ではなく width) を利用 2 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT 3
Tips: 関数の使用法? 関数名 で使用法を記したウェブページが開く ページの下のほうに 大抵の場合使用例が掲載されている 使用法既知の関数のマニュアルをいくつか読んで 慣れておこう 4
任意の領域の切り出し 入力が multi-fasta ファイル (hoge4.fa) で リストファイル (list_sub2.txt) で指定した複数領域を切り出したい場合 2 5
任意の領域の切り出し こんな感じの結果が得られます 入力 2: list_sub2.txt 6
任意の領域の切り出し 入力 : hoge4.fa 入力 2: list_sub2.txt 妥当ですよね 7
FastQC と同じ結果を得る 00 万リード 207bp からなる 3 乳酸菌 RNA-seq データの FastQC 解析結果のうち 例えば 4 の Overrepresented sequences と同じ結果を subseq と table 関数を使って得ることができます 3 2 4 8
FastQC と同じ結果を得る 頻出する配列をリストアップ 2 トップは CCCCGGTATA という 50 塩基の配列で 4,383 回出現 Percentage は.4383% 全部で 00 万リードなので妥当 オリジナル 07 bp のうち最初の 50 bp で解析している 2 result_without_nogroup.html 9
Overrepresented seq. subseq 関数を使っています やってみましょう 2 20
Overrepresented seq. 完璧に同じ結果を得られていることが分かります 2
Contents R 基礎 ( 初級 ) おさらい コード内部の説明 ( ファイルの読み込み 行列演算の基礎 ) リアル RNA-seq カウントデータ ( 数値行列データ ) R 各種パッケージ ( 中級 ): 代表的なパッケージの利用法 ( パッケージのインストール法 ) 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 22
GC 含量計算部分の説明 右のサイドバーを下に移動させると GC 含量計算部分を見られる 2 23
GC 含量計算部分の説明 fasta オブジェクトを出発点として 配列全体の GC 含量 (57.68%) を得るところの説明です 24
GC 含量計算部分の説明 黒枠部分を再度コピペしたのち fasta オブジェクトの中身を表示させたところ 25
GC 含量計算部分の説明 alphabetfrequency 関数は 塩基ごとの出現回数を返す 26
GC 含量計算部分の説明 DNA 配列上の M は A or C R は A or G などというルールがあるようです http://en.wikipedia.org/wiki/fasta_format 27
GC 含量計算部分の説明 dim 関数は行列の行数と列数を返す alphabetfrequency 関数出力結果は 4 行 8 列からなることが分かる キーボードの上下キーを上手に利用して最小限の労力でキータイプ ( あるいはコピペ ) すべし! 28
GC 含量計算部分の説明 任意のサブセットを取得可能 2:3 や c(,4) などをうまく利用 29
GC 含量計算部分の説明 黒丸中の数値は contig_ 中の A の数が 4 個 赤丸中の数値は contig_4 中の T の数が 0 個であるということ rowsums 関数は行ごとの和を返す 30
GC 含量計算部分の説明 rowsums 関数の入力として ACGT のみのカウント数を与えているが その結果 ( 返り値 ) は 配列中に N などを含まない場合は実質的にコンティグごとの配列長と同じ 3
GC 含量計算部分の説明 オブジェクト CG 中には 配列 ( コンティグ ) ごとの C と G のカウント数が格納されている オブジェクト ACGT 中には 配列ごとの A, C, G, T のカウント数が格納されている 例えば 49 塩基からなる contig_4 中に ACGT の 4 種類の塩基が 49 個 CG の数は 25 個あることを意味する sum 関数は ベクトルの要素の和を返す 32
GC 含量計算部分の説明 ここでは sum 関数を用いて配列全体の総和で GC 含量計算をしているが 2sum 関数を用いずに CG/ACGT とやると コンティグごとの GC 含量を得られる 例えば contig_ は CG の数が 6 個で ACGT の数が 24 個 それゆえ GC 含量は 6/24 = 0.6666667 となる 2 33
配列ごとの GC 含量計算 sum 関数を用いずに CG/ACGT とやって コンティグごとの GC 含量を得るための項目 記述内容がほぼ同じであることが分かる 34
配列ごとの GC 含量計算 出力ファイル (hoge.txt) 中の一番右側の列が配列ごとの GC 含量です 35
配列ごとの GC 含量計算 ACGT 列は 4 種類の塩基のみの出現数 Length 列は配列長情報を表す 配列長は ACGT 以外の全てを含むので 2 つの数値の差分 (Length - ACGT) が N などの ACGT 以外の塩基のトータルの出現回数ということになる Length ACGT という関係 36