case1 - nibb-unix/gitc202502-unix GitHub Wiki
実践演習1 RNA-Seq 解析結果の集計(1)
（以下の実践演習はローカルPC上で行うこと）
UNIXコマンドプロンプトで gitc/data/9_ex/rnaseq に移動せよ。
RNA-Seq 解析のデータを使って、Unix やR の基本コマンドを使ってもう一歩先の解析を行ってみよう。
この下の results ディレクトリには、12サンプルのRNA-Seq データをそれぞれ大腸菌ゲノムにマッピングした結果から、
htseq-count コマンド(https://htseq.readthedocs.io/en/release_0.11.1/count.html)を使って、
遺伝子ごとにその領域に重なるリード数を単純にカウントした結果(ecoli.*.htseq)が格納されている。UNIX コマンドを使って、この結果を整理してみよう。
なお、ここで使うサンプルデータは、NCBI GEO データベースの GSE59468 からデータを間引いて作成したもので、「UNIX 基本コマンド」の講習で用いたものと同じものである（データのダウンロードについては、補足を参照）。
大腸菌の２つの系統を、それぞれ２つの異なる条件で培養して、それぞれについて３つの反復をとった、計12 個のサンプルのデータが含まれている。
12 個のファイルには連番が振られており、各番号はそれぞれ以下の系統・条件の組合せに対応している。
| 番号 | strain | condition | replicate | SRA accession |
|---|---|---|---|---|
| 1 | MG1655 | MPOS | 1 | SRR1515282 |
| 2 | MG1655 | MPOS | 2 | SRR1515283 |
| 3 | MG1655 | MPOS | 3 | SRR1515284 |
| 4 | MG1655 | Urine | 1 | SRR1515285 |
| 5 | MG1655 | Urine | 2 | SRR1515286 |
| 6 | MG1655 | Urine | 3 | SRR1515287 |
| 7 | CFT_073 | MPOS | 1 | SRR1515276 |
| 8 | CFT_073 | MPOS | 2 | SRR1515277 |
| 9 | CFT_073 | MPOS | 3 | SRR1515278 |
| 10 | CFT_073 | Urine | 1 | SRR1515279 |
| 11 | CFT_073 | Urine | 2 | SRR1515280 |
| 12 | CFT_073 | Urine | 3 | SRR1515281 |
以下コマンド先頭の $ はUnix のプロンプトを表すので、それ以降を改行なしで入力すること。
1.
各 *.htseq ファイルには、それぞれのサンプルにおける遺伝子ごとのリード数が記載されている。
この12個のファイルを、以下のようにすべて横に並べて1つのファイルを作りたい。
ファイルを行単位で結合するUNIX コマンドとして、paste がある。以下を実行してみよう。
$ paste results/ecoli.[1-3].htseq | less
paste は引数の順にファイルを結合するので、引数の順番に注意する必要がある。以下のワイルドカードを使った２つのコマンドの出力を比べてみよう。ただし、echo は受け取った引数をそのままの順で表示するコマンドである。
A) $ echo results/ecoli.*.htseq
B) $ echo results/ecoli.?.htseq results/ecoli.1?.htseq
ファイルが1, 2, 3, ...の順に並ぶのはどちらか。なお、コマンドB)は、以下のようにより簡潔に書くこともできる（試してみよ）
$ echo results/ecoli.{?,1?}.htseq
この結果を用いて、paste コマンドによって htseq の出力ファイルを 1, 2, 3, ...の順に結合したファイルを作成し、結果を ecoli.count_all.tmp1 というファイル名で保存せよ。
2.
ecoli.count_all.tmp1 は、b0001 などの遺伝子名が奇数列に繰り返し出現するため、冗長である。そこで最初の列だけ残して、あとの遺伝子名の列は除きたい。すなわち、 1,2,4,6,8,...,22,24 列目のみを残したい。
これは awk を使っても行えるが、cut の方がより簡潔に書ける。cut は -f オプションによって指定された列のみを出力する。たとえば、cut -f 1,3,4 file は、file からタブ区切りで 1,3,4 番目の列を抜き出して表示する。
cut コマンドを使って上記の処理を行い、結果を ecoli.count_all.tmp2 というファイルに保存せよ。
3.
ecoli.count_all.tmp2 の最後の５行は、特定の遺伝子にきちんと対応づけられなかったリードに関する情報が記載されている。tail コマンドで確認しよう。
$ tail ecoli.count_all.tmp2
これらの行はいずれも __で始まっており、例えば、__nofeature は遺伝子領域以外にマップされたリードの数、__ambiguous は複数の遺伝子にまたがる領域にマップされたために１つの遺伝子には対応づけできなかったリードの数を表している。これらの行を、grep コマンドを使って除き、結果を、ecoli.count_all に格納せよ。
結果として、以下のようなファイルができるはずである。
