やったこと
1. velvet を使って、 mRNA-seq のシーケンスデー
velvet_1.2.03 ディレクトリに移動してコンパイルを行う。
n0rr:velvet nor$ cd velvet_1.2.03
n0rr:velvet_1.2.03 nor$ make MAXKMERLENGTH=63
MAXKMER-LENGTH=63のオプションを入れることによって、元々の最
大kmer 31を63に変更することが出来る。
これが完了すると、
n0rr:velvet_1.2.03 nor$
/Users/nor/analtools/velvet/velvet_1.2.03/velveth と打った時にvelveth の使い方が出てくる。
次にoasesのインストール
velvetのリンクからoasesのページに移動出来る。
http://www.ebi.ac.uk/~zerbino/oases/ ここからダウンロー ド。以下のコマンドで、/Users/nor/analtools ディレクトリ内
に oases ディレクトリを作って、ダウンロードした
oases_0.2.05.tgz をoases ディレクトリ内に移動させる。
n0rr:velvet_1.2.03 nor$ cd ..
n0rr:velvet nor$ cd ..
n0rr:analtools nor$ mkdir oases n0rr:analtools nor$ cd oases/
n0rr:oases nor$ mv /Users/nor/Downloads/oases_0.2.05.tgz ./oases_0.2.05.tgz
もちろん移動はFinderでやってかまわない。
n0rr:oases nor$ tar zxvf oases_0.2.05.tgz 解凍して、
n0rr:oases nor$ cd oases_0.2.05
oases_0.2.05 ディレクトリ内に移動しコンパイルを行う。
n0rr:oases_0.2.05 nor$ make
'VELVET_DIR=/Users/nor/analtools/velvet/velvet_1.2.03/' 'MAXKMERLENGTH=63'
ここで、velvetのディレクトリを指定し、また、最大kmer長
の変更について宣言する。アセンブリしたいシーケンスが 入っているディレクトリ移動し、
n0rr:20120109run_Sample_lane3 nor$
/Users/nor/analtools/velvet/velvet_1.2.03/velveth ./mac_vel-vet63 63 -fastq s_3_1_sequence.txt
n0rr:20120109run_Sample_lane3 nor$ cd mac_velvet63/
n0rr:mac_velvet63 nor$
/Users/nor/analtools/velvet/velvet_1.2.03/velvetg ./
-read_trkg yes -scaffolding no
つづいてoasesを使う。 oasesはスプライスバリアントなどの
多型もまとめたmRNA-seq解析用の拡張ツールであり、
transcripts.faを出力する。
n0rr:mac_velvet63 nor$
/Users/nor/analtools/oases/oases_0.2.05/oases ./
ここまでで、contigs.fa とtranscripts.fa が得られた。これら はショートリードをアセンブリして得られたcDNAのリスト になっている。これらnファイルの統計をとってみよう。まず はcDNA数を数えてみる
nor$ more contigs.fa | grep ">" | wc -l
それぞれ、65111本と29764本の配列が入っていることがわ かった。続いて、配列長の統計を取ろう。awk の”length” を 用いる。
nor$ more contigs.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr ">"
"\n" | awk -F "@" '{print length($2)}' | head -n4
ここで、head -n 4 というのはあたま４行だけ出力させること
を示す。塩基長測定がうまく行ってるか調べてみよう。
nor$ head -n 4 contigs.fa
としてファスタのヘッダーを確認する。
最初の０はファスタの”>” が”\n”になったものなので無視す る。次の1176とファスタのヘッダーのlength_1114とは食い 違っている。実際に正しいのは1176である。というのも、
ヘッダーに書いてある長さは実際の長さよりkmer-1だけ短く なっているからである。このことを確かめるために最も短い
contigを見てみよう。今回もっとも短いのはlength_63のもの
である。確かめるには
nor$ grep length_62_ contigs.fa
とすればよい。62だと出てこないが、63にすると出てくる。
最も短いcontigをみてみる。
nor$ more contigs.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr ">"
"\n" |grep _length_63_ | awk -F "@" '{print length($2)}' | head -n4
nor$ more contigs.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr ">"
"\n" |grep _length_63_ | awk -F "@" '{print ">"$1"\n"$2}' | head -n4
納得していただけただろうか。表示されている文字をecho コ マンドで出力して数えてみてもいい。
nor$ more contigs.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr ">"
"\n" |grep _length_63_ | awk -F "@" '{print ">"$1"\n"$2}' | head -n4
nor$ echo
"CTGCAGCAGTTAAACCAAAAGCCTCGGCAAAAGCTCCAGG AAAAGGTGTTAAGAAACCAACAGCAAAGGTGCTGCGAAAC
CTGCAGCAGTTAAACCAAAAGCCTCGGCAAAAGCTCCAGGA AAA" | wc
これでこの方法でうまく行くことがわかった。早速出力して Rに読み込ませよう。今回は手を抜いて荒いやり方をする。
nor$ more contigs.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr ">"
"\n" | awk -F "@" '{print length($2)}' | tr "\n" "," >
length_contigs.txt
nor$ more transcripts.fa | sed -e '/>/s/$/@/g' | tr -d "\n"| tr
">" "\n" | awk -F "@" '{print length($2)}' | tr "\n" "," >
length_transcripts.txt
これをテキストエディタで開いて編集する。開くと 0,1176,1372,... と続いている。ここに、
”lcon<-c”を追加する。0は消してしまう。
さらに、ファイルの一番下に行って、”,” を消して”)”を置く。
こうすることで、
lcon<-c(1176,1372,1125,2877,205,...,125,125)
というテキストができ、これはこのままRにコピーペーストで きる。すべてを選択し、
コピーしてRにペーストする。
summary(lcon) で要約をみよう。
hist(lcon)
でヒストグラムが描ける。
hist(log(lcon))
で対数をとるのもいい。
MASSパッケージを読み込むと違った書き方もできる。
library(MASS) truehist(log(lcon))
contigs.fa とtranscripts.faを比べてみよう。transcipts.faの長 さデータ、length_transcripts.txt を同様に処理した。
以上である。