ハーバード留学研究2年目編010：ちょいしょぼNGS解析覚書2

そんなに重要な解析ではないが、QCの一環として、それぞれ別々の蛍光タンパク（Tdimer2とCerulean）を発現する細胞から、とってきたRNA-seqデーター(74.bam, 75.bam)から本当に前情報と同じ蛍光タンパクを発現しているかどうかチェックしてみる。

前回行った解析手法を利用してみる。Tdimer2とCeruleanのfastaqファイル及びindexファイル作成は前回に準じる。

１）BAM fileからunmapped sq抽出
$ samtools viwe -f4 -b 74.bam | samtools sort -74-unmp
$ samtools viwe -f4 -b 74.bam | samtools sort -74-unmp

$ samtools index 74-unmp.bam
$ samtools index 75-unmp.bam

2)Tdimer アラインメント

$ bwa aln -b Td.fa 74-unmp.bam -> 74_T.sai
$ bwa aln -b Td.fa 75-unmp.bam -> 75_T.sai

3)SAI -> SAM変換

$ bwa samse Td.fa 74_T.sai 74-unmp.bam > 74_T.sam
$ bwa samse Td.fa 75_T.sai 74-unmp.bam > 75_T.sam

4)SAM->BAM変換

$ samtools view -bS 74_T.sam > 74_T.bam
$ samtools view -bS 75_T.sam > 75_T.bam

5)ソート
$ samtools sort 74_T.bam 74_T_sort
$ samtools sort 75_T.bam 75_T_sort

6)インデックス作成

$ samtools index 74_T_sort.bam
$ samtools index 75_T_sort.bam

7)結果表示

$ samtools idxstats 74_T_sort.bam
$ samtools idxstats 75_T_sort.bam

てな具合である。
主力結果は

74が
Tdimer 1395 1117 17
* 0 0 15976919

75が
Tdimer 1395 0 0
* 0 0 17205327

のようにでてくる。
このうち2つ目のカラムの上が、マップされたリード数だから（詳しくはsamtools使い方参照）、74のunmapped sequence のうち、1117リードがTdimerの配列にマップされた、一方で75のunmapped sequence のうちTdimerの配列にマップされたものはなかったということになる。

Ceruleanへのアラインメント及び結果表示も同様。

74が
Ceru 825 7 0
* 0 0 17205454

75が
Ceru 825 4451 114
* 0 0 15972354
となっている。74のほうは4リードが、75のほうは4451リードがunmapped sequenceからCeruleanの配列にマップされたことになる。

定性的ではあるが、これで74がTdimer, 75がCeruleanを発現していることがわかり、前情報とサンプルの中に発現している蛍光タンパクは同じであると考えられる。ふー。。

ハーバード留学研究2年目編009：BWAとBowtieどちらがよい？

ちょっと最近、RNA-seqアナリシスがBWAをつかってうまくいきそうだったので、Bowtieから乗り換えようかと思ったりして、悩んでいた。

http://arxiv.org/abs/1301.5187のサイトに面白いことが出ていて、

Bowtie2GP > Bowtie > Bowtie2 > BWAという感じらしい。



Bowtie2使うのが一番賢そうな感じだろうか？

ハーバード留学研究2年目編008：懲りずにおこなうRNA-seqデータ解析覚書。

共同研究でRNA-seqデーター解析をしており、性懲りもなくデーターをもらってきて勉強がてら解析してみることにした(*)。

Pair-endのデーターなので

> bowtie -S -p 2 mm9 -1 X_1.fastq -2 X_2.fastq X.sam (**)

でいいはず。
しかし、ほとんどがunmappedになってしまう。

ん？とりあえず、トラブルシュートを検討します。

とりあえず片側だけでやった方が、unmappedがすくなそうなので(***)、

ためしてみる。

> bowtie -S -p2 mm9 X_1.fastq X_1.sam

としunmapped sequenceを抽出してみることとする。

★Sam->Bam変換

> samtools view -bS X_1.sam > X_1.bam

★unmappedシーケンスの抽出

> samtools viwe -f4 -b X_1.bam | samtools sort -X_1-unmp (****)
> samtools index X_1-unmp.bam

