samtoolsの使い方色々
samtoolsの使い方をまとめます
#インストール
基本的にはcondaを使ってインストールするのが良いです。
condaはOSを考えなくていいし、何よりウェブサイトに飛ばなくてもいいうえにコンパイルしなくて良い。
コンパイル時に必要なツールがイモ蔓式に出てきて大変。。。ということもない。
condaをインストールしていない人は是非この機会にanacondaを使ってみてください。
Anaconda | The World's Most Popular Data Science Platform
こっからインストールページに飛べます。
conda install -c bioconda samtools
でインストール完了。
linuxとかWSLとかだとたまにsamtoolsが使えませんと出てくるので、その時は
conda create -n samtools_env -c bioconda
で仮想環境を作ってから
conda activate samtools_env
で仮想環境に入りましょう。
#使い方
samtoolsは基本的にsam, bam, fastaファイルに対応しています。
例えばゲノムのfastaファイルから一部分を切り取りたいときは
samtools faidx "chr?:?-?" genome.fasta > gene.fasta
というようにすれば切り取れます。
SAMからBAMへの変換は
samtools view -bS xx.sam > xx.bam
とすれば完了です。
BAMのソートは
samtools sort xx.bam > xx_sorted.bam
とすればバッチリです。
IGVでマッピング結果を見たいときには
samtools index xx_sorted.bam
としてindexファイルを作ってからIGVを起動するとよいです。
samtoolsはマッピングの状態（ちゃんとマップされてるか、逆位になってないかなど）に従って一部のマッピングのみを切り取ってくれます。
Explain SAM Flags (broadinstitute.github.io)
に従って、samtools -fの後にflagを立てましょう。
例えばread paired かつread reverse strandの場合は1 + 16で17のflagを立てればその条件で抽出してくれます。
samtools view -f 17 xxx.sam -o output.sam
上のリンクに飛んで目的の条件にチェックを入れてflagを表示するのが一番確実でしょう。
RNA-seqアセンブリについて
年度初めということでまた書き始めようと思います。いつまで続くのでしょうか。
さて、これまでひたすらゲノムが読まれていないニッチな生物でRNA-seq解析を行っていました。
流れとしては図の通り。多分ゲノムが読まれていない生物でのRNA-seq解析で一番困るのが図の④の部分のde novo assemblyだと思います。ここでゲノム読まれてないからって言って自分で仮想遺伝子（コンティグ）を作っていくのですが、これがなかなか難しいです。
いろんなツールが混在してる上、パラメーターも色々設定しないといけないしで、結局どうすればいいのか、何を使えばいいのかというのが結構困ると思います。
アセンブリツールについて、一番のオススメはやっぱりTrinity。かなり高精度かつ冗長性の少ないコンティグを作ってくれます。おそらく一番のスタンダードのため、アップデートは頻繁にしてくれるし、オプションの説明も丁寧でわかりやすく、いろんな環境で動きます。
ただ問題はめちゃめちゃメモリ食います。最低でも50GB以上のメモリがないとほとんど動かないんじゃないでしょうか。あとセットアップがめんどくさいです。公式HPに書いてありますが、trinityrnaseq/trinityrnaseq: Trinity RNA-Seq de novo transcriptome assembly (github.com)
jellyfishとか普段あんまり使わないツールも自分で入れなきゃいけません。condaで入れると関連ツール全部入れてくれますが、バージョン指定しても最新版は入んなかった気がします。（入ってv2.4くらいまで）
なので遺伝研スパコンのsingularity経由とかで使うのが一番いいじゃないんでしょうか。
続いてrnaspades。
こちらは微生物ゲノムアセンブリのために作られたspadesをRNA-seq解析用にちょっといじったものらしいです。体感では（Trinityと比べると）かなり高速なうえに正確性もいい気がします。Trinityほどの信頼性は薄い？（感覚ですが）のでevidentialGeneとかcd-hit-estとかでアセンブリをマージすると丁度いいかもしれません。環境もいろんなところで動く印象があるので割と気に入ってます。
次にidba_tran。これは開発されたのがちょっと古くて、最近アップデートもあんまりないしで廃れてしまうのでは、と危惧していますが一番のお気に入りです。
まず高速かつ省エネ。これは計算資源が限られているバイオ系の研究者には一番重要です。Trinityも見習ってほしい。。。どんな環境でも動くうえに、正確性もそこそこなので是非使ってみて欲しいです。
loneknightpy/idba (github.com)
次にtrans-abyss。
Trans-ABySS | Genome Sciences Centre (bcgsc.ca)
これは超有名ゲノムアセンブリ用ツールのAbyssのRNA-seqバージョン。これはリファレンスの作成には正直あまりお勧めはしません。。。
冗長性を諦めてその分考え得る全てのコンティグを排出する、という印象なので（多分SOAP-denovo Transもそうだと思います）出来上がってくるコンティグの数がめちゃめちゃ増えます。
これは決して悪いことではなく、ゲノムがわかっている生物のタンパク質データベースの作成とかスプライシングのされかたとか見るのにはめちゃめちゃいいです。あと異なるk-merで作った複数のアセンブルをマージしてくれる機能がついているので、k-merによる影響を減らしたいときにはめちゃめちゃ便利。活躍する場面は多いとは思います。
最後にvelvet/Oases
これは正直どっちつかずという感じです。よくもないし、悪くもないし、ただちょっと古いのかなあ？というくらいでしょうか。
Velvet and Oases (transcriptome) - Bioinformatics in BioMed (google.com)
DNA-seqで昔メジャーだったツールなようで、今はあえて使う必要はないかもしれません。
以上、五つのRNA-seqアセンブリツールを紹介しましたが、最近では一つのアセンブリツールではなくていくつかのアセンブラーからの結果を統合（マージ）するツールも出てきてます。
オススメはEvidentialGene、Oyster River Protocol、Consembleの三つです。特にOyster River Protocolはセットアップから全て勝手にやってくれる上に、最近ではnatureの論文でも使われてきており、今結構熱いツールです。
次回紹介しようと思います。
久保緑本 レビュー
久保緑本読み終わりましたー
日曜日の集中力が地に堕ちていたので結局月曜日までかかってしまった。
まあ自分の今やっている分野にはあまり必要はないかなー、という感じだった。Rの勉強がてらとか言ってたけどRは別の本で勉強してから読んでね、とのことでしたし。
ただ今後必ず必要になりそうな予感。
内容としてはデータ解析の実践が中心となっている。統計学の基礎を踏まえたうえで、じゃあ実際に解析するにはそれらをどのように扱えばいいのか、現実世界に即したモデルをうまく作るにはどのような考え方をしてどのようなツールを使えばいいのか、とのこと。
まずモデリングの基本になる一般化線形モデルGMLを紹介して、そのパラメーターの決め方として最大対数尤度ではなくAICを使うとover fitting的なのを避けられるよ、てことを解説。
次に、でも現実世界では観測されてない個体差、陽の当たり方とか、気温のわずかな違いとか、そういうのがあるよね、それを考慮するにはどうすればいいの、で出てきたのが一般化線形混合モデル。
これで観測されてない個体差を決めることができる。でも計算量すごいし解析的に解くにはちょっと厳しすぎるよね、ということでMCMCサンプリング法の登場。こいつを使って得られたサンプルとベイズ統計モデルを組み合わせれば、複数パラメーターの尤度推定を簡単にできるじゃん。
こいつをGLM組み込めないかな、GLMをベイズ化して考えればめちゃめちゃ複雑なパラメーターの組み合わせでも尤度を決定できるじゃん、と話は進み、この応用系である階層ベイズモデルや空間構造への組み込みで話は終わり。
ムツカシイ理論の部分をうまくかわしながら、（詳しくは文献とかツールのマニュアルを見てね、と）基本の考え方と適用法をさらっと説明してくれてるので、とても読みやすいしわかりやすい。統計学の基礎を修めて、一度自分でデータをとった人がこれから初めての解析を始めるぞお、というときに必須の本だと思います。
一度目を通しただけですが何らかに応用できそうであればまた読みたいですね。
一回コロナとかKaggleにあるデータ引っ張ってきてPythonなりＲなりでデータ解析を自己満でやりたいと思ってるのでまたお世話になるかもですねー
選ばれるのには理由がある
お久しぶりです。最近は応用情報技術者試験の対策だったり細胞の分子生物学（鈍器）と格闘したりしてました。
ダメもとで応募してた基生研の情報解析講習会の選考通ったので、NGS解析もまた本腰入れてやっていこうと思います。
ということでRを学んでおります。Rを学んで思ったこと。
いやあ、Rって便利！
Pythonじゃ必要なパッケージを呼び出して、テキストファイルを読み込んで、またさらに別のパッケージを呼び出して、matplotlibでもいちいちxとyを指定して、メソッド使って要素追加して、、、とするのをRじゃ２、３行で終わってしまう。
あとRのデータ構造が表形式（行列）に適しているからわざわざNumpyとか使わんでも操作が直観的にわかる。
バイオ系でRばっかり使われているのにはちゃんと理由があるんですねー、という感じ。plotで一発で図形描いてくれるなんてありがたすぎるんですよね。
ということでＲの使いやすさに感動しているこの頃です。もちろんPythonは機械学習、深層学習が非常に得意な言語なので、データ解析の最先端では使っていくことにはなると思うけど。scikit-learnとかね。でもpythonはやっぱりcsvとかテキスト形式にあんまり適応してない。ＲはBioconductorに代表されるように、そういったデータ処理する人々向けにわかりやすく作られているので、そりゃＲ使った方が効率いいよな。
Ｒの勉強は基生研のGithubとか久保緑本：生態学データ解析 - 本/データ解析のための統計モデリング入門 (kuboweb.github.io)を使ってやってます
やっぱりLinux
自分で機械学習させることもないのでしばらくPythonは趣味にとどめておくことにする。ただ機械学習自体は、さらに深層学習自体は何らかに使いたいという願望はあるので、継続はする。
バイオ系でPython使えるイメージがまだまだないので、今後新しいツールが出てきたらそのソースコードから見てみて、イメージをつかみたいところ。
調べたらPerlとRが使えそう。Perlは使用者が比較的少ないし、まだ未発達ということもあり、PythonのNumpyの方がPerlよりもいいですよ、という声もあるので微妙なところだが。でもBioPerlの有用性も、特に配列を用いた計算においては有用ということらしいので、んまあ時間があったら取り組みたい。
最も優先順位が高いのはやっぱRです。Pythonと比べて実装してるバイオ系のツールがケタ違いに多い。Pythonと同じオブジェクト指向型言語ということなので、勉強し始めにそれほど苦労しないだろうということを期待しつつ、こちらもやっていきたい。
で、話は変わり今の問題に対処するために何をすべきかの方向性が定まってきた。
前回挙げた選択肢の中で、特にメタゲノム解析ツール、blobtools（をはじめとする種々のツール）とマッピング・アラインメントを駆使したコンティグ整理、の二つをやっていこうと思う。
blobtoolsのコンセプトは、「シーケンスの際のContaminant(他の生物のゲノム)をスクリーニングする」というもの。実際結構他の生物種のコンタミは多いらしく、blobtoolsの論文の中には他の生物種の遺伝子が混じってただけなのにろくにスクリーニングしないで水平伝播だと勘違いされていた例も多いらしい。
こんなのがまかり通ってたらデータベース自体も汚染されるし、もっと注目されるべき問題なんだろうなあ。
blobtoolsの論文内にあった使えそうなツールをいくつか列挙する。
Anvi’O、t-SNE、ProDeGe、あとQC-Blind。PubMedで調べればもっと出てきそうだけどひとまずこんなもんで。さらっと見た感じQC-Blindがリファレンスゲノムもいらないらしくて一番期待。
これで何ともなんなかったらメタゲノム解析ツール使おう。
メタゲノムデータのtaxonomy assignmentを行う k-SLAM - macでインフォマティクス (hatenablog.com)
この記事参照かな。
最後にblobtoolsのもと論文を読んだから少し書き残しとく。
必要なファイルはnBlastとDiamond Blastxで作ったtsvファイル（シーケンスに相同配列検索が掛けられたもの）、アッセンブルした後のFASTAファイル、元リードをアッセンブルした後のBAMファイルである。
こいつらを準備したら下のworkflawA（リファレンスゲノムがないもの）に従って処理していく。
もとファイルを集めたblobDBファイルを作って、それをpartition及びplotしてもらい、次にBAMファイルとこのpartitionの結果を合わせてリードの分割を行って終了。
この処理は二回以上行うのが推奨ぽい？
メタゲノムcontigのカバレッジ、GC、taxonomy情報を可視化して分析できる BlobTools - macでインフォマティクス (hatenablog.com)
に実行したときのこととかエラーのことが詳しく書かれているので実際にやるときにはこれを参照されたい。
RNA発現解析
次にマッピングしたデータを処理する。
参照配列の塩基にどれだけのリードがあったかを集計・整理する。リードを集計・整理できたら（リードカウント数として確認出来たら）そのカウント数をTPMなどで補正してそれぞれの出現量とする。
さらに出現量の違いを距離として入力したサンプル間で系統樹を書いたりなどの処理が可能となる。
また発現量の違いが大きい遺伝子を抽出すればそれらの生物的な関連をジーンオントロジーやパスウェイに当てはめて考察できる。
まず必要なことは参照配列のgffデータに対して遺伝子ごとに通番のgene_idを挿入する。これはのちの作業のインデックス化に必要となる元データになる。
pip install bcbio-gffでbcbio-gffをインストールしたのちに元のgffファイルから作る。
with open('gffファイル名',"w") as (name):
でgffファイル名を開き、GFF.parseメソッドを用いてfor構文で一行ずつ読み出す。
そのさいにfeaturesをfor構文で読みだしてfeature.typeからgeneとpseudogeneに一致するものの、gene_idを0001から連番でつけていく。f.qualifiersで名づけたら次にsub_featuresも同様に、さらにそのsub_featuresのgene_idも同様に名付けていく。
最後にGFF.writeメソッドで出力して終了。
これでgene_idを連番にした参照ファイルを作ることができた。
次にこれからリードカウント数を数えるためにsbreadパッケージをダウンロードしてfeatureCountsを使って作業を行った。
The Subread package (sourceforge.net)
のLatest Version2.0.1のリンクからダウンロードした。
ダウンロード後はsudo apt install tarで
tar -xvcf
でtar.gzファイルを解凍して、
目的ファイルのあるディレクトリまで下りたらmake -f Make.Linuxと打ち込んで完了。あとはsudo apt install subreadでfeatureCountsが使える状態となった。
featureCountsはオプションが色々とあり、-pでペアエンドの指定、-T （数字）でスレッド数の指定、-t exonでexonを数えることを指定、-g （アトリビュート）でグループ化の際のアトリビュートを指定できる。また-a （アノテーションファイル）でアノテーションファイルを指定する。これでアトリビュートごとのカウント数を得ることができる。
使い方の詳細は公式サイト
http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf
を参照のこと。
一通り処理が終わったらリードカウント数のデータがまとめて入ったファイルがtsv形式で出てきた。RNA発現量解析はtsvデータをいじって行うらしい。
とまあこんな感じでいったんデータの集計は終わり。
次からこのtsvファイルを使ってTPM値を使ってデータの補正を行う。
あと今参考にしている本がさらに参考にしてたサイト見つけたのでのっけとく。
これさえ見ればpython使ったNGS解析は一応できるんじゃ、てほど丁寧。
動画付き。
情報解析講習会ビデオ＜2019年度第1回PAGS・DDBJ・DBCLS合同情報解析講習会＞ ｜ 先進ゲノム支援（先進ゲノム解析研究推進プラットフォーム） (genome-sci.jp)
ここで話は変わるが宿主と寄生生物の遺伝子の分割、SVMを用いたwebツールを使えば一発ととある論文に書いてあったので楽観視してこれまで調べなかったのだが、今日そのサイトをみたらサービス終了してた。。。
そもそもその論文自体が2005年に発行されたものだから仕方がないっちゃないけども、また論文探しから始めないといかん。
一応使えそうなツールとして、まだしっかりと論文を読んではいないがblobtoolsなるものが(運がいいことに)pythonのツールとして2018年にGithubで出されてた。
これ使えばなんとかなるかもしれないが、さすがに機械学習をさせたいので、Webツールに頼りたい。。。
論文をまた探してみるか。
ubuntuでAnaconda3を入れた話
はい、とうとうubuntuでAnaconda3を入れました。
別にubuntuではAnacondaの縮小バージョンminicondaでいいかな、と思ってたけど、minicondaからだとどうもpythonがバージョン2しか稼働しないぽく。。。
するとpython3以降からのみ対応しているパッケージが全部minicondaによって邪魔されてインストールできないor稼働しない。。。
新しいバージョンのminicondaをインストールする手もあったんだろうけど、linux置いてるところからjupyternotebook使えるのも便利だし、そこもう統一しようかということで入れちゃいました。
割と四苦八苦したので書いておく。
まずネットに公式サイトからLinux用のAnacondaをインストールしてくださいって書いてたんだけど、そのままインストールするだけじゃ多分無理。
ubuntuの置いてあるファイルがWIndowsの環境から独立しているので、インストールしたところでubuntuのとこまで正常に持っていけない。
なのでubuntuのコマンドから
wget (ダウンロードのurl)
でダウンロードするのが無難。
無事ダウンロードできたら
bash Anaconda3-2020.11-Linux-x86_64.sh
でダウンロードしたAnacondaファイルを開く。このときの2020.11は今後Anacondaのバージョンアップに伴い変化するので、そこは適宜変更する。
ファイルを開いたらあとはubuntuの画面に出てくる通りにエンターとかyesとか押してけばインストールは完了。
でも大体はここでパスが通ってない。
Thank you for installing
みたいな文字が出て打ち込み画面が現れたら
conda -V
でcondaにパス通っているか確認。パス通ってなかったらエラーが出てくる。
なのでubuntuのユーザーディレクトリにある.bashrcのファイルを開いて自分で直接パスを書き込む。ubuntuなので直接書かなきゃいけないのがめんどくさいね。WIndowsなら環境変数でいけるのに。ちなみに.bashcとかそういうファイルを開くためにAtomをインストールしとくとよい。このアプリ見やすくてお勧め。
.bashrcの最後の行にexport PATH=~/anaconda3/bin:$PATH
と打ち込んでsaveしたらubuntu に戻って
source ~/.bashrc
と入力してbashrcファイルを読み込んで終了。
これでパスも通ったので
conda -V
と打てばバージョン情報が出てくるはず。
あとubuntuからjupyter notebookも使えて、
jupyter notebook --no-browser
と打てばリンクを返してくれるからそれをweb検索エンジンに打ち込めばhomeディレクトリ内のフォルダが既に入ってる状態のjupyternotebookを返してくれる。優秀。
anacondaが無事入ったのですかさずpip install biopython、pip install bcbio-gff、conda install hisat2でインストール。無事入りましたとさ。
こういうもので壁にぶち当たったり沼にハマったときのひとつの解決策として、確かにそのエラーなりの解決策を考えるのも重要だけど、そんなに時間もないわけで。
もう沼から抜け出せねえ、てなったらいっそそのパッケージなりごと消してもう一回インストールしなおす、というのも重要かもしれん。根本からぶっ壊す、みたいな。
