Comments (36)
I come back to report the assembly was quite succesful.
However, now I am trying to assemble reads from a different experiment (traditional RNAseq, single end reads, downloaded from SRA) and face a similar issue to here:
(RNA_Bloom) [lab_eh_2@NeuroServer denovo_transcriptome]$ rnabloom -fpr 0.01 -sef /home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz -ntcard -t 30 -outdir output
RNA-Bloom v1.4.3
args: [-fpr, 0.01, -sef, /home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz, -ntcard, -t, 30, -outdir, output]
name: rnabloom
outdir: output
WARNING: Output directory does not exist!
Created output directory at `output`
K-mer counting with ntCard...
Running command: `ntcard -t 30 -k 25 -c 65535 -p output/rnabloom @output/rnabloom.ntcard.readslist.txt`...
Parsing histogram file `output/rnabloom_k25.hist`...
Unique k-mers (k=25): 262.333.277
Unique k-mers (k=25,c>1): 129.744.533
K-mer counting completed in 31.585s
Bloom filters Memory (GB)
====================================
de Bruijn graph: 0.57971644
k-mer counting: 2.2937248
paired k-mers (reads): 0.57971644
paired k-mers (frags): 0.57971644
screening: 0.57971644
====================================
Total: 4.612591
> Stage 1: Construct graph from reads (k=25)
Sampling read lengths (l>=25, n=1000) from each file...
min Q1 M Q3 max file
28 50 51 51 51 /home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz
Weighted [Q1,M,Q3]: [50, 51, 51]
Absolute [min,max]: [28, 51]
Paired k-mers distance: 15
Max. tip length: 26
Parsing `/home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz`...
Parsed 25.694.006 sequences in 2m 14s
DBG Bloom filter FPR: 0.998 %
Counting Bloom filter FPR: 1.01 %
Reads paired k-mers Bloom filter FPR: 0.335 %
> Stage 1 completed in 2m 23s
> Skipping Stage 2 for single-end reads.
> Stage 3: Assemble transcripts for "rnabloom"
ERROR: null
java.lang.NullPointerException
at rnabloom.RNABloom$SingleEndReadsIterator.<init>(RNABloom.java:4679)
at rnabloom.RNABloom.assembleSingleEndReads(RNABloom.java:4805)
at rnabloom.RNABloom.assembleTranscriptsSE(RNABloom.java:5379)
at rnabloom.RNABloom.main(RNABloom.java:7104)
How do I specify the reads when I only have 1 file?
from rna-bloom.
Don´t worry, fixed it!
Thanks!
from rna-bloom.
Finally finished succesfully, thanks a lot for your help @kmnip !
from rna-bloom.
As mentioned previously, you need to trim the UMIs and barcodes entirely from the reads with umitools extract
. After that, only right (_2.fastq
) reads contain the sequences that can be assembled.
You have several options here. You can align the right reads (with bwa/minimap2) to either your CDS, reference transcripts, or the reference genome. For example, the gene for "Bgal" should be GLB1 (chr3:32,996,617-33,097,146 in GRCh38).
After alignment, you can extract the aligned reads with samtools fastq
.
samtools fastq -F 0x4 alignment.bam > extracted_reads.fastq
You can assemble the extracted reads with RNA-Bloom:
java -jar RNA-Bloom.jar -sef extracted_reads.fastq -ntcard -t THREADS -outdir OUTDIR
from rna-bloom.
A couple of doubts:
-Read2 should not contain any UMI or Barcode right? I mean, my average insert size is 400 bp, read2 is 88 bases long...unlikely any of them included an insert small enough to get to the other side ,right?
-If I align my reads to the CDS (thats the only thing I have by now) I will keep only those reads that map to the CDS and will discard reads mapping outside...precisely what I am looking for: those reads that map to the 3'UTR....or am I missing something?
from rna-bloom.
Based on the sample data from 10x, I don't think the right read (98-bp) contains any UMI/barcodes at all. But if you believe that yours do, please do remove them before alignment/assembly.
The insert sizes probably peak at 400 bp, but they should vary a bit.
Yes, you are correct about the CDS missing the UTR. That's why I suggested reference transcript and genome as alternatives.
from rna-bloom.
The problem is that being a transgene (a sequence inserted artificially in the mouse genome in order to be expressed by select cells so we can identify them), the sequence is not present in any reference....I only know the CDS. I need a way to select the reads that map to the CDS and use those to "build" a transcript by adding more reads (specially in the direction 3' of the CDS)...thats why I initially thought Kollector could do it.
from rna-bloom.
You can just assemble the trimmed reads and then extract the transcripts that align to your CDS?
from rna-bloom.
Mmm I guess so. You mean:
-Discard read1
-Trim any UMI, cell barcode, etc from read2
-Perform de novo transcriptome assembly with ALL read2 using RNA-bloom
-BLAST the assembled transcriptome against my CDS
-Pick candidate transcripts and identify the 3'UTR from those long enough
?
from rna-bloom.
Yes, that seems like your only option because your left read most likely contain only UMIs + cell barcodes.
Kollector would not work in this case because its recruitment step expects paired end reads where one paired read would anchor on a given sequence.
from rna-bloom.
Alright.
Hope RNA-bloom will not crash!
Then my next question is: I have 2 samples (=libraries, from different animals) should I supply RNA-bloom with the 2 read2 files, or better do it one by one?
from rna-bloom.
Let me know if you run into any issues.
If they are different animals, then you should assemble them separately.
from rna-bloom.
Thanks a lot, Ka!
Will keep you posted
from rna-bloom.
Unfortunately, I ran into something...
(RNA_Bloom) [lab_eh_2@NeuroServer RNA_Bloom]$ rnabloom -fpr 0.01 -ser /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_unpaired_norRNA_2.fq.gz /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_val_norRNA_2.fq.gz -ntcard -t 30 -outdir output
RNA-Bloom v1.4.3
args: [-fpr, 0.01, -ser, /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_unpaired_norRNA_2.fq.gz, /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_val_norRNA_2.fq.gz, -ntcard, -t, 30, -outdir, output]
name: rnabloom
outdir: output
K-mer counting with ntCard...
Running command: `ntcard -t 30 -k 25 -c 65535 -p output/rnabloom @output/rnabloom.ntcard.readslist.txt`...
Parsing histogram file `output/rnabloom_k25.hist`...
Unique k-mers (k=25): 1.990.801.340
Unique k-mers (k=25,c>1): 802.505.180
K-mer counting completed in 13m 57s
Bloom filters Memory (GB)
====================================
de Bruijn graph: 4.399367
k-mer counting: 14.187311
paired k-mers (reads): 4.399367
paired k-mers (frags): 4.399367
screening: 4.399367
====================================
Total: 31.784777
> Stage 1: Construct graph from reads (k=25)
Sampling read lengths (l>=25, n=1000) from each file...
min Q1 M Q3 max file
61 87 88 88 88 /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_unpaired_norRNA_2.fq.gz
78 87 88 88 88 /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_val_norRNA_2.fq.gz
Weighted [Q1,M,Q3]: [87, 88, 88]
Absolute [min,max]: [61, 88]
Paired k-mers distance: 52
Max. tip length: 63
Parsing `/home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_unpaired_norRNA_2.fq.gz`...
Parsed 100.549 sequences in 1.42s
Parsing `/home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_val_norRNA_2.fq.gz`...
Parsed 345.166.566 sequences in 1h 6m 3s
Parsed 345.267.115 sequences from 2 files.
DBG Bloom filter FPR: 0.999 %
Counting Bloom filter FPR: 1.01 %
Reads paired k-mers Bloom filter FPR: 0.25 %
> Stage 1 completed in 1h 6m 50s
> Skipping Stage 2 for single-end reads.
> Stage 3: Assemble transcripts for "rnabloom"
ERROR: null
java.lang.NullPointerException
at rnabloom.RNABloom$SingleEndReadsIterator.<init>(RNABloom.java:4670)
at rnabloom.RNABloom.assembleSingleEndReads(RNABloom.java:4805)
at rnabloom.RNABloom.assembleTranscriptsSE(RNABloom.java:5379)
at rnabloom.RNABloom.main(RNABloom.java:7104)
(RNA_Bloom) [lab_eh_2@NeuroServer RNA_Bloom]$
´´´
Any ideas whats happening here?
from rna-bloom.
It looks like you have found a bug.
You can bypass it by using both -ser
and -sef
, i.e.
rnabloom -fpr 0.01 \
-ser /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_unpaired_norRNA_2.fq.gz \
-sef /home/LAB/lab_eh_2/AER_Data/snRNAseq/C1/trimmed/Clean_rRNA/A_S1_L001_R2_001_val_norRNA_2.fq.gz \
-ntcard -t 30 -outdir output
Since your reads are not strand-specific, it doesn't matter whether you specify your reads with -ser
or -sef
.
I will fix this in the next release. Sorry about that!
from rna-bloom.
Alright, good to know it was not something worse. I was already touching the JAVA options.
So you are saying the problem was I specified 2 input read files with only one argument, right? It should have been -ser X -sef Y, not -ser X Y. But then that means the bug is that you cannot, so far, provide more than 1 read file per "modality" -left, -right, -sef, -ser, correct?
Also, I was wondering if I could use the -mergepool option here if I only have single end reads: how should I specify the files from different samples in the readslist.txt
file? Given I need to put them in the 4th column, do I just leave the columns 2 and 3 empty?
Sample1 sef_file ser_file
Sample2 sef_file ser_file
from rna-bloom.
No, the bug was that -ser
and -sef
had to be used together.
In in version 1.4.3, the -pool
option doesn't support single-end reads (yet):
https://github.com/bcgsc/RNA-Bloom/releases/tag/v1.4.3
support for single-end reads in -pool will be implemented in a future release
It will be available in the next release though.
from rna-bloom.
Ok, thanks for the info!
Is it possible to use the -mergepool parameter if I provide several manually assembled transcriptomes?
I am running it now with the first sample, fingers crossed!
from rna-bloom.
Just wondering, is it normal that it will take more than 9 hours to assemble a transcriptome from 345.267.115 sequences using 30 threads? (Currently at "Reducing redundancy in assembled transcripts...").
from rna-bloom.
At this stage, it is reducing the redundancy in the assembled sequences to generate *.transcripts.nr.fa
, but you can already use the sequences in *.transcripts.fa
.
345 million reads is not a small amount data and assembly of single-end reads is less efficient than paired-end assembly. So, that is kind of expected.
from rna-bloom.
Alright, thanks Ka!
Regarding output till now, I guess the transcripts.short.fa
file contains the assembled transcripts smaller than certain threshold, right? Is it taking those into account too for the redundancy reduction?
from rna-bloom.
Yes, the default length threshold is 200.
Those short sequences are not included in the redundancy reduction process.
from rna-bloom.
Yes, that was the same bug as before. Sorry for the inconvenience.
You need to split your FASTQ file in half and use -ser
and -sef
options together.
For example, you can split the file like so:
zcat SRR7039666_3_trimmed.fq.gz | split -l 51400000 - partition_
mv partition_aa SRR7039666_3_trimmed.01.fq
mv partition_ab SRR7039666_3_trimmed.02.fq
gzip SRR7039666_3_trimmed.01.fq SRR7039666_3_trimmed.02.fq
rnabloom -fpr 0.01 -sef SRR7039666_3_trimmed.01.fq.gz -ser SRR7039666_3_trimmed.02.fq.gz -ntcard -t 30 -outdir output
from rna-bloom.
Ok, that was the smart move.
I just ran rnabloom with the same reads file twice like:
rnabloom -fpr 0.01 -sef /home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz -ser /home/LAB/lab_eh_2/AER_Data/Chamessian/3/trimmed/SRR7039666_3_trimmed.fq.gz -ntcard -t 30 -outdir output
Is that crazy? The file is small so the assembly already finished, can I use those results or better not?
from rna-bloom.
It is okay if noise is not an issue.
Since the read file is relatively small, I suggest you re-run the assembly with the split files.
from rna-bloom.
Alright!
Thanks!
from rna-bloom.
I just realized that I messed up the number for the example split
command previously.
25,694,006 reads * 4 lines/read = 102,776,024 lines
Split that in half would be 51,388,012 lines per partition.
So, something like -l 51400000
would yield exactly two files.
from rna-bloom.
Hi Ka,
Not sure if I should open a new issue or keep adding here, because this is a new error not discussed previously. If you prefer, I will close this and post as new issue. In the meantime:
I have been de novo assembling several public datasets. In this case from a RNAseq experiment (not single cell). The only difference I can see is that the size of the reads files is way bigger than previous (I have merged together 33 SRR single end files and then split in 2 as recommended above: each reads file is about 39GB). And I get this:
(RNA_Bloom) [lab_eh_2@NeuroServer trimmed]$ rnabloom -fpr 0.01 -sef *_00_trimmed.fq.gz -ser *_01_trimmed.fq.gz -ntcard -t 30 -outdir output
RNA-Bloom v1.4.3
args: [-fpr, 0.01, -sef, Chongtam_00_trimmed.fq.gz, -ser, Chongtam_01_trimmed.fq.gz, -ntcard, -t, 30, -outdir, output]
name: rnabloom
outdir: output
K-mer counting with ntCard...
Running command: `ntcard -t 30 -k 25 -c 65535 -p output/rnabloom @output/rnabloom.ntcard.readslist.txt`...
ERROR: Error running ntCard!
Not a very informative error...any ideas??
from rna-bloom.
Hi,
Can you please run this command and see whether you got an error from ntCard?
ntcard -t 30 -k 25 -c 65535 -p output/rnabloom @output/rnabloom.ntcard.readslist.txt
Thanks!
from rna-bloom.
Mmmm this is what I get...
(RNA_Bloom) [lab_eh_2@NeuroServer trimmed]$ ntcard -t 30 -k 25 -c 65535 -p output/rnabloom@output/rnabloom.ntcard.readslist.txt
gzip: Chongtam_00_trimmed.fq.gz: unexpected end of file
PID 3694 exited with status 1
from rna-bloom.
That means your fastq file Chongtam_00_trimmed.fq.gz
is truncated.
from rna-bloom.
Ok.
I made it by concatenating several single end fastqs, counting lines and dividing in half as previously advised. How could the file be truncated?
from rna-bloom.
That error actually came from gzip
. ntCard was simply using it to decompress your FASTQ file. So, your gzip file is corrupted according to gzip
.
Try this and see if it is really truncated:
zcat Chongtam_00_trimmed.fq.gz | tail -n 4
You can try decompressing and then compressing the file again.
from rna-bloom.
Result:
(RNA_Bloom) [lab_eh_2@NeuroServer trimmed]$ zcat Chongtam_00_trimmed.fq.gz | tail -n 4
gzip: Chongtam_00_trimmed.fq.gz: unexpected end of file
@SRR13563559.53035757 NB501946:387:HNHTGBGXC:1:23107:17045:17168 length=70 CACATATTCATAATCAGGGTCTTGCAGCAATTGAGAAGTGTGGAGAAGCTGCTCTAAATGACCCCA
+
AA
(RNA_Bloom) [lab_eh_2@NeuroServer trimmed]$
How did this happen??
from rna-bloom.
That definitely looks truncated.
Check your original FASTQ file that contains the read SRR13563559.53035757
.
Otherwise, you may have had multiple commands running at the same time and they write to the same file.
from rna-bloom.
You were right, the original reads file:
(in202109) [lab_eh_2@NeuroServer Chongtam]$ zcat SRR13563559.fastq.gz | tail -n 4
@SRR13563559.68742091 NB501946:387:HNHTGBGXC:1:11311:18817:13247 length=70
TAGGATTGTTTGTACGGAATAGAGAAAAGTGTGGCAGAAAAGCGGAACGCGCAAAGAGAACGTATATGGT
+SRR13563559.68742091 NB501946:387:HNHTGBGXC:1:11311:18817:13247 length=70
AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
So I am reconstructing again, will report back!
Thanks!
from rna-bloom.
Related Issues (20)
- Transcriptome assembly from ONT direct RNA sequencing seems to be incomplete HOT 9
- conda install url not working HOT 1
- Inputting multiple long-read files at once HOT 3
- Java out of memory exception HOT 3
- Exception in thread "main" java.lang.OutOfMemoryError: Java heap space HOT 4
- RNABloom not creating transcriptome files HOT 3
- Errors (reading file and I/O) for one of 9 pooled assemblies HOT 11
- can rnabloom run from stage 3? HOT 2
- question about the output HOT 2
- Can output be used for DE analysis? HOT 2
- Question about assemble long-read sequencing data with short-read polishing HOT 3
- RAM requirements. What is the minimum? HOT 2
- Consensus transcriptome de novo with multiple samples HOT 1
- reduce redundancy in direct RNA long-read only assembly? HOT 1
- error in assembling ONT DRS long reads HOT 1
- assembling nanopore sequenced bacterial transcriptomes HOT 1
- Improving Sensitivity for Similar Isoform Assembly and Handling PacBio HiFi RNA Reads Without Polishing HOT 3
- is the assembly deterministic? HOT 3
- Long reads genome guided HOT 1
- Best practices for pre-processing prior to assembly HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from rna-bloom.