104ページkallistoの出力ファイルを読み込む、について

[Dsaito-gunma]

学生です。教本でNGS解析の勉強をしています。

kallistoとsleuthの処理を行っているのですが、 104ページのkallistoの出力ファイルを読み込む、のコマンドを入力したところ、

In check_target_mapping(tmp_names, target_mapping) : intersection between target_id from kallisto runs and the target_mapping is empty. attempted to fix problem by removing .N from target_id, then merging back into target_mapping. please check obj$target_mapping to ensure this new mapping is correct.

と出てしまい、先に進めません。 sleuth_prep() Error #111ものぞいたのですが、解決しませんでした。 エラーの意味が分からず、困っています。 R初心者の質問ですみませんが、よろしくお願いいたします。

[Dsaito-gunma]

@RolyPolyCoily お世話になります。質問してから1週間になりますが、いかがでしょうか。お忙しいところ恐縮ですが、よろしく回答お願い申し上げます。

[RolyPolyCoily]

お返事しておらず申し訳ございません。 kallistoの出力とtarget2gene.txtの互換性がないことに原因があるかもしれません。 どれか適当なサンプルについて下記コマンドの出力、 head abundance.tsv さらに下記の出力 head target2gene.txt を教えていただけますでしょうか。 また、たとえば下記コマンドではどのようにメッセージが出るかもご確認いただけますと幸いです。 so <- sleuth_prep(s2c,extra_bootstrap_summary=T) お手数おかけしますがよろしくお願い申し上げます。

[Dsaito-gunma]

@RolyPolyCoily お忙しいところお返事いただき、ありがとうございます。 明日、出力結果を報告いたします。今後ともどうぞよろしくお願いいたします。

[Dsaito-gunma]

@RolyPolyCoily いつもありがとうございます。先日教えていただいた出力ですが、 SRR1550989のサンプルのhead abuncance.tsvは target_id length eff_length est_counts tpm ENST00000631435.1 12 8 0 0 ENST00000415118.1 8 8 0 0 ENST00000448914.1 13 6.33333 0 0 ENST00000434970.2 9 5 0 0 ENST00000632684.1 12 8 0 0 ENST00000633010.1 16 6.8 0 0 ENST00000633009.1 20 7.44444 0 0 ENST00000632524.1 11 7 0 0 ENST00000633353.1 31 8.30769 0 0

head targrt2gene.txtは target_id ens_gene ext_gene ENST00000387314 ENSG00000210049 MT-TF ENST00000389680 ENSG00000211459 MT-RNR1 ENST00000387342 ENSG00000210077 MT-TV ENST00000387347 ENSG00000210082 MT-RNR2 ENST00000386347 ENSG00000209082 MT-TL1 ENST00000361390 ENSG00000198888 MT-ND1 ENST00000387365 ENSG00000210100 MT-TI ENST00000387372 ENSG00000210107 MT-TQ ENST00000387377 ENSG00000210112 MT-TM

また、最後のコマンドですが、

so <- sleuth_prep(s2c,extra_bootstrap_summary=T) reading in kallisto results dropping unused factor levels ....... normalizing est_counts 64663 targets passed the filter normalizing tpm merging in metadata summarizing bootstraps ....... ＞ となっております。不足な情報がございましたらお申し付けください。 お手数おかけしてばかりで申し訳ございませんが、よろしくお願いいたします。

[Dsaito-gunma]

@RolyPolyCoily 追加情報です。 先程のエラーメッセージはUbuntu環境で行ったものでしたので、一応、Macでもやってみました。 下記のように、同様のエラーが出ますので、OSの違いによるものではなさそうです。

警告メッセージ: 1: check_target_mapping(tmp_names, target_mapping, !is.null(aggregation_column)) で: intersection between target_id from kallisto runs and the target_mapping is empty. attempted to fix problem by removing .N from target_id, then merging back into target_mapping. please check obj$targetmapping to ensure this new mapping is correct. 2: `select()was deprecated in dplyr 0.7.0. Please useselect()instead. This warning is displayed once every 8 hours. Calllifecycle::last_warnings()` to see where this warning was generated.

ちなみに、lifecycle::last_warnings()を入力してみると、下記メッセージが出ます。

message: `select_()` was deprecated in dplyr 0.7.0. Please use `select()` instead. Backtrace: 1. sleuth::sleuth_prep(...) 2. dplyr::select_(target_mapping, "target_id", aggregation_column) 3. dplyr:::lazy_deprec("select", hint = FALSE)

[RolyPolyCoily]

@Dsaito-gunma 最初のご質問にある「先に進めない」という部分ですが、それ以降のコマンドが実行不可になりますでしょうか。 よく見ると本件はエラーではなく警告のようです。abuncance.tsvにあるtarget_idにはピリオドに続いて１桁の数字がついていますが、target2gene.txtにはありません。この部分を .N と形容して、それを除外してマッチングさせた、という内容のようです。もしコマンド実行可能でしたら、この警告は無視して操作を進め、書籍と同様の解析結果が得られるかどうかをご確認いただけますでしょうか。 もしコマンド実行不可になるようでしたら、どのようなメッセージが表示されるか共有いただければ幸いです。

