Closed devonorourke closed 6 years ago
warrenlr
commented
6 years ago Building the Bloom filter is usually the fast part..
I suggest you use the writeBloom.pl utility in the tools folder. This should give you some information about the health of the Bloom::Filter module you compiled (and you can try running tests on your system using smaller genomes). Keep me posted. Rene
devonorourke
commented
6 years ago It looks like the program installed by our admin generates this error when I run the writeBloom.pl script:
-bash: ./writeBloom.pl: /home/martink/perl/5.10.0/bin/perl: bad interpreter: No such file or directoryPerhaps that's an indication that the compilation wasn't successful?
warrenlr
commented
6 years ago Above is caused by error while looking for interpreter that doesn't exist. This was in error, I had forgotten to replace the 1st line.
In commit ed8ad09 I replaced the shebang line by /usr/bin/env perl Hopefully it points to the same version as the one used to build the Bloom::Filter
devonorourke
commented
6 years ago Hi Rene, While I'm waiting on my system administrator to edit the shebang I ran a local install of v1.8.5 and it completed without issue. I went from 11,654 scaffolds to 11,638; I'm surprised that was all that was tightened up but perhaps I can now look back and determine what settings are best to modify. Thanks
warrenlr
commented
6 years ago the power of links is in scaffolding using various distance constraints, iteratively.
1) Running links multiple times (when working with large long read dataset / genomes are very big)
eg.
./LINKS -f $2 -s ont.fof -b links1 -d 1000 -t 10 -k $1 -l $3 -a $4 ./LINKS -f links1.scaffolds.fa -s ont.fof -b links2 -d 2500 -t 5 -k $1 -l $3 -a $4 -o 1 -r links1.bloom ./LINKS -f links2.scaffolds.fa -s ont.fof -b links3 -d 5000 -t 5 -k $1 -l $3 -a $4 -o 2 -r links1.bloom ./LINKS -f links3.scaffolds.fa -s ont.fof -b links4 -d 7500 -t 4 -k $1 -l $3 -a $4 -o 3 -r links1.bloom ./LINKS -f links4.scaffolds.fa -s ont.fof -b links5 -d 10000 -t 4 -k $1 -l $3 -a $4 -o 4 -r links1.bloom ./LINKS -f links5.scaffolds.fa -s ont.fof -b links6 -d 12500 -t 3 -k $1 -l $3 -a $4 -o 5 -r links1.bloom ./LINKS -f links6.scaffolds.fa -s ont.fof -b links7 -d 15000 -t 3 -k $1 -l $3 -a $4 -o 6 -r links1.bloom ./LINKS -f links7.scaffolds.fa -s ont.fof -b links8 -d 30000 -t 2 -k $1 -l $3 -a $4 -o 7 -r links1.bloom
2) Running links iteratively in a single command (works best when genomes are small/RAM not limiting)
./LINKS -f $2 -s ont.fof -b links1 -d 1000,5000,7500,10000,12500,15000,30000 -t 10,5,5,4,4,3,3,2 -k $1 -l $3 -a $4
*the value of -t will determine the #kmer pairs extracted and incidentally the RAM used. this will vary from one dataset to the next. I also recommend correcting ONT reads if you can, it will allow you to choose higher k values and increase the specificity. If not, no big deal (but are stuck choosing k15-17).
devonorourke
commented
6 years ago Hi again René,
I've gained access to a computer that has lots of memory and recently ran a LINKS job following what I think emulates the parameters you've suggested above with the beluga.fa script. Perhaps I'm mistaking some critical parameters, but I believe I've set things in the job as you've described...
The job failed despite allocating 1.2 Tb of memory, which makes me think that I'm doing something wrong right at the outset - do most LINKS jobs require this much memory? The job failed when scaffolding, after building the initial bloom filter. The command was:
Running: /pylon5/mc3bg6p/devono/software/links_v1.8.5/LINKS [v1.8.5]
-f /pylon5/mc3bg6p/devono/alignments/mylu/mylu_bwa/mylu.fa
-s /pylon5/mc3bg6p/devono/assembly/links/mylu/mylu.fof
-m
-d 1000
-k 15
-e 0.1
-l 5
-a 0.3
-t 10
-o 0
-z 500
-b mylu1000
-r
-p 0.001
-x 0The job stopped after about 12 hours. The tail of the .log file suggested that the issue was exceeding my 1.2 Tb memory allocation. There appear to be many, many kmers extracted. This has me wondering if I should be altering something about that above command?
Extracted 7481580 15-mer pairs, from all 5008903 sequences provided in /pylon5/mc3bg6p/devono/assembly/links/mylu/mylu.fof
Extracted 2556027163 15-mer pairs overall. This is the set that will be used for scaffolding
=>Reading sequence contigs (to scaffold), tracking k-mer positions: Wed Jul 4 21:11:58 EDT 2018
Contigs (>= 500 bp) processed k=15:
418slurmstepd: error: Job 3427612 exceeded memory limit (1258291556 > 1258291200), being killed
slurmstepd: error: Exceeded job memory limitThis was to be the first of several iterations, as you've suggested. I wonder if the first pass, given the parameters used above, is the most memory intensive of all? This was my overall intention:
LINKSPATH=/pylon5/mc3bg6p/devono/software/links_v1.8.5/LINKS
INFASTA=/pylon5/mc3bg6p/devono/alignments/mylu/mylu_bwa/mylu.fa
FOFPATH=/pylon5/mc3bg6p/devono/assembly/links/mylu/mylu.fof
$LINKSPATH -f $INFASTA -s $FOFPATH -d 1000 -t 10 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu1000.fa -s $FOFPATH -b mylu2500 -d 2500 -t 5 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu2500.fa -s $FOFPATH -b mylu5000 -d 5000 -t 5 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu5000.fa -s $FOFPATH -b mylu7500 -d 7500 -t 4 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu7500.fa -s $FOFPATH -b mylu10000 -d 10000 -t 4 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu10000.fa -s $FOFPATH -b mylu12500 -d 12500 -t 3 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu12500.fa -s $FOFPATH -b mylu15000 -d 15000 -t 3 -k 15 -l 5 -a 0.3 -r mylu1000.bloom
$LINKSPATH -f mylu15000.fa -s $FOFPATH -b mylu30000 -d 30000 -t 2 -k 15 -l 5 -a 0.3 -r mylu1000.bloomAny suggestions on how to modify parameters to reduce memory usage would be greatly appreciated. My advisor suspects that the issue is the number of unique kmers, and that this is likely due to the high error rate in the Nanopore data itself. Perhaps this speaks to your recommendation of trying to error correct these reads before running LINKS. I'm curious if you've used LINKS with much ONT data before, and if so if there are suggested best practices to follow prior to executing this program.
Thanks René!
warrenlr
commented
6 years ago You'll have to stay under 1 TB.
This is the tricky part, because it varies based on long read quality, coverage and the size of your genome.
To achieve this, you will have to tune -t Yes, you seem to have a lot of ONT reads, hence the large #kmer pairs set extracted 2,556,027,163.
There are guidelines in the readme re: how to curb memory usage. Essentially starting off with a high step size (-t) for short distances.
I recommend using -t 150 with -d 1000 to begin. Monitor the mem. usage and adjust -t accordingly. As -d increases, -t must decrease (otherwise you'll end up with too few pairs for scaffolding over larger kmer distances). Also, you might have to adjust -k. The sweet spot will be somewhere k 15-19 for your 2Gb genome (assuming raw ONT reads).
devonorourke
commented
6 years ago Thanks René,
The latest run using -t 150 and -d 1000 completed; this resulted in merging just one pair of contigs among around 70,000. I'm wondering where to go from here.
The job extracted about 171M 15-mer pairs overall; in addition, there were about 500k unique kmer pairs. What's surprising to me is that there was only one instance in which a pair of contigs could be joined. Of the 72784 contigs, it appears there is only one instance identified in the *assembly_correspondence.tsv file where an overlap exists:
LINKS_scaffold_ID LINKS_contig_ID original_name orientation(f=forward/r=reverse) number_of_links links_ratio gap_or_overlap(-)
scaffold65692 54658 contig_54657 f 8 0 -80Perhaps you have a suggestion of what the next parameter changes should be. In your previous post you mentioned this decision might be informed by memory consumption of the first run. Unfortunately I haven't found an programatic way to monitor the memory usage on this cluster as the run progresses; however I can pull together a summary report for the job. Here's a sample:
JobID JobName Elapsed MaxRSS MinCPU AveCPU AveDiskRead MaxDiskRead AveDiskWrite MaxDiskWrite AveVMSize MaxVMSize ReqMem TotalCPU
------------ ---------- ---------- ---------- ---------- ---------- -------------- ------------ -------------- ------------ ---------- ---------- ---------- ----------
3496821 LUlinks 06:29:11 840Gn 06:29:52
3496821.bat+ batch 06:29:11 6832K 00:00:06 00:00:06 38.42M 38.42M 37.98M 37.98M 390808K 390808K 840Gn 00:06.123
3496821.ext+ extern 06:29:12 1200K 00:00:00 00:00:00 0.00M 0.00M 0 0 107904K 140452K 840Gn 00:00.002
3496821.0 LINKS 06:29:04 233672268K 06:28:20 06:28:20 73874.25M 73874.25M 9046.68M 9046.68M 0 234208308K 840Gn 06:29:46The bloom filter itself isn't particularly large:
Bloom filter specs
elements=677791916
FPR=0.001
size (bits)=9745012672
hash functions=9
...
=>Writing Bloom filter to disk (mylu1.bloom): Fri Jul 6 01:33:23 EDT 2018
Storing filter. Filter is 1218126584 bytes.And the reads themselves shouldn't change in terms of their total memory allocation moving forward (assuming I use the same reads):
=>Reading long reads, building hash table: Fri Jul 6 01:33:24 EDT 2018
Reads processed k=15, dist=1000, offset=0 nt, sliding step=150 nt:
... (a lot of "Reads processed from file lines ) ...
Reads processed from file 46/46, /pylon5/mc3bg6p/devono/reads/nanopore/mylu/albout/run4/fastq_runid_e51c6e4c8f8ae0451a5b38a5123c82bde7f5b36a.fastq:
5008903
Extracted 501246 15-mer pairs, from all 5008903 sequences provided in /pylon5/mc3bg6p/devono/assembly/links/mylu/mylu.fof
Extracted 171359519 15-mer pairs overall. This is the set that will be used for scaffoldingI'm still unclear how to evaluate my known total memory consumption (~230 Mb I think) with the processes that use up the memory (presumably the extracted kmers and ONT data?). I'm curious what recommendations you might make for the next iteration and am wondering how altering the -t and -d parameters will change the next iteration of memory usage. Do they scale linearly?
Thanks for your consideration
warrenlr
commented
6 years ago Glad it's working in your hands. fyi the bloom filter only contains your draft assembly kmers, so it is not expected to be very large.
for monitory resource usage (time and peak mem) you may use: /usr/bin/time -v -o mylu1000.time ./LINKS -f XX blah
Based on what you've told me: -t 10 exceed mem and -t 150 uses 230GB (not MB), you may wish to start -d 1000 (-t 100 or -t 75). You'll just have to try. Then, half -t for each doubling of -d (that's a very rough solution); It depends a lot on your read set size and length distribution.. (you have wish to run abyss-fac on your ONT reads)
note: we rarely start at -d 1000 for large genomes (unless we have low read coverage).
Now that you have it working, it is only a matter of parameterizing (-d/-t first, then -k then -a/-l). I'll be away for 3 weeks starting today; good luck!
Hi Rene, I've tried submitting a LINKS job using a few different parameters (including your recent suggestion by modifying the
-tand -dparameters, but unfortunately my observations remain the same: the program is stuck on building the bloom filter for at least 24 hours. Perhaps it needs even longer before any additional messages are written to the.logfile?I have access to a variety of machines on our computer cluster with anywhere from 128 Gb to 995 Gb of RAM. I've been using the more readily available smaller machines to run LINKS so far, but I don't get a sense that the issue is memory at the moment.
I'm working with a 2 Gb genome and trying to close up an existing assembly with Nanopore reads. I have about 30 Gb of Nanopore data thus far. The two commands I've tried which left me with the same sort of holding pattern on the bloom filter portion were:
(default parameters)
then following your suggestions, as well as modifying the kmer size
Perhaps I need to provide more memory? I've attached a typical .log file that is generated below. Thank you for any information you can provide which might help me get over the bloom filter step.