============== -*- outline -*- ============== bioinfo_tools 2021/03/22 Kuninori Morimoto <kuninori.morimoto.gx@renesas.com> ============================================= * 何これ？ 某研究者から bioinformatics 関連のツールを作って欲しいと 依頼を受けて作っているやつ * ライセンス GNU General Public License v2 * MAC ユーザー MAC の sed (BSD sed) の改行の扱いが相当ややこしいので MAC ユーザーは GNU sed を使ってください > brew install gnu-sed * mo_dfast_to_cdspickup.sh https://dfast.ddbj.nig.ac.jp/ http://kazumaxneo.hatenablog.com/entry/2019/10/10/073000 https://github.com/nigyta/dfast_core 上記の Dfast は A) RARx-METABAT__x/protein.faa B) RARx-METABAT__x/genome.gbk A) に説明を添えた B) を作成する mo_dfast_to_cdspickup は B) の CDS の /note の中に特定のキーワードを含んでいるものを A) からピックアップするために、中間ファイルである mo_cdspickup_xxx.xls を作る (xxx は指定したキーワード)。 以下は "AmoA" と言うキーワードを全フォルダ (= RAR*) からピックアップする例 > ls RAR1-METABAT__1 RAR2-METABAT__2 RAR3-METABAT__3 ... > mo_dfast_to_cdspickup.sh AmoA RAR* 上記例の場合、実行後に mo_cdspickup_AmoA.xls と言うファイルができる。 確認のためこのファイルを MS Office 等で開き余計なものが含まれていないか等を確認する。 余計なものが含まれている場合は、その行を削除する。 確認・編集後、「CSV 形式」または「テキスト（タブ）形式」で保存する。 保存したファイルを使い、次の mo_cdspickup_to_faa.sh を使用する事で、 特定のキーワード含んだ新しい protein.faa を作成する * mo_cdspickup_to_faa mo_dfast_to_cdspickup により作成された mo_cdspickup_xx ファイルを元に 新しい protein.faa を作成する 以下は "AmoA" と言うキーワードを全フォルダ (= RAR*) からピックアップした mo_cdspickup_AmoA.xls （テキスト（タブ）形式）を元に全フォルダ内にある protein.faa からキーワードを含む情報だけをピックアップする。 > ls RAR1-METABAT__1 RAR2-METABAT__2 RAR3-METABAT__3 ... > mo_cdspickup_to_faa.sh mo_cdspickup_AmoA.xls RAR* 上記の例の実行後、mo_protein.faa が作成される * mo_analysis.sh ------------------- MGA_2|NED1_00010 MGA_3|NED1_00020 K06925 MGA_4|NED1_00030 K00259 ... MGA_11|NED2_00110 MGA_12|NED2_00120 K04065 MGA_13|NED2_00130 K01480 ... MGA_12|NED3_00110 K01284 MGA_13|NED3_00120 MGA_14|NED3_00130 ... ------------------- 上記のようなファイルから、Core, Singleton 等の情報を解析する > mo_analysis.sh sample.txt summary Total genes (a)全遺伝子数（入力ファイルの行数） Genes with K No (b)全遺伝子に付与された K No の数 rate (genes) (b) * 100 / (a) Uniq K No 全 K No をユニークにした数 Total Core (x)全コアの数 rate (genes) (x) * 100 / (a) rate (K No) (x) * 100 / (b) Total Singleton (y)全シングルトンの数 rate (genes) (y) * 100 / (a) rate (K No) (y) * 100 / (b) 各ゲノムフォルダ all_genes.tsv 各ゲノムの全遺伝し数 all_k_no.tsv 各ゲノムに付与された K No の数 uniq_k.tsv all_k_no の内容をユニークにしたもの 出現個数付き singleton.tsv Singleton の Kxxx summary.tsv Total genes (1)全ゲノム数 (all_genes.tsv の行数) Genes with K No (2)Kxxxx の数 (all_k_no.tsv の行数) rate (genes) (2) * 100 / (1) Uniq K No Kxxx をユニークにした数 Singleton (3)シングルトンの数 rate (genes) (3) * 100 / (1) rate (K No) (3) * 100 / (2) core フォルダ uniq_core.tsv 全コア 出現個数付き summary.tsv Total core (x)全コア数（uniq_core.tsv の行数） SCG (9/9) (y)9 ゲノム中 9 個がシングルの数 rate (9/9) (y) * 100 / (x) SCG (8/9) (z)9 ゲノム中 8 個がシングルの数 rate (8/9) (z) * 100 / (x) ... * mo_kno_cnt ------------------- MGA_2|NED1_00010 MGA_3|NED1_00020 K06925 MGA_4|NED1_00030 K00259 ... MGA_11|NED2_00110 MGA_12|NED2_00120 K04065 MGA_13|NED2_00130 K01480 ... MGA_12|NED3_00110 K01284 MGA_13|NED3_00120 MGA_14|NED3_00130 ... ------------------- 上記のようなファイルで、指定された K No が各ゲノムに何回 登場するか数える > ./mo_kno_cnt.sh K21572 user_ko.txt NED1 5 NED2 5 NED3 4 NED4 4 NED5 7 NED6 7 NED7 7 NED8 10 NED9 1 * mo_kno_name 指定された K no の名前を取得する カテゴリ、シンボル、名前をタブ区切りで表示 複数のカテゴリがある場合は複数行で表記 カテゴリ名などは KEGG ファイルから取得する KEGG ファイルは１日１回更新する (以下はわかりやすくするためにタブを揃えているが、本当はタブは１つ) > mo_kno_name.sh K00001 09100 Metabolism 09101 Carbohydrate metabolism 00010 Glycolysis / Gluconeogenesis [PATH:ko00010] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09101 Carbohydrate metabolism 00620 Pyruvate metabolism [PATH:ko00620] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09103 Lipid metabolism 00071 Fatty acid degradation [PATH:ko00071] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09105 Amino acid metabolism 00350 Tyrosine metabolism [PATH:ko00350] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09108 Metabolism of cofactors and vitamins 00830 Retinol metabolism [PATH:ko00830] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09111 Xenobiotics biodegradation and metabolism 00625 Chloroalkane and chloroalkene degradation [PATH:ko00625] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09111 Xenobiotics biodegradation and metabolism 00626 Naphthalene degradation [PATH:ko00626] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09111 Xenobiotics biodegradation and metabolism 00980 Metabolism of xenobiotics by cytochrome P450 [PATH:ko00980] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] 09100 Metabolism 09111 Xenobiotics biodegradation and metabolism 00982 Drug metabolism - cytochrome P450 [PATH:ko00982] E1.1.1.1, adh alcohol dehydrogenase [EC:1.1.1.1] * mo_genome_function.sh File style ---------- MGA_3|NED1_00020 K06925 ... ---------- or ---------- YTPLAS18_00010 K02313 ... ---------- 指定されたファイルを元に K ナンバーのカテゴリや名前、 EC 経路を表記し、各ゲノムがその K ナンバーを持っているかどうかの一覧を作る ${FILE}_genome_function.tsv と言うファイルに出力する 表計算ソフトで開ける > mo_genome_function.sh ./user_ko.txt user_ko_genome_function.tsv ${FILE}_genome_function.xls と言うファイルも出力するが 中身は実は HTML。開く時何か言われるけど無視して開ける。 EC や K ナンバーがリンクになっててクリックすると KEGG の Web ページが開ける 重複する K ナンバーがある場合、-1 オプションをつけると最初のものだけを表示する 重複する K ナンバーすべて表示 > mo_genome_function.sh ./user_ko.txt 重複する K ナンバーの最初の１つだけを表示 > mo_genome_function.sh -1 ./user_ko.txt [FIXME] KEGG の Web ページは EC 番号をふって無い K ナンバーも EC として認識するけど、KEGG ファイル自体はふって無い このスクリプトはそれを考慮できない 時間が経てば解決する問題かも？ * mo_ec_to_kno.sh EC 番号を入れると、必要な K ナンバーの一覧を返す > mo_ec_to_kno 2.2.1.2 K00616 K13810 * mo_ec_pass.sh 指定された EC を通れるかどうかを示す. OK/NG で表記した後、根拠も示す オプションで指定する場合 > mo_ec_passh.sh user_ko.txt 2.2.1.2 OK : 2.2.1.2 ... 複数の EC を経由する場合は "-" でつなげる > mo_ec_passh.sh user_ko.txt 4.4.1.22-1.1.1.284-3.1.2.12 NG : 4.4.1.22-1.1.1.284-3.1.2.12 ... スペースで区切る事で、複数の経路を同時に指定もできる > mo_ec_passh.sh user_ko.txt 2.2.1.2 4.4.1.22-1.1.1.284-3.1.2.12 分かりやすくファイルに落としておいて、-f で指定もできる --- Formaldehyde_to_Formate --- 1.2.1.46 4.4.1.22-1.1.1.284-3.1.2.12 ------------------------------- > mo_ec_passh.sh user_ko.txt -f Formaldehyde_to_Formate NG : 1.2.1.46 NG : 4.4.1.22-1.1.1.284-3.1.2.12 .... KEGG ファイルには EC が付いてないものの、KEGG の Web ページでは （なぜか）EC に含まれる K ナンバーも対象に入れたい場合、 mo_ec_pass.conf ファイルを作成すると追加できる ---- mo_ec_pass.conf ---- 1.14.13.25:K16160,K16162 ... ------------------------- この場合、結果には(custom) と表示される > mo_ec_passh.sh user_ko.txt -f GenomeX_to_GenomeY NG : w.x.y.z OK : a.b.c.d-1.14.13.25 (custom) .... * mo_add_prefix_to_file.sh mo_add_prefix_to_file PREFIX FILES PREFIX: 頭につけるプレフィックス FILES: プレフィックスをつけるファイル群 > ls Bin_Bin_0.faa Bin_MAXBIN__017_sub_1.faa Bin_METABAT__100.faa ... Bin_Bin_24_sub_1.faa Bin_MAXBIN__017_sub_2.faa Bin_METABAT__101.faa ... Bin_Bin_41_sub.faa Bin_MAXBIN__018_2.faa Bin_METABAT__102.faa ... ... 指定されたファイルにプレフィックスと連番をつける 上記のファイル群の頭に「LAS_xxx-」と言うプレフィックスをつける場合 > mo_add_prefix_to_file.sh LAS ./*.faa > ls LAS_000-Bin_Bin_0.faa LAS_000-Bin_MAXBIN__017_sub_1.faa ... LAS_0xx-Bin_Bin_24_sub_1.faa LAS_0xx-Bin_MAXBIN__017_sub_2.faa ... LAS_0xx-Bin_Bin_41_sub.faa LAS_0xx-Bin_MAXBIN__018_2.faa ... ... * mo_stream_dfast.sh 指定されたファイルを順番に dfast にかける この時、ファイルのプレフィックスを --locus_tag_prefix に使用する mo_add_prefix_to_file.sh を使用している事を想定 > ls LAS_000-Bin_Bin_0.faa LAS_000-Bin_MAXBIN__017_sub_1.faa ... LAS_0xx-Bin_Bin_24_sub_1.faa LAS_0xx-Bin_MAXBIN__017_sub_2.faa ... LAS_0xx-Bin_Bin_41_sub.faa LAS_0xx-Bin_MAXBIN__018_2.faa ... > mo_stream_dfast.sh ./*.faa LAS_000-Bin_Bin_0.faa ^^^^^^^ これが使用される * mo_dfast_statistics_rna_table.sh DFAST の出力 *.fa の「statistics.txt」「rna.fna」を表にする *.fa のある場所で実行すると dfast_statistics_rna_table.xls に出力する > cd DFAST_OUT > ls aaa.fa bbb.fa ccc.fa ... > mo_dfast_statistics_rna_table.sh > ls dfast_statistics_rna_table.xls