[Dsaito-gunma]

@RolyPolyCoily ご連絡ありがとうございます。 Ubuntuでやった時は、ファイルが見つからない、というエラーが出ていました。 このため、Macで行い、先に進めてみました。(この機会にMacに切り替えます) 私が行ったところ、アイソフォームレベルではENST00000116717(GADD45A)、ゲノムレベルではENSG00000285395(XYLT1)が、p値が最小と出ました。このため、書籍とは異なる結果かと思います。 また、p値が低い順に並べ替えた上で作成したヒートマップが添付のようになっており、p値の並べ替えが充分に反映されていない図になっています。コマンドは書籍およびGithubの通りに入力しました。可能でしたら、この点もコメントいただけませんでしょうか。

heatmap.pdf

[RolyPolyCoily]

@Dsaito-gunma お返事できておらず大変申し訳ございません。こちら原因等を調べるのに時間がかかっております。お時間いただけますと幸いです。 また、p値の並べ替えが十分にされていない点ですが、GADD45AのENST番号はENST00000370986.9である可能性はございますでしょうか。お示しいただいた番号はENSGの可能性がございます。

[Dsaito-gunma]

@RolyPolyCoily ご連絡ありがとうございます。また、お手数かけまして申し訳ございません。GADD45AのENST番号はENST00000370986.9であり、ENSG00000116717でした。ご迷惑をおかけし、大変失礼しました。

[RolyPolyCoily]

@Dsaito-gunma おまたせしており申し訳ありません。同じ現象を再現できず検証できておりませんが、もとのリファレンスデータに原因があるのではないかと考えております。お手数ですが、書籍90ページの４つ目の枠にあるコマンドで...cdna.all.fa.gz を再取得していただきたいです。 また同じ場所に CHECKSUMS というファイルもあるかと思います。こちらも get CHECKSUMS により取得お願いします。 ２つのファイルをダウンロードしたら、Terminalで以下のコマンドを実行してください。 sum Homo_sapiens.GRCh38.cdna.all.fa.gz 返される数値が、CHECKSUMS（catなどで表示できます）に記載されているもの（例：48231 66314 Homo_sapiens.GRCh38.cdna.all.fa.gz）と同じであることを確認してください。

その後、kallistoによる再インデックス化〜Sleuthでの発現差解析を実施し、書籍と同じ結果が得られるかご確認いただけますでしょうか。よろしくお願いいたします。

[Dsaito-gunma]

@RolyPolyCoily ご回答いただき、ありがとうございます。また、私の解析の時間がかかってしまい、申し訳ございませんでした。 CHECKSUMSを取得し、数値の確認をしましたが、Homo_sapiens.GRCh38.cdna.all.fa.gzと同じ数値であることを確認しました。 その上で解析したところ、アイソフォームレベルでは書籍と同様の結果で、ヒートマップやvolcanoplotも同じものが作図できます。 しかし、遺伝子レベルでは、ENSG00000059804 SLC2A3のp値が最小となっており、ENSG00000169710 FASNのp値は5番目に小さい値として出てきました。不足な情報がございましたら、お申し付けいただけますと幸いです。

[RolyPolyCoily]

@Dsaito-gunma ご連絡ありがとうございます。アイソフォームレベルでは再現できたのに遺伝子レベルでは再現できないというのは、また別の背景がありそうですね。原因を考えてみますのでしばらくお待ちください。…なお解析はご自身の都合で進めていただければと思います。私のことは気にしないでください。

[RolyPolyCoily]

@Dsaito-gunma だいぶおまたせしてしまい申し訳ございません。こちら私自身の環境、研究室内の別の者の環境等で解析した結果をみたところ、ご指摘のとおりSLC2A3が最小のP値を示すようです。執筆当時の解析環境が動かなくなっており、原因は定かではありませんが、少なくともご連絡いただいた解析結果で間違いないようです。

[Dsaito-gunma]

@RolyPolyCoily ご連絡ありがとうございます。私の環境が問題ないようで、安心しました。しかし、解析環境の変化で結果が変わった、とのことですが、それですと、自分の解析が正しいか、なかなか判断が難しそうですね。その点、どのように考えればよろしいでしょうか。漠然とした質問で申し訳ございません。

[RolyPolyCoily]

@Dsaito-gunma ご心配とお手数おかけしました。今回は本当に解析環境の違いに起因するのか不明ですが、ツールやデータセットのバージョンが違うことで完全に同じ処理結果が得られないことはあります。ただし結果の解釈が変わるほど違う結果を返すことは多くないと思います。たとえば今回もSLC2A3がトップに出ましたが、qvalをみると有意ではありません。したがってgene-levelでは、有意に発現差が見られた遺伝子は無かったという結論は変わりません。 逆に言えば、確実に検出されるであろう既知の変動遺伝子（ポジコンのようなもの）があれば、自身のデータセット＆解析手法でもその遺伝子を検出できるかを確認することで、自分の解析が正しそうか判断する材料になると思います。