PDF版はこちら

先に ~/gitc/data/HN に移動せよ

1. etec_1.cut.fq, etec_2.cut.fq をsingle end readのデータと見なしてbowtie2でetec をリファレンスとしてマッピングし、結果をファイルetec_bowtie2_single.sam に出力せよ。その際、リードファイルはカンマ区切りで複数指定できることを使え。

2. 出力ファイルの行数を、etec_bowtie2.samと比較せよ。

3. それぞれのファイルの先頭20行をheadで出力して比較し、以下の点について違いを論ぜよ。

   a) ヘッダ行

   b) 出力されるリードの並び

   c) フラグ値

   d) リードがマップされた位置

   e) ペアとなるリードがマップされた位置

1. これらのオプションはどういう意味を持っているか。
2. このコマンドを、出力ファイルをetec_bowtie2_X200.sam として実行せよ。
3. 出力ファイルの行数は、etec_bowtie2.samと比べて変化したか。
4. ファイルの内容を以下のコマンドで比較し、どこが変わったか検討せよ。ただし、diffは2つのファイルを行ごとに比較して異なる行を出力するコマンドで、'<'で始まる行は最初のファイル、'>'で始まる行は2番目のファイルのみに出現する行を示す。また、less の-Sオプションは、長い行を折り返さずに表示することを指示する。

$ diff etec_bowtie2.sam etec_bowtie2_X200.sam | less -S

reference  ATGA-TGGTGTCGA
read       ATGGGTGGAG--GA

1. トランスクリプトNM_001145811.1に関する行のみを抜き出し、sox6_tr1.gtfとして保存せよ。以下2)-4)はこのファイルを対象に調べよ。
2. このトランスクリプトにはいくつのexonが含まれているか。
3. このトランスクリプトに関する情報の各行を、開始位置（リファレンス配列上左端の位置）が転写される向きの順に並ぶように並べ替えよ。
4. このトランスクリプトのCDSの長さの和を計算せよ。
5. sox6.gtfには何種類のトランスクリプトが含まれているか。（ヒント：transcript_id のカラムを抜き出し、ユニークな行の数を数える。ユニークな行はsort -u を用いて抽出できる）

1. samtools stats コマンドは、BAMファイルから様々な統計情報を取り出すことができる。このコマンドを etec_bowtie2_sorted.bam に対して実行し、結果をlessで見てみよう。この結果には様々な情報が含まれているが、行の先頭の文字列で情報を区別して取り出すことができるようになっている。ここでは先頭が COV で始まる　coverage（リファレンス配列の各位置にリードが何重にマップされたか）の分布の情報を取り出そう。Unixコマンドを使った取り出し方は、結果ファイルの中に書かれている。この結果をetec_bowtie2_cov.txtというファイルに保存しよう。

2. 次に、これをscpを使ってローカルの計算機に転送しよう。新しいターミナルを開き、必要に応じてディレクトリを作ってそこに移動する。scpコマンドで、bias5上のファイル（ホームディレクトリからのパスで指定）をカレントディレクトリ(.)に転送する。

3. Rを立ち上げて、転送したファイルを読み込み、２カラム目(coverage)と３カラム目(position数)のプロットを作成しよう。プロットのtypeはlineで、Y軸は対数をとって表示すること。

1. エディタを使って以下のスクリプト(exec_bowtie2.sh)を作成する。

   #!/bin/sh
   #PBS -l ncpus=4
   cd ${PBS_O_WORKDIR}
   bowtie2 -p ${NCPUS} -x etec -1 etec_1.cut.fq -2 etec_2.cut.fq -S etec_bowtie2_2.sam
   

2. qsubで実行する。

   $ qsub exec_bowtie2.sh
   

3. qstat -u ユーザ名 を実行する。何も出力されなければ、ジョブは終了している。

1. bowtie2を実行するコマンド $ bowtie2 -x etec -U etec_1.cut.fq,etec_2.cut.fq -S etec_bowtie2_single.sam

2. 行数のカウント $ wc etec_bowtie2.sam etec_bowtie2_single.sam

行数はどちらも100009 行で同じ。なお、行数のうち9行はヘッダ行、残りの100000行がマッピング結果で、各リードについて必ず1行のマッピング結果がある。

3. 各ファイル先頭20行の表示 $ head -20 etec_bowtie2.sam etec_bowtie2_single.sam

headで先頭20行を表示した結果から、etec_bowtie2.sam (以下pairedと呼ぶ)と、 etec_bowtue2_single.sam (以下singleと呼ぶ)との間に以下のような違いが観察される。

a) ヘッダ行は@PG行のCL（コマンドライン）のみが異なり、あとは同じ。

b) 1カラム目のリード配列名が、pairedでは各リードにつき2行続けて出力されるのに対して、singleでは1行ずつしか出力されていない。pairedでは2つのファイルが同時に読み込まれ、対応するリードが対として扱われているのに対して、singleでは各ファイルが独立なものとして順次処理されている。

