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ツールに頼りたい。。。
論文をまた探してみるか。