と遅くなってきたので、今日はここまで。。
PCが遅い。。(涙）

(2014年9月3日追記）

★多色蛍光マウスなので、CAG－GFPの配列にあててみることとする。

詳しくは、STAP細胞関係のゲノムデータを解析してみた7619にある通りなのだが、

１）CAG-GFPの配列をとってき、fast形式のの配列をつくる（CAG.fa）（ここはSTAP細胞関係のゲノムデータを解析してみた7619を参照のこと）

２）インデックスをつくる。

>bwa index -a bwtsw cag.fa

３）次にアラインメントするが、ここはSTAP細胞関係のゲノムデータを解析してみた7619の通りだとなぜかうまくいかなかった。

そこでお手本のように引数を渡す手法でなく、分割しておこなった。
本当はbowtieが使えるとよいのだけれど、bamファイルでの入力ができない。

まずアラインメント
>bwa aln -b CAG.fa X_1_unmp.bam -> X_1-cag.sai

次にsai -> sam変換
>bwa samse CAG.fa X_1-cag.sai X_1_unmp.bam > X_1-cag.sam

sam -> bam変換

>samtools view -bS X_1-cag.sam > X_1-cag.bam

ソート
>samtools sort X_1-cag.bam X_1-cag_sort

インデックス作成

>samtools index X_1-cag_sort.bam

で無事作業は終わっていると思う。

最後にアラインメントされたところのみをみると、

>samtools view -F4 X_1-cag_sort.bam

のように出力できる。
ということで今日はここまで。

（2014年9月18日＆10月18日）

出力は見にくいので、

>samtools view -F4 X_1-cag_sort.bam > X_1-cag.txt

みたいにファイル出力して、あとでエクセルで見た方がよいかも。

また単にあたっているかどうか、見るだけなら、

>samtools idxstats X_1-cag_sort.bam > X_1-cag.txt

でもよいかもしれない。

列番号 　　　　説明 　　　　　例
1 リファレンスの配列名 CAG
2 リファレンスの配列長 6118
3 マップリード数 　　　　311
4 非マップリード数 0　　

のように何本あたっているかが表示される。詳細はsamtoolsの使い方参照。


(*)いつぞやのSTAPデータでは、メモリ不足でこりました。今回は多少Ubuntuを高速化したりしてリトライです。何とかネットブックでも遊べそうな勢いです。何でSTAPデータはあんなに巨大だったのか。。
(**)Xはサンプル名です。
(***)何かやり方を間違っているのだろう。それでも30％くらいある。
(****)詳しくは以下のサイトが参考になる。

１）STAP細胞関係のゲノムデータを解析してみた7619

２）sam形式ファイルsam format

ちなみにsamtoolsのオプションの説明はsam形式ファイルsam formatを参考に。

-S はinput fileがsam形式であることを指定している。
-f は許可する数字，この場合4なのでマップされていないリードを許可します。
-Fは許可しない数字，この場合194（128+64+2）なのでpaied endでマップされているリードは許可しない

その他のオプションの説明はsamtoolsの使い方
より

Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help

ハーバード留学研究2年目編007：留学開始1年5か月たって

もうそろそろ留学開始1年5か月。。というか一年半。「一丁前な留学」的な期間になってきた（*）

ようやくイメージングの系が何とか軌道に乗り始め、また新たに立ち上げた遺伝子スクリーニング系もそこそこ動き始めてきた。

ここまでくると、ある程度の期間きちんとやれば、何らかのものは得られそうである。このフェーズの研究は、真理の探究といった側面が強くなり、楽しい。

そしてちょっと時間をかけてみたいという欲もでてくるので、何とか来年度のグラントが取れるとよいなと思う。もしくは政治力・交渉力そして英語力が必要か？

1年目の終わりに立てた目標を振り返ってみると、

１）1-6か月目：プラットフォームを作る（スクリーニングとブラッシュアップ）
　　　　　イメージング、バイオインフォマティックス等のツールを確保、
　　　　　ターゲット遺伝子確定など。

２）6-7か月目：奨学金応募。できれば中長期的なfunding確保。

３）6-10か月目：論文の基礎となるFig1-2のデーターを確保。

であった。

１）に関しては前にも述べたように、新規イメージングプローブを2種開発、遺伝子スクリーニング系も立ち上げた。

イメージングプローブに関してはin vivoへの導入も何とか可能そうであるが、当初目標をちょっと下方修正し、あと1-2か月でどんな形のストーリーでどいったデーターを必要とするかを見極められたらよしとしたい。

遺伝子スクリーニング系は、当初考えていたよりはバックが出やすいという欠点はあるが、そこそこ動きそうである。また新しくつくったイメージングプローブとの組み合わせで、ちょっとエレガントな遺伝子スクリーニングも作れそうで、これも楽しみである。

なお日本でさんざんやっていたことだったので、恥ずかしい限りだが、超簡単なジェノミックDNAの回収＆PCRの系の確立に最初手間取ってしまい、トラブルシュートに追われた(**)。ターゲット遺伝子の確定までに1か月でたどりつけるかどうか？これもちょっと下方修正しないといけないかもである。

２）に関しては、とりあえず今年は3つはだせそうである。1個はだしたので、あと2個は最低。また中長期的なものが出せるかどうか、情報収集につとめたい。またここ数か月の間にでるであろうデーターとともにボスもしくはコラボレーターとの相談というカードも考慮にいれていくのがよいかもしれない。

３）については、6月末にあと3か月中には何とかしたいといっていた。期限まであと一か月。何とか目星をつけるところまで行きたい。Fig1-2については年末をめどといったところだろうか？ちょい下方修正。。


(*)医局派遣の留学なんかだと、1年半くらいで帰国されるケースも多い。ここで帰国しても、「あっつ、留学されたんですね。英語ぺらぺらですよねー」と勘違いコメントをされるような気がする。。残念ながら、1年半いても英語はそんなにうまくなりません（笑）

(**)キットを使わなくてよいので、いったん確立してしまえば楽です。

ハーバード留学研究2年目編006：留学開始1年と2か月たって

留学2年目がこの間始まったと思ったら、もう2か月たってしまった。
一年目に取り組んでいたプロジェクトの進みが遅く、第二、第三の矢として始めたプロジェクトも何とか第一段階はクリアしたものの結果を得るまでにまだまだかかりそうだ。

1年目の終わりに立てた目標を振り返ってみると、

１）1-6か月目：プラットフォームを作る（スクリーニングとブラッシュアップ）
　　　　　イメージング、バイオインフォマティックス等のツールを確保、
　　　　　ターゲット遺伝子確定など。

２）6-7か月目：奨学金応募。できれば中長期的なfunding確保。

３）6-10か月目：論文の基礎となるFig1-2のデーターを確保。

であった。

１）に関しては新規イメージングプローブを2種開発、遺伝子スクリーニング系も立ち上げと、第一段階としては順調ではないかと思う。特に新規イメージングローブはボスに気に入ってもらえて、今度出す予算の中のプロジェクトとして組み込んでもらった。

しかしながら、結果についてはシビアに求められる感じで、特にin vivoのデータへの要求が高い。新規プローブのため validationを少し慎重にしているのと、最近分かったのだけれど、ここでよく扱っている動物モデルが思った以上に扱いづらく、in vivoへの導入に時間がかかっているところが現時点の問題である。

また評価してもらっている気もするのだが、ちょっと長めの数年単位のプロジェクトを提案するとあまりいい顔をされないのが微妙である（*）。この辺２）ともかかわるが、手放しで給料をくれるというほどの評価でもないのかとちょっと落胆する。

２）に関してはいくつか出せるものを積極的に探していかないといけないのだけれど、とりあえずは去年落ちたグラントに、テーマを変えて応募してみることにする。新規遺伝子スクリーニングなのだけれど、ちょっと複雑だねとボス受けは悪いのがまた悩みの種である。

３）については、あと3か月中には何とかしたいというか、今までの経験上、何らかの形で第一、第二、第三のどれかは伸びてくると思うので、3か月中には何とかできるのではないかと思っている。

慣れてきた半面、英語で思うようにコミニュケーションが取れないところなど、ほおっておいた弱点の部分がボディブローのように効いてきている気がする。

(*)早く帰れっていうことか。。(笑)