アグリな凡人が奮闘する話

大抵インフォ、時々バイオ、稀にアグリな日記。

泥沼

今日NGS解析に進捗がほんの少し生まれたので書いておく。

まず渡されたbamファイルからunmappedリードを取り出す。

共生特異的なリード結果を無菌のものにマッピングした結果である。

ここからunmappedリードを取り出したことによって実際共生菌のみのリードを取り出すことができた。

なので次にこのunmappedリードを共生特異的ランの全リードデータから引けば、実際にランのみで発現していたリードを取り出すことが可能になる。

なので次にやる作業はbamファイル(リード)同士の引き算。

bamファイル、つまりバイナリファイル同士の比較をすればよいのだが、この比較の方法がいまいちわからん。

コマンドプロンプト上でcmpコマンドを使えば1バイトずつ比較して違う場所を返してくれる。(cmp -l ファイル名 ファイル名)で

でもあくまで違うファイルを比較するというだけで、違う部分を抽出して新しくファイルにし直すだとか、そもそも1バイトずつ比較しているので、全体からその行を検索しているとかのわけではない。

バイナリ形式じゃなければファイル間の比較もできるのか?でもそもそもSAMファイルを置いておけるほどパソコンの容量が無い。cmpやdiffを使う方法は有効と思われたけどやめておこう。

 

tsvファイルをいじっていく方が簡単なのか?

ということで。

共生菌Tsulasnella Colasporaのアノテーションデータがウェブからダウンロードできたので、せっかくなのでアノテーションしたtsvファイルも作ってみようかと思う。

JGI Genome Portal - Home (doe.gov)

 

データをダウンロードした後すべてのゲノムが入ったgffファイルがあったので、これをもとにアノテーションしてみる。

 

tsvファイルの作り方はAtomから作った本の内容をそのまま、ちょっとファイル名変えただけで。

ちなみにダウンロードしたすぐのgffファイルはgz形式で圧縮されていたのでアプリで解凍して使った。

いけるかな、と思ってたけど、甘かった。。。

そもそもTulasnellaの遺伝子のアノテーションCDS,exon,inferred_parentしかない。mRNAがあればいいんだけど、そんなのない。(仮にmRNAがあればgene_idからidを使ってきれいにtsvファイルとして保存できる)

どうしよと思ってしょうがないからこれらの領域使ってやったみたけど当然うまくできない。アノテーションファイル作れなかったら解析うまく進まないしなあ。。。

ということでpythonからTulasnellaのtsvファイルをいじるのはあきらめた方がいいかも。ゲノムが全部FASTAファイルで入ってるからIGVとか使ってマッピングはできるかもしれんが。もうちょい粘ってみる。

 

さて、今後の方針をまとめると

①unmapped readの抽出を再度元データにマッピング

②Zengさんの論文を読んでpythonで適用できないか考える(望み薄)

③メタゲノム解析の手法を取り入れる

④blobtoolsを使ってみる(望み薄)

⑤orchidとTsulasnellaのアノテーションファイルからtsvファイルを作り、発現変動を視覚化・機械学習させて分割させる。(望み薄)

こんな感じ。だんだん何のためにpythonを使っているのかわからなくなってきた。もともとのイメージとしてはPythonで配列データを一つずつ読み込んで、その類似性を見つけてもらってから、主成分分析とかSVMで解析できる、という感じだったけども、

そもそも配列データ量が多すぎてPythonじゃそんなことできないのか。

発現データ解析はRの方が適しているという話だし、Rやるか検討します。