バイオインフォマティシャンではありませんが、RNAseqで発現解析ぐらいはできないとねと、いろいろガチャガチャLinuxでやってます、やらされてます。日進月歩の世界で、すぐにスタンダードが変わってしまいますが、まあ、慣れるより慣れろという感じでしょうか。
・染色体の名前
良く引っかかったのがゲノム情報上にある染色体番号の表記。>Chr1だったり、>chr1だったり、>1だったりとまちまち。人のプログラムとか使おうとすると、そのあたりでよく引っかかり、sedコマンドで手動でよくいろんなテキストファイルを書き換えてました。
【Linuxコマンド】sedで文字列を置換する方法 | 侍エンジニア塾ブログ（Samurai Blog） - プログラミング入門者向けサイト
igv viewerはエイリアスで、複数の染色体名でも頑張って表記してくれます。これは便利。
https://software.broadinstitute.org/software/igv/LoadData/#aliasfile
ミトコンドリアゲノムとか、気づかずに最後まで解析終わってから、「あれれ、なんかミトコンドリア由来の遺伝子全部０になってる。。。」と気づいてやり直しですからね。シロイヌナズナでは染色体番号がChr1, Chr2, Chr3, Chr4, Chr5, ChrC, ChrMだったり、1, 2, 3, 4, 5, Mt, Pt.だったりして、まちまちです。Plastidか、Chloroplastかで名前が違うのか。
Arabidopsis – Genomics Virtual Lab – Queensland
シロイヌナズナの場合は、TAIR10の次に、Araport11という2016年にアップデートされたアノテーションがあり、これを使いたいわけなのですが、Ensemble plantsでとってきたtair10のゲノムのfastaファイルとは染色体の名称が少々違うようで。。。ゲノム配列にミトコンドリア・葉緑体も含めるのか、含めないのか、カウントする時にどうするのかなど、いろいろ判断が必要なようです。。。しかもこのAraport 11のGTFファイル、RのGenomicFeaturesライブラリのmakeTxDbFromGFFで読ませると、いくつかのTranscriptでストップコドンの位置が転写産物の中にない(stop codons that cannot be mapped to an exon)という警告が返ってきました。makeTxDbFromGRangesのほうだと、drop.stop.codonsのオプションがあるから、エラーにはなりませんでした。
Updated Col-0 Genome Annotation (Araport11 Official Release) Updated Jun 2016 | Araport
・STAR
マッピングは以前Tophat2を使っていましたが、STARに乗り換え。ラボに16CPUのワークステーションがありまして、メモリーも大量に搭載していて、サクっとマッピングは終了します。時間は1/10ぐらいでしょうか、超高速。ところが、データをよく見ると、怪しげなものがいっぱい。。。STARはパラメータの初期設定がヒト寄りになっているのでいろいろ合わない様。下記の通りイントロンの最大サイズが1000kbになってますが、シロイヌナズナのゲノムでは、99％以上が1kb以下です。
--alignIntronMax 1000000 (default)
maximum intron length--alignMatesGapMax 1000000 (default)
maximum genomic distance between mates
Large Introns of 5 to 10 Kilo Base Pairs Can Be Spliced out in Arabidopsis
50bp シングルエンドとリードが短いせいもあって、遺伝子をまたぐような長いイントロンを含むようにマッピングされてました。探してみるほかにも困ってる人がいた！
https://groups.google.com/forum/#!topic/rna-star/cqqaMC_LoDo
ということで、上記のパラメータを小さくしとくのが大事なようです。10000ぐらいでとりあえず満足な結果が得られました。
・トリミング
した方がいいのか、しない方がいいのか。。。Smart seq 2で少量のmRNAをたくさん増幅したせいもあって、プライマーの配列がずいぶんリードの中に含まれています。気にはなるのですが、、、trim_galoreの-a オプションの配列指定(-stringency 5ぐらいで)で取り除いたら乗り除いたでリードが短くなってしまって、ゲノムのあっちこっちにマッピングされるようになるようです（multi loci）。50bp SEですからね。。。マッピングされたデータを見てみるとpoly Tとかpoly Aとかのリードもあっちこっちにマッピングされ、、、ナズナのゲノムには結構Tがつながってるところがあるんですね。。。これもPRINSEQで取り除けますが、illuminaのアダプターもあるし、いったい何回トリムすればいいんやら。。。
・STAR 再び
ということで、何とかならんかなとSTARに戻り、まずは、outFilterMultimapNmaxを1にして、複数個所にマッピングされるようなリードを採用しないことに。まあ、発現解析しようというときは、こっちの方が安全かなあ。
--outFilterMultimapNmax 20 (default)
max number of multiple alignments allowed for a read: if exceeded, the read is consideredunmapped
次に、
--outFilterMatchNminOverLread
default:0.66
float: outFilterMatchNmin normalized to read length (sum of mates’ lengths for paired-end reads)
のパラメータの0.66を0.9にして、マッピングの条件を厳しくしました。トリムしてないリードを直接マッピングすることで、アダプターだろうと、プライマーだろうと、ゲノムにない配列をすべてはじくことにしました(リードの50bp中、40bp以上がゲノムにマッピングされていないといけないので、実質的にはじく)。マッピングされたリードの数はだいぶ減ってしまいましたが、それでもデータ解析をしてみると、結果の傾向は変わらずで、igv viewerで見ても変なピークが減ってだいぶきれいになったので、とりあえずこれで進めることに。ちなみに、75bpほどリードがあると、プライマー配列を除いても十分にマッピングできるほどの配列があるため、このパラメータは変更なしの方がよさそうです。
たかが発現解析、されど発現解析。原理を理解して、生データ見て、試行錯誤しながら条件だしして、きりがないからどっかでケリをつけて先に進む感じは結局ほかの実験と大して変わらないですね。。。
追記）・トリミング再び
プロのバイオインフォマティシャンに相談すると、トリミングはCutadaptの方が便利とのこと。複数のオプションもいっぺんに処理できるし、マルチコア対応で処理も早い。Smart-seq2のプライマー配列と、その相補配列、poly T, poly A, Illuminaのアダプター配列をいっぺんにトリミングしたところ、まあまあきれいになりました。マッピングではじいてもよいけど、それはやはり本筋ではないようで。。。おかげでマップされたリード数も増えました。
・トランスジーンにマッピングする
遺伝子組み換え生物の発現解析なので、トランスジーンにマップされるリードも数えたいところ。ベクターの配列から、mRNAが作られる配列のところのFASTAファイルを作り、ゲノムインデックスを作らせてSTARでマッピングしましたが、Segmentation fault (core dumped)のエラーに。Googleでいろいろ調べてみると、サイズが小さいゲノムの場合は、index作るときにパラメータを変えないといけなかったようです。。。知らなかった。bowtie 2の時はそんなんなかった。。。
For small genomes, the parameter--genomeSAindexNbases needs to be scaled down, with a typical value of min(14, log2(GenomeLength)/2 - 1).
1000bp満たないので、4にして、genome indexを再度作りなおしたらうまくいきました。超高速のソフトウェアだから、一瞬で終わるのかなと思いきや、結構時間かかって、bowtie 2と大して変わらなかったような。。。アルゴリズムって不思議です。
・情報を共有する。
マッピングされたデータを共有したい時があります。すでに離れた人が残したデータを解析するとなると、やっぱりその人にも見てもらいたいわけで。かといって、元ファイルアップロードして、igv viewerをインストールしてもらって。。。というのはハードルがまだ高いです。
そういう時に役に立ったのが、UCSC Genome Browserでした。まず、ヨーロッパにミラーのサーバがあって（上記リンク）、ファイルのアップロード・読み込みが早い。対応ゲノム一覧の系統樹には植物は現れませんが、"enter species or common name"のボックスにaraThaとタイプするとシロイヌナズナのゲノムがちゃんと出てきます。なんだか隠しコマンドみたい(笑)。
Viewerに飛んで(go)、下のところの「add custom tracks」のボタンをクリックすると、自分のファイルをアップロードできるようになります。ブラウザ経由でアップロードしてると遅いし、いろいろ止まりそうになるので、僕の場合は、ハイデルベルク大学のファイルサーバー(heiBOX)にアップロードして、公開リンクを取得して、そのリンクを貼り付け方法で行けました。数十MBのbedgraphファイルが10個ぐらいあってもへっちゃらの様です。その後、genome browserのアドレスをメールすれば、ほかの人にもデータ見てもらえました。
もちろん永久保存というわけではなく、何日か誰もデータにアクセスしないと、custom tracksは自動的に削除されるようです。
目下、NGS解析した結果をまとめて論文執筆中。その手の専門家からみたらトンチンカンなことやってるんじゃないかとひやひやです。。。まあ、とりあえずはやったことを人よりは丁寧に説明して、あとは査読者・読者の判断を仰ぐ感じですかね。バイオインフォマティクスの方は人材不足の様で、あっちこっちで求人広告見ますから、食い扶持なくなったら、そっち方面で生きていくのもありかもですね。思った以上に微調整がいろいろ大変で、こりゃ人手かかりますわ。
