UNIX 講習会演習問題

UNIX 講習会演習問題復習問題１ UNIX 基本コマンド 1. ssh コマンドを使い bias4.nibb.ac.jp にログインして、unix15 ディレクトリに移動しカレン
トディレクトリを確認します。加えてカレントディレクトリ(~/unix15)以下にある全てのフ
ァイルを確認しましょう。オプション"-R"を使うとディレクトリ配下全て再帰的に表示されま
す。
2. ~/unix15/sprot ディレクトリ内の FASTA ファイル全てを~/unixtest/FASTA-15 ディ
レクトリにコピーしましょう。ディレクトリがない場合は新規に作成します。
3. FASTA-15 ディレクトリに移動して、コピーした FASTA ファイル全てのタイトル行のみを
抜き出して less コマンドで確認しましょう。ただし、ここではすべてのファイルに配列は 1
つずつ入っています。確認方法は、head コマンドを使用する方法と、grep コマンドを使用す
る方法の２種類を試してみましょう。それぞれ結果にファイル名が付加されてしまうので、フ
ァイル名を表示しないオプションを man コマンドで調べましょう。
4. 同じように、FASTA ファイル全ての配列行のみを less コマンドで確認しましょう。こちら
は、tail コマンドのみを使用します。
5. ~/unixtest/FASTA-15 にコピーした FASTA ファイルから HUMAN,MOUSE,DROME に関す
る ( ファイル名に含まれる ) Multi-FASTA 形式のファイルをそれぞれ human.fasta,
mouse.fasta, drome.fasta として~/unixtest/の下に cat コマンドを使用して作成しま
しょう。
6. FASTA-15 ディレクトリを削除しますが、削除する前にこのディレクトリ全体の圧縮アーカ
イブを~/unixtest/FASTA-15.tar.gz として tar コマンドで作成しましょう。作成できた
ら、圧縮アーカイブの中身を確認してから FASTA-15 ディレクトリ全体を削除します。
7. 圧縮アーカイブ FASTA-15.tar.gz は自分自身のみ読み書きできるよう chmod コマンドを
使用して(グループとその他は読み書きできないよう)アクセス権を変更します。
8. 圧縮アーカイブ FASTA-15.tar.gz を、保存領域(~/save)へ移動しましょう。
1
復習問題２テキスト処理とスクリプトこの演習では、ecoli.htseq ファイルを使用する。ファイルは下記のパスにある。
~/unix14/text/
ecoli.htseq には、アノテーションテーブルを使って、各遺伝子にマッピングされたリード
数をカウントしたものが書かれている。
１） ecoli.htseq に記載されていて、カウントされた値が１００回以上である行を標準出力
に表示せよ。
２）そのうち、必要なのは b**** という文字列で始まる行である。それらだけを取り出して、
ecoli.temp というファイルに保存せよ。
３） ecoli.temp を、カウントされた回数の多い順にソートし、ecoli.htseq.sorted とい
うファイルに保存せよ。
４）１～３）の操作を一つのコマンドでできるように、htsort.sh というシェルスクリプト
を作成せよ。引数は１つとし、下記のようなコマンド書式でファイル名を与えられるよ
うにし、またソート済みのデータはファイルではなく標準出力に出力すること。
$ ./htsort.sh ecoli.htseq
５）４）のスクリプトにおいて、実行後に中間ファイル（ecoli.temp）を消したい。どうす
ればよいか。
2
復習問題３ qsub を使ったスクリプト講義「スクリプト 2」の最後に作成した「for 文を使った sam → bam 変換」スクリプトを、
下記のように改変してください。
・for でなくアレイジョブを使う
・sam -> bam 変換し、できた bam ファイルをソート
・ソートされた bam ファイルに対してインデックスを作成
・ファイル名は任意につけて別名保存すること（script3.sh を上書きしない）
for 文を使った場合と、アレイジョブを使った場合とで決定的に違うのは何か
qstat の結果から考えてみてください
3
実践演習１ RNA-Seq 解析結果の集計大腸菌の RNA-Seq の例題をもう一度考える。~/unix14/rnaseq に移動しよう。この下の
results ディレクトリには、大腸菌ゲノムにマッピングした結果から htseq-count を使って
各遺伝子領域に重なるリード数を集計した結果(*.htseq)が格納されている。この結果を基に、
UNIX コマンドを使って、簡単な解析結果の集計を行ってみよう。
なお、ここで使ったサンプルデータには、大腸菌の２つの系統を、それぞれ２つの異なる条件
で培養して、それぞれについて３つの反復をとった、計 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
1) 各 *.htseq は遺伝子ごとのリード数が記載されている。この１２個のファイルを、以下の
ようにすべて横に並べて１つのファイルを作りたい。
(ecoli.1.htseq)
(ecoli.2.htseq)
(ecoli3.htseq)
b0001
11
b0001
0
b0001
7
b0002
117
b0002
9
b0002
131
b0003
33
b0003
1
b0003
31
b0004
44
b0004
4
b0004
58
b0005
3
b0005
0
b0005
2
...... (ecoli12.htseq)
ファイルを行単位で結合する UNIX コマンドとして、paste がある。以下を実行してみよう。
$ paste results/ecoli.[1-3].htseq |less
4
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.paste_all に格納せよ。
2) ecoli.paste_all は、b0001 などの遺伝子名が奇数列に繰り返し出現するため、冗長であ
る。そこで最初の列だけ残して、あとの遺伝子名の列は除きたい。すなわち、1,2,4,6,...,22,24
列目のみを残したい。これは awk を使っても行えるが、cut の方がより簡潔に書ける。cut は
-f オプションによって指定された列のみを出力する。たとえば、cut -f 1,3,4 file は、
file からタブ区切りで 1,3,4 番目の列を抜き出して表示する。
また、ecoli.paste_all の最後の５行は遺伝子領域にマッチしなかったリード数が記載され
ているが、これらの行を除きたい。grep コマンドを使って除くことを考えよ。これらの処理を
行った結果を、ecoli.count_all に格納せよ。
3) 発現解析でよく行われる標準化の方法に RPKM (Read Per Kilobase per Million mapped reads)が
ある。これを計算するには各遺伝子の長さの情報が必要である。遺伝子アノテーションファイ
ル ecoli.gtf には各遺伝子（exon）の開始位置と終了位置の情報が含まれている。awk を使
ってこれから遺伝子の長さを計算しよう。3 列目が exon である行について、開始位置（4 列目）
と終了位置（5 列目）を使って長さを計算し、10 列目にある遺伝子名ともに出力する。
正しく実行すると、以下のような出力が得られるはずである。 "b0001"; 66
"b0002"; 2463
"b0003"; 933
"b0004"; 1287
"b0005"; 297
.....
ダブルクォート(“)とセミコロン(;)が余分なので、取り除こう。これには sed コマンドが使え
る。sed はファイルの簡単な編集を行うプログラムで、‘s/置換する文字列/置換後の文字列/g’
によって、ファイル中の指定した文字列を一斉に置換することができる。置換後の文字列を空
にすることで、特定の文字列を除去することもできる。また、置換する文字列には正規表現が
使えるので、〔ダブルクォートまたはセミコロン〕という指定は正規表現を使って書ける。そこで、以下のコマンドによって、file 中のすべてのダブルクォートとセミコロンを除くこと
5
ができる。 $ sed ‘s/[“;]//g’ file
これによって、以下のような出力が得られるはずである。
b0001 66
b0002 2463
b0003 933
b0004 1287
b0005 297
....
これをファイル ecoli.gene_length に保存せよ。
4) ecoli.count_all と ecoli.gene_length を結合して一つのファイルにしよう。実は、
ecoli.count_all は遺伝子名でソートされているが、ecoli.gene_length は遺伝子名でソ
ートされていない（ecoli.gtf が位置によってソートされているため）。従って、そのまま
paste すると正しく結合されない（確認せよ）。そこで、まず ecoli.gene_length を遺伝子
名によってソートし、ecoli.gene_length_sorted というファイルに保存しよう。
ecoli.count_all と ecoli.gene_length_sorted を paste コマンドで結合することもで
きるが、また遺伝子名の行が余分にできてしまう。ここでは別のコマンド join を使った結合を
してみよう。join は２つのファイルに共通の列（ここでは１列目の遺伝子名）があり、かつフ
ァイルがいずれもそれらの列に関してソートされているとき、それらの列をキーとして２つの
ファイルを対応づけて結合するコマンドである。
$ join -1 1 -2 1 ecoli.gene_length_sorted ecoli.count_all
引数の -1 1 -2 1 は１番目のファイルの１列目と２番目のファイルの１列目をキーとして対
応づけることを意味しているが、これはデフォルトの動作なので、省略しても同じ結果になる。
paste はキー（ここでは遺伝子）の並びが２つのファイルで完全に一致していないと、行が合
わなくなって意味がなくなるが、join はキーに重複や抜けがあって２つのファイルの間でずれ
ていても、同じキーを持つ行を正しく対応づけてくれる。非常に強力なコマンドである。
なお、join はタブ区切りのファイルをスペース区切りに変更してしまう。見やすさを維持する
ためにタブ区切りに直したい場合は、sed ‘s/ /
/g’ file によって、タブ区切りに変換
できる。ただし、sed コマンドの最初の文字列はスペースで２番目はタブである。コマンドラ
インでタブを入力するには Control-v を押してから Tab キーを押す。または、この場合タブ
の代わりに¥t を入れても同じ意味になる。
これによって以下のようなファイルができるはずである。
6
b0001
66
11
0
7
4
b0002
2463
117
9
131
b0003
933
33
1
31
7
17
0
0
3
51 20 161
20 1
0
32 34 2
10 4
9
2
7
30
9
1
7
9
1
4
.....
この結果を ecoli.result_all に保存しよう。
5) ecoli.result_all のようなタブ区切りテキストは、Excel などの表計算ソフトや R などの
統計解析ソフトに読み込んで処理できる。これ以降はおそらくそのようにした方が早いが、
UNIX の練習のため、もう少し UNIX コマンドを使った解析を続けよう。
RPKM を計算するには、あと各サンプルでマップされたリード数の和を計算する必要がある。
これも awk で計算することはできるが、やや面倒なのでここでは省略しよう。
ecoli.result_all は 1 行目に遺伝子名、2 行目に遺伝子の長さが入っており、3,4,5 列目、
6,7,8 列目、9,10,11 列目、12,13,14 列目が、それぞれ４つの条件の組合せについての 3 回ずつの
反復実験の結果に対応している。この 3 回の反復実験の結果を平均し、さらに遺伝子の長さで
割った上で 1000bp あたりの値に標準化しよう。1 列目に遺伝子名、2, 3, 4, 5 列目に各条件の反
復実験の平均値を入れる。この処理を、awk を使って行い、結果を ecoli.result_mean に
格納せよ。
7
実践演習２ TopHat/Cufflinks を使った RNA-Seq 解析の準備以下は、 Trapnell et al. Nat. Protoc. 7, 562-578 に掲載された、TopHat/Cufflinks の実行手順の最初
の部分を Bias 上で実行する際の手順を示している。以下のサイトからこの論文を取得すること。
http://www.ncbi.nlm.nih.gov/pubmed/22383036
また、この実習は/scratch ディレクトリ以下で行う。このディレクトリ上のファイルは 1 ヶ
月後に自動的に消去されるので、作成したスクリプトなど、残しておきたいファイルはホーム
ディレクトリ上にコピーしておくこと。
1) ディレクトリ /scratch/ユーザ名に移動し、さらに test_tophat ディレクトリを作成
してその下に移動しよう。
以下のファイルを使う。
/usr/local/data/unix14.extra/GSE32038_simulated_fastq_files.tar.gz
このファイルのシンボリックリンクをカレントディレクトリに作成せよ（ln –s コマンド）。
これを tar コマンドでカレントディレクトリに展開せよ。全部で 12 個の FASTQ ファイルが作成
される。
2) ファイル名から GEO のアクセッション番号（GSM794483 など）を取り除いたシンボリック
リンクを作成しよう。
$ ln –s GSM794483_C1_R1_1.fq.gz C1_R1_1.fq.gz
のようなコマンドを 12 回繰り返してもよいが、ここではスクリプトで処理してみよう。
ファイル名からアクセッション番号を取り除くには、
$ echo GSM794483_C1_R1_1.fq.gz | sed 's/GSM[0-9]*_//'
とするとよい（実行してみよ）。以下のスクリプトはこれを使ってすべてのファイルのシンボリ
ックリンクを作成する。
for f in GSM*
do
newf=`echo $f | sed ‘s/GSM[0-9]*_//’`
ln -s $f $newf
done
3) Illumina iGenomes は、様々な生物種のゲノム配列のアセンブリと、マッピングコマンド用の
インデックス、およびアノテーションをパッケージ化して、すぐに解析を始められるようにし
たものであり、Bias 上では/bio/db/igenomes 以下に格納されている。まず Drosophila のゲノムデータと Bowtie2 用のインデックスについて、カレントディレクトリに
8
シンボリックリンクを作成する。 $ ln –s /bio/db/igenomes/Drosophila_melanogaster/Ensembl/BDGP5.25/Seque
nce/Bowtie2Index/genome.* .
コマンドは改行せずに 1 行で入力すること。ファイル名が非常に長いが、Tab キーを適宜押して、
ファイル名補完の機能を使って効率的に入力する。また、コマンドの最後のドット（.）はリン
ク先がカレントディレクトリであることを表すので、忘れないこと。同様にして、Drosophila の遺伝子アノテーション（GTF ファイル）について、カレントディレク
トリにシンボリックリンクを作成する。 $ ln –s /bio/db/igenomes/Drosophila_melanogaster/Ensembl/BDGP5.25/Annot
ation/Genes/genes.gtf .
4) これで論文 570 ページ以下の PROCEDURE を実行する準備が整った。まず６つのサンプルにつ
いて tophat を実行するが、Bias では大きな計算ジョブはコマンドを直接打って実行するのでは
なく、qsub を使って実行させるのが基本である。ここでは、アレイジョブを使って 6 本並列に
走らせてみよう。また、tophat と cufflinks は一連の処理として直列に実行させるものなので、一
つのジョブにまとめよう。
以下のスクリプトでは、6 種類のサンプルの名前(Condition 1, 2 と Replicate 1, 2, 3 の組合せで
つけられている)を配列変数 $files に格納し、各ジョブはこの配列を参照して異なるファイル
をとって並列実行している。配列の添字は 0 から、アレイジョブの変数 $SGE_TASK_ID は 1
から始まるので、添字は$SGE_TASK_ID から 1 引いている。
tophat、cufflinks ともに-p オプションを指定して、8 スレッドを使ったマルチスレッドで実行し
ている。このことを SGE に伝えるため、qsub のオプションとして-pe smp 8 を指定している。
#!/bin/bash
#$ -t 1-6
#$ -pe smp 8
#$ -cwd
files=(C1_R1 C1_R2 C1_R3 C2_R1 C2_R2 C2_R3)
fname=${files[$SGE_TASK_ID-1]}
tophat -p 8 -G genes.gtf -o ${fname}_thout genome ${fname}_1.fq.gz ${fna
me}_2.fq.gz
cufflinks -p 8 -o ${fname}_clout ${fname}_thout/accepted_hits.bam
このスクリプトを作成して qsub によって実行せよ。ただし、7 行目から始まる tophat コマン
ドは、次の行と改行なしで続けて打つこと。
ここまでで PROCEDURE の 2 までが終わっている。結果が確認できたら、続けて 3 以降を実行し
ていこう。
9
演習問題解答編復習問題１ UNIX 基本コマンド 1.
$
$
$
$
ssh [email protected]
cd unix15
pwd
ls -R
2. 以下は、先に~/unixtest ディレクトリに移動してから操作する場合となります。
$ cd ~/unixtest
$ mkdir FASTA-15
$ cp ~/unix15/sprot/*.fasta FASTA-15
3.
$
$
$
$
$
cd FASTA-15
man head
man grep
head -1 -q * | less
grep -h '^>' * | less
4. tail コマンドの引数に "-n +行数" を指定すると、先頭から数えて、指定された行数以
降を表示します。
$ tail -n +2 -q * | less
5.
$ cat *HUMAN.fasta > ../human.fasta
$ cat *MOUSE.fasta > ../mouse.fasta
$ cat *DROME.fasta > ../drome.fasta
6. tar でアーカイブを作成する場合、対象ディレクトリのひとつ上で対象ディレクトリを指定
すると、解凍した際に対象ディレクトリが作成されてその下にファイル群が配置される。
反対に対象ディレクトリ内でカレントディレクトリ(.)を指定すると、解凍時にカレントディレ
クトリ内にファイル群が作成されてしまうので注意すること。
$
$
$
$
cd ..
tar zcvf FASTA-15.tar.gz ./FASTA-15
tar ztvf FASTA-15.tar.gz
rm -rf FASTA-15
7. chmod u+rw,g-rwx,o-rwx FASTA-15.tar.gz でもよい。
$ chmod 600 FASTA-15.tar.gz
8.
$ mv FASTA-15.tar.gz ~/save
10
復習問題２テキスト処理 1) $ awk '$2>=100{print}' ecoli.htseq
2) $ awk '$2>=100{print}' ecoli.htseq | grep ‘^b’ > ecoli.temp
3) $ sort –k 2,2nr ecoli.temp > ecoli.htseq.sorted
4)
filename=$1
awk ‘$2>=100{print}’ $filename | grep ‘^b’ > ecoli.temp
sort –k 2,2nr ecoli.temp
5) 最後に rm ecoli.temp を加える。
復習問題３ qsub を使ったスクリプト #!/bin/sh
#$ -cwd
#$ -t 1-12
fn=`basename result/ecoli.${SGE_TASK_ID}.sam .sam`
samtools view -bS result/${fn}.sam > result/${fn}.bam
samtools sort result/${fn}.bam
result/${fn}_sorted
samtools index result/${fn}_sorted.bam
11
実践演習１ RNA-Seq 解析結果の集計 1) ファイルが 1,2,3,...の順に表示されるのは B)の方。 $ paste results/ecoli.{?,1?}.htseq > ecoli.paste_all
2) grep –v はマッチする行を除いて出力する。
$ cut –f 1,2,4,6,8,10,12,14,16,18,20,22,24 ecoli.paste_all | grep -v
‘^_’ > ecoli.count_all
3) 長さは（終了位置）—（開始位置）＋１で計算される
$ awk ‘$3==”exon” {print $10, $5-$4+1}’ ecoli.gtf | sed ‘s/[“;]//g’ >
ecoli.gene_length
4) 行の先頭からアルファベット順に並べたい場合は sort のオプションは不要。また、join の
キーがどちらも 1 列目であるときは join のオプションは不要。
$ sort ecoli.gene_length > ecoli.gene_length_sorted
$ join ecoli.gene_length_sorted
ecoli.count_all | sed ‘s/ /¥t/g’ >
ecoli.result_all
5)
$
awk
‘{
print
$1,
($9+$10+$11)/3/$2*1000,
($3+$4+$5)/3/$2*1000,
($12+$13+$14)/3/$2*1000
>ecoli.result_mean
12
($6+$7+$8)/3/$2*1000,
}’
ecoli.result_all

Download Report