I am trying to assemble a genome (90x coverage, mean read quality: 30). The system on which I am working has 128 memory threads and 188 GB RAM. I am trying to assemble using the following command: hifiasm -o assembly -t 16 -f0 _______.sra.fastq, but it is getting killed. Even the log file is empty. I even tried to decrease the number of threads and even used -f0, but nothing is working out. How do I figure this out?
•
written
by
pranavdatar01 •
2 answers
The data that you are using (SRA accession SRX12366993) consists of PacBio continuous long reads generated on a Sequel I system. These reads are noisy and have a lower base accuracy compared to PacBio HiFi reads. The hifiasm assembler is designed specifically for PacBio HiFi reads, which explains why your command fails with no overlaps detected or the process getting killed due to excessive memory usage during error-prone overlap detection.
The original assembly in the associated publication (DOI: 10.3389/fevo.2022.850817) used both PacBio continuous long reads and 10x Genomics linked reads. For the PacBio data, they employed the Canu assembler, which is suitable for noisy long reads. They first converted PacBio BAM files to FASTA format using DEXTRACTOR, filtering out sequences shorter than 1,000 base pairs and those with quality below 0.8.
Since you wish to assemble using only the PacBio reads, I recommend using Canu with default parameters, as done in the publication. Your system has sufficient resources (188 GB RAM and 128 threads), but start with fewer threads to monitor memory usage. Assuming you have the reads in FASTQ format, first convert them to FASTA if needed, then run:
canu -p assembly -d canu_assembly genomeSize=1.5g -pacbio-raw your_reads.fastq
Replace your_reads.fastq with your actual file name. The -pacbio-raw flag specifies noisy PacBio reads. Monitor the process with tools like top or htop to check for memory issues. If memory exceeds limits, reduce the number of threads via the maxThreads parameter or subsample further.
If your data is still in BAM format from SRA, use samtools fastq to extract it first.
Kevin
The Data and the assembler: The SRA accession SRX12366993 is from a PacBio Sequel I system. This platform produces CLR (Continuous Long Reads), which have an error rate of ~10–15%. Hifasm is mathematically designed for HiFi reads (errors <1%). When it encounters noisy CLR data, it fails to find overlaps, bloats the memory usage trying to find them, and gets killed by the system (memory failure)
Match the Tool to the Platform: You should use eg pysradb (or your SRA discovery script) to check the model and platform before starting. If Platform == "HiFi": Use hifiasm or LJA. if Platform == "CLR" OR "ONT": Use Flye or Canu. These tools are "noise-tolerant" and will actually generate an assembly where hifiasm will simply crash.
Manage Coverage (The "Less is More" Rule) : High coverage is not always better for noisy reads. 90x coverage of 15% error reads creates a massive amount of computational "chatter." Subsample to ~40x–50x. This reduces the RAM requirement significantly and often results in a more contiguous assembly because the "noise" is reduced.
Log in to answer this question.
-
Canu Nanopore Genome assembly optimization
written by SomeOne •
-
Deduplication
written by tanmoymukherjee07 •
-
BLASTn produces empty output + output is 'Killed'
written by Shraddha •
-
Space Required to run 15 GB File in python
written by anasjamshed
-
AssemblyHaplotypesMultiThreadPlugin fails
written by kathrynm •
-
The memory is not enough when run cellranger to align the fq using the two parameters --localcores…
written by whaiyu06 •
-
memory allocation in metaSpades
written by Maorkn •
-
Platanus assembler : unable to start contig assembly, any problem with the command?
written by johnnytam100
-
MetaSpades can't access memory? Crashes when requesting 45 Gigs on a 128 Gig RHEL6 machine
written by jnowacki
-
GATK GenotypeGVCFs 8 weeks runtime
written by esticca •
How large do you expect the genome to be? For some human genomes that I tried it (i.e. genome size ~3Gbp), peak memory usage was ~160GB RAM (HiFi data, 50x coverage).
In my case, I ran
hifiasmlike this:I am expecting the genome to be around 1.5 Gb.
What kind of dataset is this? How many reads are there and what is the total data size?
As already noted by @Panos, the process is likely running out of memory (are you the only user on this system, if not other processes could reduce the available memory). Reducing the number of threads (you were alredy using a reasonable number to begin with) can only do so much.
These are PacBio HiFi reads (total reads 5,237,005), and the total data size is 255 GB.
I am the only user on the system.
Why don't you subsample your HiFi reads and try running hifiasm again? Maybe this will help if you have memory problems. However, I haven't done this myself recently, so I cannot really suggest a tool.
Also, it's weird that the log file is empty... Do you really get nothing in the log? In that case try I'd do some sanity checks. For example, download some public PacBio data from a relatively small genome (e.g. Drosophila) and run hifiasm with that data set...
I did subsample the data to 30x coverage. I did run hifiasm on it, the log file did show some contents, but hifiasm ended again, stating that there were no overlapping reads :(. I am currently trying on 60x and above subsamples now.
No overlapping reads? Have you checked for contamination?
Yes, I did check for adaptor contamination, and did not find any. I have kept all the reads with quality score of more than 30.
I don't mean adapter contamination. I mean contamination of your sample with DNA from other organisms. In my previous lab, for example, we were sequencing insect samples and some of them were contaminated by plant DNA; what the insects were eating. Some other times there's contamination at the sequencing facility by other samples that are processed along yours.
And if it's not a contamination, it could also be a bad DNA prep; relatively fragmented DNA, instead of the high molecular weight DNA that is typically needed for long read sequencing.
All valid points. This is a SRA dataset. If it is covered by a publication hopefully it is not a bad one. @Pranav if you can post the SRA accession # we can take a look at the metadata.
The SRA accession is SRX12366993
Are you trying to reproduce the assembly result?
If not, the complete assembly is already available at NCBI: https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_023646645.1/
Looks like there are two data types and you are using the smaller of the two datasets. https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=766478
Other data is 10x long range genome sequence (method now deprecated by 10x genomics), so that data may be harder to work with since it will require
longrangersoftware.I am trying to reproduce the assembly from the raw pacbio reads.
I want to reconstruct the assembly using the pacbio reads only.
Is there a paper associated with this data? If the original publication/submitters used both types of data to make their assembly then there must be a reason for a it. This does not appear to be PacBio HiFi data, which is superior quality and can lead to good assemblies on its own.
Were you able to make any progress using subsamples?
This is the paper: https://doi.org/10.3389/fevo.2022.850817 (Genome of Varanus salvator macromaculatus (Asian Water Monitor) Reveals Adaptations in the Blood Coagulation and Innate Immune System)
I tried running hifiasm on the subsamples. The job was getting killed for all the subsamples above 30X, and not on 30X. However, hifiasm was not able to find out any overlaps on the 30X subsample.
I then went on using the -f38 flag on the whole 91X sample, which essentially means that any k-mer appearing fewer than 38 times will be ignored. However, hifiasm was not able to find out any overlaps. I tried to decrease the -f below 38, but it was getting killed below -f38 :(
For assembly reproducibility, I would try to preserve not only the raw reads and final assembly, but also the exact assembler version, command-line parameters, input checksums, reference files if reference-assisted assembly is used, intermediate outputs, logs, memory/CPU information, and final assembly checksums.
In many cases, the issue is not only whether the assembly can be completed, but whether the exact run can be reviewed or reproduced later under the same conditions.