Closed Ci-TJ closed 3 years ago
Hi Qin,
Thanks for using Terminus. It seems Terminus panicked while parsing the salmon_infRep20_out/quant_out/SRR6868519/aux_info/eq_classes.txt.gz
file. My guess is that Terminus is not getting the format it expects from the equivalence class file. If this is made from the experiment SRR6868519
, I can try to reproduce the error when I get a chance. Before that can you please tell me what version of Salmon are you using and which index the salmon is run with.
@hiraksarkar Thank you very much for your timely response. I use mouse Ensemble release 99 cdna(ftp://ftp.ensembl.org/pub/release-99/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz) and ncrna(ftp://ftp.ensembl.org/pub/release-99/fasta/mus_musculus/ncrna/Mus_musculus.GRCm38.ncrna.fa.gz) to build index by salmon(v 1.2.1), which built index with k=31 and other parameters by default. I will send my index to your email. By the way, I have update my version of salmon to v 1.3.0 in my current conda env. Does terminus call salmon during runtime?
Hi @Ci-TJ ,
It does not call Salmon internally and salmon(v 1.2.1) should ideally work. I would definitely check SRR6868519
when I get a chance and get back to you if I can reproduce the result.
Hi @Ci-TJ ,
Can you please re-run the salmon with the following command in your script
salmon quant -i /home2/ymwang/linqin/RNA-seq/fasta/release99/mouse_rna_salmon_k31_index -l A -1 ../download_GSE112055/$line".sra_1.fastq" -2 ../download_GSE112055/$line".sra_2.fastq" -p 32 --validateMappings -o quant_out/$line --numGibbsSamples 20 --seqBias --gcBias -d
Notice the use of -d
instead --dumpEq
, which does not provide rich equivalence classes (equivalence classes with weights), a requirement for terminus
.
Hi @hiraksarkar It's no problem in group, but I got trouble in collapse.
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ cat salmon_quant_terminus.sh
### quant for terminus
pwd
salmon -v
terminus --version
date
for line in $(head -n 4 ../Run_GSE112055)
do
salmon quant -i /home2/ymwang/linqin/RNA-seq/fasta/release99/mouse_rna_salmon_k31_index -l A \
-1 ../download_GSE112055/$line".sra_1.fastq" -2 ../download_GSE112055/$line".sra_2.fastq" -p 32 --validateMappings -o quant_out/$line --numGibbsSamples 20 --seqBias --gcBias -d && echo $line Done!
terminus group -d quant_out/$line -o group_out && echo group $line OK!
done
date(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ ll group_out/SRR6868519/
total 159240
drwxrwxr-x 2 ymwang ymwang 118 Nov 29 00:13 ./
drwxrwxr-x 6 ymwang ymwang 94 Nov 29 00:38 ../
-rw-rw-r-- 1 ymwang ymwang 438614 Nov 29 00:13 collapsed.log
-rw-rw-r-- 1 ymwang ymwang 97092240 Nov 29 00:13 delta.log
-rw-rw-r-- 1 ymwang ymwang 62018524 Nov 29 00:11 die_roll.log
-rw-rw-r-- 1 ymwang ymwang 72122 Nov 29 00:13 groups.txt
-rw-rw-r-- 1 ymwang ymwang 3427455 Nov 29 00:11 infrv.log
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ terminus collapse -h
terminus-collapse
analyze a collection of per-sample groups, and produce a consensus grouping.
USAGE:
terminus collapse [OPTIONS] --dirs <dirs>... --out <out>
FLAGS:
-h, --help Prints help information
-V, --version Prints version information
OPTIONS:
-c, --consensus-thresh <consensus-thresh> threshold for edge consensus [default: 0.5]
-d, --dirs <dirs>... direcotories to read the group files from
-o, --out <out> prefix where output would be written
-t, --threads <threads>
number of threads to use when writing down the collapsed quantifications [default: 8]
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ terminus collapse -d group_out/SRR6868519/ group_out/SRR6868520/ group_out/SRR6868521/ group_out/SRR6868522/ -o collapse_out
experiment name
thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: Os { code: 2, kind: NotFound, message: "No such file or directory" }', src/util.rs:267:38
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ terminus collapse -d $PWD/quant_out/SRR6868519/ $PWD/quant_out/SRR6868520/ $PWD/quant_out/SRR6868521/ $PWD/quant_out/SRR6868522/ -o collapse_out
experiment name
# targets : 137934
did serialize eq classes : true
# boot : 20
thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: Os { code: 2, kind: NotFound, message: "No such file or directory" }', src/main.rs:224:63
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$ RUST_BACKTRACE=1 terminus collapse -d $PWD/quant_out/SRR6868519/ $PWD/quant_out/SRR6868520/ $PWD/quant_out/SRR6868521/ $PWD/quant_out/SRR6868522/ -o collapse_out
experiment name
# targets : 137934
did serialize eq classes : true
# boot : 20
thread 'main' panicked at 'called `Result::unwrap()` on an `Err` value: Os { code: 2, kind: NotFound, message: "No such file or directory" }', src/main.rs:224:63
stack backtrace:
0: backtrace::backtrace::libunwind::trace
at /cargo/registry/src/github.com-1ecc6299db9ec823/backtrace-0.3.46/src/backtrace/libunwind.rs:86
1: backtrace::backtrace::trace_unsynchronized
at /cargo/registry/src/github.com-1ecc6299db9ec823/backtrace-0.3.46/src/backtrace/mod.rs:66
2: std::sys_common::backtrace::_print_fmt
at src/libstd/sys_common/backtrace.rs:78
3: <std::sys_common::backtrace::_print::DisplayBacktrace as core::fmt::Display>::fmt
at src/libstd/sys_common/backtrace.rs:59
4: core::fmt::write
at src/libcore/fmt/mod.rs:1076
5: std::io::Write::write_fmt
at src/libstd/io/mod.rs:1537
6: std::sys_common::backtrace::_print
at src/libstd/sys_common/backtrace.rs:62
7: std::sys_common::backtrace::print
at src/libstd/sys_common/backtrace.rs:49
8: std::panicking::default_hook::{{closure}}
at src/libstd/panicking.rs:198
9: std::panicking::default_hook
at src/libstd/panicking.rs:217
10: std::panicking::rust_panic_with_hook
at src/libstd/panicking.rs:526
11: rust_begin_unwind
at src/libstd/panicking.rs:437
12: core::panicking::panic_fmt
at src/libcore/panicking.rs:85
13: core::option::expect_none_failed
at src/libcore/option.rs:1269
14: terminus::main
15: std::rt::lang_start::{{closure}}
16: std::rt::lang_start_internal::{{closure}}
at src/libstd/rt.rs:52
17: std::panicking::try::do_call
at src/libstd/panicking.rs:348
18: std::panicking::try
at src/libstd/panicking.rs:325
19: std::panic::catch_unwind
at src/libstd/panic.rs:394
20: std::rt::lang_start_internal
at src/libstd/rt.rs:51
21: main
22: __libc_start_main
23: <unknown>
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus]$
Do I misunderstand with input of collapse?
Hi @Ci-TJ ,
There is a directory dependency between group
and collapse
, So your group
command should write the output in the same location where you want to write your collapse output in,
For example if you have three experiments, SRR1852860-81
terminus group -m 0.05 -d quant/SRR1852860/ -o term_out/SRR1852860
terminus group -m 0.05 -d quant/SRR1852861/ -o term_out/SRR1852861
terminus group -m 0.05 -d quant/SRR1852881/ -o term_out/SRR1852881
Your collapse command will be as follows,
terminus collapse -d quant/SRR1852860 quant/SRR1852861 quant/SRR1852881 -c 1.0 -o term_out
In your case, you can do the following,
terminus collapse -d group_out/SRR6868519/ group_out/SRR6868520/ group_out/SRR6868521/ group_out/SRR6868522/ -o group_out
@hiraksarkar Thank you for you help during this time! I have reproduced your codes, and I have a question if I can feed directly the outpout of collapse to fishpond for differential expression analysis. In any case, I will try later.
Best
Hi @Ci-TJ
Glad to hear. You can use the code directly in DE-seq2 by following the tutorial, https://combine-lab.github.io/terminus-tutorial/2020/running-terminus/
Yeah, I know how to use the DESeq2 for DE analysis, but fishpond is special and easy-to-use for salmon; If outputs of terminus can be fed directly to fishpond, it's very wonderful!
Hi @Ci-TJ
I just had a brief conversation with Mike Love (who wrote fishpond), Can you please put an issue there with how you wanna try fishpond. He might add some added code to make that work. Feel free to tag me too.
Thanks
Hi @hiraksarkar @mikelove I have tried to import the outputs of teminus directly by tximeta, but failed.
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus/fishpond]$ cat ../salmon-terminus.sh
### quant and terminus
echo ==================================================================================
echo ==================================================================================
pwd
salmon -v
terminus --version
date
for line in $(cat ../Run_GSE112055)
do
salmon quant -i /home2/ymwang/linqin/RNA-seq/fasta/release99/mouse_rna_salmon_k31_index -l A \
-1 ../download_GSE112055/$line".sra_1.fastq" -2 ../download_GSE112055/$line".sra_2.fastq" -p 32 --validateMappings -o quant/$line --numGibbsSamples 20 --seqBias --gcBias -d && echo $line Done!
terminus group -m 0.5 -d quant/$line -o term_out && echo group $line OK!
done
# for collapse
ls quant
terminus collapse -d quant/SRR6868519 quant/SRR6868520 quant/SRR6868521 quant/SRR6868522 quant/SRR6868523 quant/SRR6868524 quant/SRR6868525 quant/SRR6868526 quant/SRR6868527 quant/SRR6868528 -c 1.0 -o term_out
date
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus/fishpond]$ ll ../term_out/SRR6868519/
total 369552
drwxrwxr-x 3 ymwang ymwang 4096 Nov 30 10:28 ./
drwxrwxr-x 12 ymwang ymwang 4096 Dec 2 19:34 ../
drwxrwxr-x 3 ymwang ymwang 55 Nov 30 10:28 aux_info/
-rw-rw-r-- 1 ymwang ymwang 50441 Dec 2 19:35 clusters.txt
-rw-rw-r-- 1 ymwang ymwang 465 Dec 2 19:35 cmd_info.json
-rw-rw-r-- 1 ymwang ymwang 423733 Dec 2 18:17 collapsed.log
-rw-rw-r-- 1 ymwang ymwang 94214845 Dec 2 18:17 delta.log
-rw-rw-r-- 1 ymwang ymwang 274633069 Dec 2 18:15 die_roll.log
-rw-rw-r-- 1 ymwang ymwang 69588 Dec 2 18:17 groups.txt
-rw-rw-r-- 1 ymwang ymwang 3424114 Dec 2 18:15 infrv.log
-rw-rw-r-- 1 ymwang ymwang 5578342 Dec 2 19:35 quant.sf
-rw-rw-r-- 1 ymwang ymwang 68 Dec 2 19:35 terminus_info.json
(R40) [ymwang @ ~/linqin/RNA-seq/Analysis/GSE112055/Terminus/fishpond]$
> library(tximeta)
> library(fishpond)
> suppressPackageStartupMessages(library(SummarizedExperiment))
> dir<- "/home2/ymwang/linqin/RNA-seq/Analysis/GSE112055"
> coldata <- read.table("coldata.txt", sep="\t", header=T)
> head(coldata)
names BioSample drug_treatment Experiment GEO_Accession Sample.Name
1 SRR6868519 SAMN08741761 None SRX3822981 GSM3056373 GSM3056373
2 SRR6868520 SAMN08741758 None SRX3822982 GSM3056374 GSM3056374
3 SRR6868521 SAMN08741757 None SRX3822983 GSM3056375 GSM3056375
4 SRR6868522 SAMN08741754 None SRX3822984 GSM3056376 GSM3056376
5 SRR6868523 SAMN08741749 None SRX3822985 GSM3056377 GSM3056377
6 SRR6868524 SAMN08741742 None SRX3822986 GSM3056378 GSM3056378
surgery X
1 Sham NA
2 Sham NA
3 Sham NA
4 Sham NA
5 Sham NA
6 TAC NA
> coldata$files <- file.path(dir, "Terminus", "term_out",coldata$names, "quant.sf")
> file.exists(coldata$files)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> head(coldata)
names BioSample drug_treatment Experiment GEO_Accession Sample.Name
1 SRR6868519 SAMN08741761 None SRX3822981 GSM3056373 GSM3056373
2 SRR6868520 SAMN08741758 None SRX3822982 GSM3056374 GSM3056374
3 SRR6868521 SAMN08741757 None SRX3822983 GSM3056375 GSM3056375
4 SRR6868522 SAMN08741754 None SRX3822984 GSM3056376 GSM3056376
5 SRR6868523 SAMN08741749 None SRX3822985 GSM3056377 GSM3056377
6 SRR6868524 SAMN08741742 None SRX3822986 GSM3056378 GSM3056378
surgery X
1 Sham NA
2 Sham NA
3 Sham NA
4 Sham NA
5 Sham NA
6 TAC NA
files
1 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868519/quant.sf
2 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868520/quant.sf
3 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868521/quant.sf
4 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868522/quant.sf
5 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868523/quant.sf
6 /home2/ymwang/linqin/RNA-seq/Analysis/GSE112055/Terminus/term_out/SRR6868524/quant.sf
> se <- tximeta(coldata)
importing quantifications
reading in files with read_tsv
1 2 3 4 5 6 7 8 9 10
Error in if (!is.na(m)) { : argument is of length zero
By the way, I can import the outputs of salmon's quant to fishpond by tximeta. What's the solution?
Can you try with tximport instead? Or with skipMeta=TRUE
? The crux is that tximeta typically works by attaching the transcript ranges to the quantification data, but after Terminus, the rows are no longer transcripts but transcript groups. I'd like to develop a solution for Terminus output but this doesn't exist yet in tximeta.
@mikelove It's working when I import outputs of salmon by tximeta with skipMeta=TRUE. I didn't try tximport, because I thought it would be work too. For now, there are some problems for processing DE annalysis for gene-level, but it's doesn't matter. I will look forward to next update of tximeta. Thanks
Hi. I just installed the Terminus by conda. I didn't get any results after I run terminus with "group". Unfortunately, I know nothing about rust..
Best, Ci