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

[nakane-scc]

Linuxにて作業をしています。

下記のように、104ページのkallistoの出力ファイルを読み込むを実行しました。

so <- sleuth_prep(s2c, extra_bootstrap_summary=T, target_mapping=t2g) すると、次のようなエラーが出てしまい、処理がされませんでした。 . get_kallisto_path(path) でエラー: 'SRR1550989' does not exist.

作業ディレクトリは、kallistoにしており、s2cを表示させると、きちんとデータはあるようです。 sample condition path SRR1550989 Healthy_Control SRR1550989 SRR1551005 Type_1_Diabetes SRR1551005 SRR1551011 Type_1_Diabetes SRR1551011 SRR1551050 Healthy_Control SRR1551050 SRR1551057 Healthy_Control SRR1551057 SRR1551071 Healthy_Control SRR1551071 SRR1551091 Type_1_Diabetes SRR1551091

どうすれば良いでしょうか。 ご多用のところ恐縮ですが、教えていただければ幸いです。

[RolyPolyCoily]

@nakane-scc こちらもお返事が遅くなり、メールでのやり取りとなってしまいました。 fastqファイルを７検体分ダウンロードできていなかったことが原因でしたので、足りない１検体分のデータを追加して再解析を実施するか、sample2conditionのファイルで欠損している検体の行を削除することでエラー回避できるかと思います。 ７検体分ダウンロードできていることを確認するステップを追記したいと思います。

[nakane-scc]

ご回答有り難うございます。

7検体分データをダウンロードして、

so <- sleuth_prep(s2c, extra_bootstrap_summary=T, target_mapping=t2g, aggregation_column='ens_gene', gene_mode=T)

を実行したところ、エラーは出なかったのですが、下記のところで止まってしまいました。

'gene_mode' is TRUE. Sleuth will do counts aggregation at the gene level for downstream normalization, transformation, and modeling steps, as well as for plotting and results. reading in kallisto results dropping unused factor levels ....... normalizing est_counts 62667 targets passed the filter normalizing tpm merging in metadata aggregating by column: ens_gene 16071 genes passed the filter summarizing bootstraps .......

10分近く待ってみたのですが、データを読み込むだけでそんなに時間がかかるのはおかしいのではないかと思い、Ctrl+Cでコマンドを停止させてみました。すると、次のような警告メッセージが出ました。