c) 2カラム目のフラグは、singleでは0, 4, 16の値をとるのに対し、pairedでは89, 73, 133などの値をとっている。シングルの場合は、1となるフラグは4(セグメントがマップされなかった)または16(逆鎖にマップされた)のいずれかのみであるのに対し、ペアの場合にはより多くの情報が格納されるため、異なる値となる。

d) 4カラム目（マップされた位置）は、異なっている場合と同じ場合とがある。5カラム目（マッピングクオリティ;MAPQ）が42になっているときには4カラム目が同じになっている点に注意しよう。MAPQが高い場合は、ユニークにマップされたことを意味しており、位置はつねに同じになるが、MAPQが低い場合は実際には複数箇所にマップされており、その中の一つがランダムに選ばれている。そこで、ペアで照合する際に別の位置にマップされたと考えられる。なお、相補鎖がマッチした場合は、10カラム目は相補鎖の配列が、11カラム目はクオリティ値が逆向きに表示されている。

e) 7-9カラム目の、ペアとなるフラグメントに関する情報が、singleの方では出力されない（すべて* 0 0となっている）。

1. オプション -I 100 -X 200 は、リード対をマップしたときのフラグメント長が100から200の 間の値であることを指示する（デフォルトは0から500）。

2. bowtie2を実行するコマンド

$ bowtie2 -x etec -1 etec_1.cut.fq -2 etec_2.cut.fq -S etec_bowtie2_X200.sam -I 100 -X 200

3. $ wc etec_bowtie2.sam etec_bowtie2_X200.sam

行数はいずれも100009行で同じ。すなわち、リード対についての条件を変えてもデフォルトではすべての行が出力されるので、行数は変わらない。条件を満たすかどうかはフラグの値で表される。

4. 2つのSAMファイルの違いを表示
   $ diff etec_bowtie2.sam etec_bowtie2_X200.sam | less -S

diffコマンドで表示したetec_bowtie2.sam (以下defaultと呼ぶ)と、 etec_bowtie2_X200.sam (以下X200と呼ぶ)とで異なる行について、一般的に以下のような特徴が観察される（若干の例外はある）。

• 異なる行においては、2カラム目（フラグ）の値がX200の方がdefaultより2小さくなっている。ペアリードの間隔や向きが正しくマップされたかどうかは、フラグの2ビット目（2進数の10、すなわち10進数の2）で示される。X200の方が間隔に対する条件が厳しいため、defaultで正しくマップされたと判定されたものが、X200では正しくないと判定されることがあり、その場合にフラグの2ビット目が1から0に変化した結果、値が2小さくなった。
• 異なる行においては、9カラム目の絶対値（フラグメントの長さ）が200より大きいか100より小さくなっている。そのような場合において、-I 100 -X 200の条件を満たさなくなるのでフラグの値が変化した。

1. $ samtools view etec_bowtie2_sorted.bam ETEC_chr:336-2798 | wc

2. $ samtools view etec_bowtie2_sorted.bam ETEC_chr:4518271-4522299 | wc

数はそれぞれ6個と195個。 マッピングクオリティが20以上という条件を付けるには、コマンドに-q 20を加える。この場合は、結果は変わらない。

$ samtools view -f 13 etec_bowtie2_sorted.bam | wc

ペアリードがあり(1)、自身がマップされていない(4)、相手もマップされていない(8)　フラグは合計13。 このビット全てが立っているフラグ値を持つリードを検索する。数は380。なお、条件を満たすフラグはいろいろあるが、実際に出力される矛盾しないフラグは以下の2通り

01001101 = 77
PAIRED,UNMAP,MUNMAP,READ1

10001101 = 141
PAIRED,UNMAP,MUNMAP,READ2

4M1I5M2D2M

1. $ grep 'NM_001145811\.1' sox6.gtf > sox6_tr1.gtf
   「.」は正規表現で「任意の１文字」を表すので、それを打ち消して「.」という文字のみにマッチさせるために「.」の前に\（バックスラッシュ）をつけている。ただし、この場合はつけなくても結果は変わらない。

2. 15個
   $ grep exon sox6_tr1.gtf | wc

3. $ sort -k 4,4nr sox6_tr1.gtf

4. 2,403 bp
   $ awk '$3=="CDS"{sum+=($5-$4+1)} END{print sum}' sox6_tr1.gtf

5. 4種類
   $ awk '{print $16}' sox6.gtf | sort -u | wc

1. bias5上で実行

   [bias5]$ samtools stats etec_bowtie2_sorted.bam | grep ^COV | cut -f2- > etec_bowtie2_cov.txt
   

2. ローカル計算機で実行

   [local]$ scp [email protected]:data/IU/etec_bowtie2_cov.txt .  （USERNAMEは自分のユーザ名を入れる）
   

3. ローカル計算機でRを立ち上げて以下を実行する。

   > setwd("DIRECTORY")　　（DIRECTORYはファイルを転送したディレクトリ名を入れる）
   > cov <- read.table("etec_bowtie2_cov.txt", sep="\t", header=FALSE) 
   > plot(cov[,2], cov[,3], type="l", log="y", xlab="coverage", ylab="number of positions")