警告メッセージ: 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$target_mapping to ensure this new mapping is correct. 2: sleuth_prep(s2c, extra_bootstrap_summary = T, target_mapping = t2g, で: 3523 target_ids are missing annotations for the aggregation_column: ens_gene. These target_ids will be dropped from the gene-level analysis. If you did not expect this, check your 'target_mapping' table for missing values.

target2gene.txtに問題があるのでしょうか。 目視した限りでは、91ページの画像の内容と同じに見えました。

なお、

so <- sleuth_prep(s2c, extra_bootstrap_summary=T, target_mapping=t2g)

上記コマンドを実行すると、次のような警告メッセージは出ましたが、一応、読み込みはできたようです。

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$target_mapping to ensure this new mapping is correct.

[RolyPolyCoily]

@nakane-scc ありがとうございます。たしかに10分もかかるのは長いですが、特段エラーメッセージ（警告メッセージは無視して構いません）が出ているわけではないようなのでデータに問題はないかと思います。こちら何度試しても読み込みが進まない状態でしょうか。

[nakane-scc]

返信、ありがとうございます。

何度か試し、１時間前後放置してみたのですが、先には進みませんでした。 今年購入したメモリ128GB、core i9のCPUを積んだパソコンを使っているので、通常の処理であればそんなに時間はかからないはずなのですが。

[nakane-scc]

試しに、sample2condition.txtのサンプル数を減らして実行してみました。

sample condition path SRR1550989 Healthy Control SRR1550989 だけをファイルに記載した状態で、 so <- sleuth_prep(s2c, extra_bootstrap_summary=T, target_mapping=t2g, aggregation_column='ens_gene', gene_mode=T) を実行すると、次のような警告メッセージは出たものの、一応、処理は終わりました。

警告メッセージ: 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_lifecycle_warnings()` to see where this warning was generated. 3: sleuth_prep(s2c, extra_bootstrap_summary = T, target_mapping = t2g, で: 3523 target_ids are missing annotations for the aggregation_column: ens_gene. These target_ids will be dropped from the gene-level analysis. If you did not expect this, check your 'target_mapping' table for missing values. 4: data.table::melt(bs_df, id.vars = "bootstrap_num", variable.name = "target_id", で: The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(bs_df). In the next version, this warning will become an error. 5: data.table::dcast(scaled_bs, sample ~ target_id, value.var = "scaled_reads_per_base") で: The dcast generic in data.table has been passed a data.frame and will attempt to redirect to the reshape2::dcast; please note that reshape2 is deprecated, and this redirection is now deprecated as well. Please do this redirection yourself like reshape2::dcast(scaled_bs). In the next version, this warning will become an error.

次に、サンプルを SRR1551005 Type 1 Diabetes SRR1551005 だけにしてプログラムを実行すると、やはり警告メッセージは出たものの、処理は終わりました。

そこで、今度は処理に成功したSRR1550989,SRR1551005の2つをテキストに記入してプログラムを実行してみると、30分が経過しても処理が終わりませんでした。

原因、対処方法がおわかりになりますでしょうか。また、警告メッセージにはどのように対応すればよろしいでしょうか。

[RolyPolyCoily]

@nakane-scc エラーが出ずに止まったままというのが不思議ですね。もし電源をつけっぱなしにできるようでしたらそのまま一晩放置してみていただけますでしょうか。

[nakane-scc]

お世話になってます。

昨日、こちらに書き込みをしながら６時間くらい放置してみたのですが、止まったままでした。再度、2サンプルだけ記載したsample2condition.txtを用いて、明日まで放置してみます。

[nakane-scc]

24時間、パソコンを放置してみましたが、下記のところで止まったままでした。どうすればよろしいでしょうか。

'gene_mode' is TRUE. Sleuth will do counts aggregation at the gene level for downstream normalization, transformation, and modeling steps, as well as for plotting and results. reading in kallisto results dropping unused factor levels .. normalizing est_counts 77499 targets passed the filter normalizing tpm merging in metadata aggregating by column: ens_gene 17838 genes passed the filter summarizing bootstraps ..

[RolyPolyCoily]

@nakane-scc ありがとうございます。こちらでも summarizing bootstrapsのあとの ....... で少々待ちますが、24時間もかかりません。同じ症状を再現できるかわかりませんが、下記の出力を教えていただけますでしょうか。 sessionInfo() packageVersion("sleuth") packageVersion("tximport") あわせて、kallistoで作成された .h5 拡張子のついたファイルのサイズや行数も教えていただきたいです。

同時並行で、condaの仮想環境を作成、その環境中で動作をみていただきたいです。 conda create -n ngs -y conda activate ngs conda install -c bioconda r-sleuth -y conda install -c bioconda bioconductor-tximport -y cd ~/Documents/expression/kallisto R library(sleuth) library(tximport) tmp=read.csv("../seq/SraRunTable.csv",stringsAsFactors=F) s2c=data.frame(sample=tmp$Run,condition=tmp$diseasestatus,path=tmp$Run) s2c$condition=gsub(" ","_",s2c$condition) s2c$path=file.path(s2c$path) t2g=read.table("../ref/target2gene.txt",header=T,stringsAsFactors=F) so=sleuth_prep(s2c,extra_bootstrap_summary=T) これでも読み込みが止まってしまうかどうか、結果を教えていただきたいです。

[nakane-scc]

お世話になっております。

ご案内いただきました作業環境は以下の通りです。

sessionInfo() R version 4.1.2 (2021-11-01) Platform: x86_64-redhat-linux-gnu (64-bit) Running under: Fedora Linux 35 (KDE Plasma)

Matrix products: default BLAS/LAPACK: /usr/lib64/libflexiblas.so.3.1

locale: [1] LC_CTYPE=ja_JP.UTF-8 LC_NUMERIC=C
[3] LC_TIME=ja_JP.UTF-8 LC_COLLATE=ja_JP.UTF-8
[5] LC_MONETARY=ja_JP.UTF-8 LC_MESSAGES=ja_JP.UTF-8
[7] LC_PAPER=ja_JP.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=ja_JP.UTF-8 LC_IDENTIFICATION=C

attached base packages: [1] stats graphics grDevices utils datasets methods base

loaded via a namespace (and not attached): [1] compiler_4.1.2

packageVersion("sleuth") [1] ‘0.30.0’

packageVersion("tximport") [1] ‘1.22.0’

また、kallistoで作成された*.h5ファイルのサイズは次の通りです。 SRR1550989 65.7M SRR1551005 70.2M SRR1551011 64.7M SRR1551050 70.4M SRR1551057 68.3M SRR1551071 66.8M SRR1551091 57.5M

condaの仮想環境について、実行してみたところ、下記のようなエラーが出てしまいました。

library(sleuth) エラー: package or namespace load failed for ‘sleuth’: オブジェクト ‘h5write.default’ は 'namespace:rhdf5' によってエクスポートされていません

お忙しいところ、何度も申し訳ありませんが、アドバイスいただければ幸いです。

[RolyPolyCoily]

@nakane-scc ありがとうございます。 仮想環境 ngs を起動した状態（コマンドラインの先頭に (ngs) と明記された状態）で、conda install -c bioconda r-sleuth -y を再実行し、エラー等が出ていないかご確認いただけますでしょうか。 また念の為、 conda install -c bioconda bioconductor-rhdf5 -y も実行していただけますでしょうか。 あわせて which R と実行したときの結果も教えていただきたいです。

上記で特段エラーなどがなければあらためてRを起動させ、library(sleuth) をお試しください。 ちなみに tximport は不要でしたね…すみません。

[nakane-scc]

教えていただいたコマンドを実行すると、以下のようになりました。

(ngs) $ conda install -c bioconda r-sleuth -y Collecting package metadata (current_repodata.json): done Solving environment: done

All requested packages already installed.

(ngs) $ conda install -c bioconda bioconductor-rhdf5 -y Collecting package metadata (current_repodata.json): done Solving environment: done

All requested packages already installed.

(ngs) $ which R ~/miniconda3/envs/ngs/bin/R

特段エラーが出なかったので、Rを起動させ、library(sleuth)を試したところ、先程と同じエラーメッセージが出て、パッケージを読み込むことができませんでした。

[nakane-scc]

FedoraでSleuthを正常に動作させるのは難しいかもしれないと思い、Docker Image(quay.io/biocontainers/r-sleuth:0.30.0--r41hdfd78af_5)を使ってみることにしました。すると、警告メッセージは出たものの、 so <- sleuth_prep(s2c, extra_bootstrap_summary=T, target_mapping=t2g, aggregation_column='ens_gene', gene_mode=T) は一応、終了いたしました。

その後、

so <- sleuth_fit(so, ~condition, 'full') LRT <- sleuth_fit(so, ~1, 'reduced')

の両方で、

fitting measurement error models shrinkage estimation 1 NA values were found during variance shrinkage estimation due to mean observation values outside of the range used for the LOESS fit. The LOESS fit will be repeated using exact computation of the fitted surface to extrapolate the missing values. These are the target ids with NA values: ENSG00000156508 computing variance of betas

というメッセージが出ました。とりあえず、ここを無視して更に先に進め、

p <- plot_bootstrap(LRT, "ENST00000503567.5", units="est_counts", color_by="condition")

では、

Error in get_bootstrap_summary(obj, target_id, units) : couldn't find target_id 'ENST00000503567.5' In addition: Warning message: In check_quant_mode(obj, units) : your sleuth object is in gene mode, but you selected 'est_counts'. Selecting 'scaled_reads_per_base'...

というエラーが出たので、Target IDの確認できた「ENSG00000059804」を使って上記コマンドを実行すると、 Warning message: In check_quant_mode(obj, units) : your sleuth object is in gene mode, but you selected 'est_counts'. Selecting 'scaled_reads_per_base'... という警告文はありましたが、処理はできたようです。

正常に動作した場合の結果がわからないので、私の行った一連の処理で問題ないのか、警告メッセージは無視せず対処すべきもので、出力された結果が間違っているのか判断できません。いかがでしょうか。

[RolyPolyCoily]

@nakane-scc お返事が遅くなり申し訳ございません。 dockerで最後まで動いたとのことでよかったです。Linuxでのトラブルはこちらで再現できないためご不便おかけしました。 警告は無視して構わないようです。ただしNAとして除外された転写産物があったことは記録したほうがいいかも知れません。今回は１つだけですが、大量に落ちている場合や、解析で注目したい遺伝子だった場合などは検討が必要かと思います。 エラーは、sleuth_prep() にて gene_mode=Tとしているので gene-levelのデータを扱っていますが、plot_bootstrap() では ENSG（gene）ではなくENST（transcript）を指定していたためエラーが出たようです。