thegenemyers / fastk Goto Github PK
View Code? Open in Web Editor NEWA fast K-mer counter for high-fidelity shotgun datasets
License: Other
A fast K-mer counter for high-fidelity shotgun datasets
License: Other
Hello,
Thanks for FastK, it is truly useful.
I'm looking at kmer coverage of contigs of a small assembly. It means that if I want to get kmer coverage/histogram for each single contig I need to create 1 fasta file per contig and apply FastK to it. It is quite involved.
I do trivial parallelisation over all files. My pipeline stop randomly as a consequence of a segfault on random contigs. Running incriminated step by itself outside of the pipeline does not allow to reproduce the issue. Restarting the pipeline from scratch does make it segfault again but not on the same files. My intuition is that there might be an issue with multiple FastK instance writing in the same temporary folder?
Here is an example of error message:
/bin/bash: line 1: 193787 Segmentation fault (core dumped) FastK -t1 -T2 seqs/folder43/contig_12.fa -P'seqs/folder43'
I am having another similar problem with random segfault. This time with Logex with
Logex -H 'result=A&.B' sample.ktab contigs.ktab
In all examples I've looked at, segfault happened when contigs.ktab was empty (well formed table with 0 kmer). Though running the same command line outside of the pipeline works without issue and indeed produce a working .hist file (even though with no kmer).
This second issue is less problematic in the sense that I can just pre-filter for empty contigs.ktab.
Additionally here are a few miscellaneous issues I encountered. I'm mostly puzzled by the first one:
If it's useful, I'm working on Ubuntu 16.04.7 LTS with gcc version 9.4.0.
Best,
Seb
Hello Gene,
Lines 648 to 650 in 25277f0
should line 648-650 of libfastk.c be
c = s1[0];
d = s2[0];
e = s3[0];
Chenxi
I want to understand how kmers are counted in FastK and how that affects the totals in merquryFK calculations.
Why do the total values in completeness.stats and qv files differ so much? What do they represent and how they relate to each other? I run merquryfk with a single genome assembled using Pacbio HiFi and HiC data, and run against an Illumina kmer dataset.
# mMelMel1_T1.qv
Assembly No Support Total Error % QV
GCA_922984935.2.subset 6005 7999890 0.0024 46.2
# mMelMel1_T1.completeness.stats
Assembly Region Found Total % Covered
GCA_922984935.2.subset all 2268391877 2268397787 100.00
From Merqury, marbl/merqury#84
The Total in QV are kmers that are 'present' in the assembly. So if there is one specific kmer found 3 times in the assembly, but never in the reads, it is counted as 3 error kmers (no suppurt). The 3 error kmers are part of the Total.
The Total in completeness are distinct solid kmers in the reads. In other words, a kmer that is present over a certain frequency in the reads is counted as one kmer. I forgot how exactly the Total is computed in MerquryFK completeness. It's likely that it is only filtering out kmers with frequency of 1, which is the default in FastK? Might be a good question for Gene.
We expected the opposite because the total for QC is ~8M whereas the total for Completeness is ~2.2B.
Synopsis: I think the suffix[] table at line 138 in io.c has a typo in the last item in the second row.
Details:
Looking to use this package to quickly identify half-covered scaffolds (or scaffold segments) in an assembly, having seen the two presentations in VGP conference calls.
Working with a mar/31 clone of the repo.
For first test I gave it a small random apple.fa file. When I tried 'FastK apple.fa', it reported "FastK: Cannot open apple.fa as a .cram|[bs]am|f{ast}[aq][.gz]|db|dam file."
That report comes from Fetch_File() in io.c. Examining that, well, it's not completely clear to me what the purpose of the suffix[] and extend[] tables are (maybe it forgives the user for giving an incomplete file extension?). But given that the loop that attempts the different extensions, in order, records the first success in i, and then assigns ftype as FASTQ or FASTA based on whether i is odd or even (and >=5), the intent appears to be that the table should alternate between fastq and fasta extensions. But at the end of the second line of suffix[] we see ".fq" twice. I infer the second of those should be ".fa", and when I make that change FastK is then able to open apple.fa without complaint.
Is there any method to convert "jellyfish count"/"KMC count" or "jellyfish dump" format to FastK db format?
As I couldn't see the FASTK's Docker container, I created it by myself. Feel free to take the file if you like.
https://gist.github.com/junaruga/dd71378add9fb8e840dad65bb050bf0d
https://quay.io/repository/junaruga/garg-fastk?tab=info
Note we are using the container at https://github.com/GargGroup/BioDivGenomics/tree/main/FASTK as a temporary workaround. After you prepare the Docker container official version, we would like to use it.
As you may know, if you create the bioconda package, then the Docker container is automatically created at https://biocontainers.pro .
E.g.
https://anaconda.org/bioconda/hifiasm
https://biocontainers.pro/tools/hifiasm
https://quay.io/repository/biocontainers/hifiasm?tab=info
Hi Gene,
I wonder if it would be possible for Logex to accept a modimizer for the output? Like say Logex 'reduced=A%51==0' source<.ktab>? I would like to get modimized homozygous kmers to use alongside het kmer pairs.
Best,
Haynes
Hi Gene,
Just digging into this. Looks great so far. Can you make Haplex output a ktab file as well as a text file? I'm wanting to get a set of het kmer pairs and then run those with some modimized homozygous kmers run through the profiler so that I can better categorize het kmer pairs. Or if there is a way to go from the text output to a .ktab that would be good as well.
Thanks,
Haynes
Hello Gene,
Logex and Fastmerge assume the ibyte
field is always 3. I ran into a very small genome (~40 Mb), which gave me 2. The histograms from Logex and Fastmerge are correct. The Kmer tables, however, are junks.
(1) Line 1026 and 1044 should be something like bst+ibyte
here instead of bst+3
?
Line 1026 in d694359
Line 1044 in d694359
(2) Line L1028 and L1050 should be something like this?
x = 0; if (ibyte==3) x = ((int) bst[0] << 16) | (bst[1] << 8) | bst[2]; else if(ibyte==2) x = ((int) bst[0] << 8) | bst[2]; else if(ibyte==1) x = bst[0];
Line 1028 in d694359
Line 1050 in d694359
(3) You write the minval
as 1. Do you think if it is worth updating it after Logex?
Line 1391 in d694359
(4) another minor point is do you think if it makes sense to grant the group read permission to files generated by Logex?
Fastmerge is pretty much the same as Logex.
Best,
Chenxi
Hi @thegenemyers,
Thanks for developting this tool, I like it a lot.
I wanted to profile a ~20 Gb gzipped fastq file using:
FastK -v -p -N../Test/FOO -P./ My_file.fastq.gz
And it goes successfully through the first three stages but then in Phase 4: Merging Profile Fragments
It breaks everytime and displays,
Phase 4 (-p option): Merging Profile Fragments
0%
FastK: Cannot open external file /current/directory//FOO.25.P0.0 in /current/directory/
*** Error Exit 1 ***
sh: error while loading shared libraries: libc.so.6: cannot open shared object file: Error 24
Could not clean up !
EDIT: Simplified path names
So it seems the -N
argument is not working correctly inb phase 4 when using the -p
profiling mode. I can confirm that the FOO profile files are partially complete in ../Test/
. I also tested FastK -v -p -N/current/directory/FOO -P/current/directory/ My_file.fastq.gz
Now the Error 24
suggests there a too many files open. But I raised my ulimit to the hard limit which is >1M, and monitoring the files attached to the process does not indicate large numbers at all. Is this a bug or am I using the command wrong?
Kind regards,
Ronald
Granted, this is a stupid mistake on my part. I created a file named sorang.telomeric.fasta but in fact it was a gzipped file. The first two byes of the file are the gzip magic number: 1f 8b.
Running FastK -v -k40 -t1 data/sorang.telomeric.fasta
produced this output:
Partitioning 1 .fasta file into 4 parts
Determining minimizer scheme & partition for sorang.telomeric
Estimate 155.896M 40-mers
Handling data in a single block
Segmentation fault
Clearly, the fault is mine. But it did take me a while to discover the error of my ways.
I wonder if some simple test could be done early on to validate that what the user claims is fasta is really fasta. The fact that the first character (1f) is outside the printable ascii range, is strong evidence. Looking at the code (but not completely understanding it), it looks like such a test could be added to the FK1 case in fast_nearest without (probably) adding much runtime overhead.
(This doesn't seem like it needs to be addressed in the short term, if at all).
I tried "FastK orange.fa.gz" where orange.fa.gz was a symlink. The result is
gunzip: ./orange.fa.gz is not a regular file
FastK: Cannot get stats for ./orange.fa
Apparently gunzip doesn't like symlinks, so it fails to create the unzipped file. The downstream code in FastK, I guess, doesn't notice that gunzip failed, but a later sanity check notices the unzipped file doesn't exist.
This would be an issue for the use case where the user has read access to a shared directory of gzipped read data or assemblies, but doesn't have write access. They can't give FastK a path to the original (because, I presume, it would try to write the unzipped file in that directory). Traditionally a symlink would be the 'right' solution, to avoid wasted disk space. But perhaps disk deduplication technology makes this less of an issue?
I suspect this is only an issue for gzip'd files. I assume for the other compressed formats you are able to decompress on-the-fly and don't need to write an uncompressed file.
Bob H
synopsis: When running FastK with -k set 32 or less, I'm seeing segfaults.
Specifically, I've seen this with one particular input file (the blue whale assembly) and k in {20,26,31,32}.
Details:
I was trying to count 31-mers in this VGP blue whale assembly:
https://s3.amazonaws.com/genomeark/species/Balaenoptera_musculus/mBalMus1/assembly_curated/mBalMus1.pri.cur.20200528.fasta.gz
My current directory had mBalMus1.pri.cur.20200528.fasta.gz as a symlink to some other directory. Then I ran
FastK -v -k31 mBalMus1.pri.cur.20200528.fasta.gz
and got this output:
Gzipped file mBalMus1.pri.cur.20200528.fasta.gz being temporarily uncompressed
Partitioning 1 .fasta file into 4 parts
Determining minimizer scheme & partition for mBalMus1.pri.cur.20200528
Estimate 2.375G 31-mers
Dividing data into 2 blocks
Using 5-minimizers with 1024 core prefixes
Phase 1: Partitioning K-mers into 8 Super-mer Files
There are 105 reads totalling 2,374,852,541 bps
Part: 31-mer super-mers ave. length
0: 1,151,752,026 74,926,939 15.4
1: 1,197,133,971 112,483,639 10.6
Sum: 2,348,885,997 187,410,578 12.5
Range 1,151,752,026 - 1,197,133,971 (3.86%)
Resources for phase: 1:27.610u 5.846s 1:01.816w 151.2%
Phase 2: Sorting & Counting K-mers in 2 blocks
Processing block 2: Sorting super-mers **Segmentation fault**
I'm using commit 7cebc7d, from a few hours ago.
I followed the example 2 for hpc operations as I got an not enough disk space error
when running fastk in single instance.
I completed the commandsfastk
and fastmerge
and end up with 25 slices. For each slice I have two files, e.g. slice10.hist
and slice10.ktab
. Now I am trying to apply the command fastcat
.
However, it throws the following error:
Fastcat -vht full slice1 slice2
Fastcat: Target name cannot have a .hist, .ktab, or .prof suffi
I tried renaming files to remove the suffix, but it doesnt make a difference and the error persists.
I also tried:
Fastcat -vt full slice1.ktab slice2.ktab
Fastcat: Target name cannot have a .hist, .ktab, or .prof suffi
What am I doing wrong?
Could you tag a release so it's easier to version control fastk?
For linked reads we want to trim 23 bases from R1 but not from R2. Right now I'm making two calls and then doing a Logex 'result = (A |+ B)' R1 R2.
This would be fine except for min counts on 2 files vs one file. Making a trim of 23 for both is okay except you do lose a fair amount of data. I can make this change and make a pull request if you are willing to take pull requests. Or i can just make this change in my forked version.
Would it make sense for profex to be able to output all profiles? I am wanting to stream it and save all non-0 entries.
Thanks for creating yet another incredibly useful looking bit of a software :) I have been searching for a program with all this functionality for a while but I have run into this issue. I have created this as a new issue as I believe it is separate to the file write to /tmp problem raised previously. I am getting corrupted memory faults when try to profile the following sequences downloadable as:
wget https://ebitutorial.s3.climb.ac.uk/graph5.fasta
Running the command:
FastK -k51 graph5.fasta -v -T1 -p
I get:
Partitioning 1 .fasta file into 1 parts
Determining minimizer scheme & partition for graph5
Estimate 1.703M 51-mers
Handling data in a single block
Using 5-minimizers with 1024 core prefixes
Phase 1: Partitioning K-mers into 1 Super-mer Files
There are 38 reads totalling 1,705,389 bps
Part: 51-mer super-mers ave. length
0: 1,703,489 81,679 20.9
Sum: 1,703,489 81,679 20.9
Range 1,703,489 - 1,703,489 (0.00%)
Resources for phase: 0.076u 0.016s 0.108w 84.7%
Phase 2: Sorting & Counting K-mers in 1 blocks
Part: wgt'd k-mers savings
0: 988,833 1.7
All: 988,833 1.7
Resources for phase: 0.312u 0.028s 0.368w 92.1%
Phase 4 (-p option): Merging Profile Fragments
That may be the /tmp dir problem. If I point it at a drive with a large amount (1 Tbp of free space) as temp dir with -P I get:
Partitioning 1 .fasta file into 1 parts
Determining minimizer scheme & partition for graph5
Estimate 1.703M 51-mers
Handling data in a single block
Using 5-minimizers with 1024 core prefixes
Phase 1: Partitioning K-mers into 1 Super-mer Files
There are 38 reads totalling 1,705,389 bps
Part: 51-mer super-mers ave. length
0: 1,703,489 81,679 20.9
Sum: 1,703,489 81,679 20.9
Range 1,703,489 - 1,703,489 (0.00%)
Resources for phase: 0.088u 0.008s 0.356w 26.9%
Phase 2: Sorting & Counting K-mers in 1 blocks
Part: wgt'd k-mers savings
0: 988,833 1.7
All: 988,833 1.7
Resources for phase: 0.400u 0.076s 0.551w 86.4%
Phase 4 (-p option): Merging Profile Fragments
Sorry for the long post but I have though the above may be helpful. This is compiled with gcc:
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609
Dear Gene, thanks for making FASTK!
Would you consider making FASTK available also via conda?
I tried to compile it on our cluster, run into some errors and I bet I am not the only one. Having it on conda would make the life of loads of people easier.
Synopsis: I have a contrived example with a short genome and some perfect reads sampled from that genome. I'm using the relative profile feature to get counts of read kmers across the genome. The results are mostly 1s with a few short stretches of higher counts. The 1s conflict with the result I expect.
Details:
The example can be found at https://github.com/rsharris/fake_reads_2021. apple.fa is a random genome of length 3K. orange.fa is 80 fake reads, each 300 bp, sampled from apple with no sequencing error, at random positions and with random orientation. Average depth โ8X.
Running this
FastK -k40 -t1 apple.fa
FastK -k40 -t1 orange.fa
FastK -k40 -t1 -p:orange.ktab apple.fa
Profex apple.prof 1 > orange_on_apple.profile
I expected to get a count profile that roughly matched the coverage profile I would get from piling up orange to apple alignments (realizing that the 40-mers absent at the end of the reads would reduce avg count to โ7X).
(orange_on_apple.profile is in that same repo)
Instead, I get mostly 1s with occasional short bursts of higher numbers. I'm at a loss to understand this. Am I missing something in how I am using the tools? I know the motivation for profiles was to profile along reads, not along genomes, so maybe there's a difference I'm not understanding.
I wonder if the fact that the kmer counts in apple itself are all 1 is relevant.
FWIW I'm using a fetch of the repo from about an hour ago.
Hello,
On the latest master f18a4e6 , I was trying to test the FastK
with a small test data. Then I got the following Segmentation fault. Do you have any idea about how to fix it? Thank you.
$ pwd
/home/jaruga/git/FASTK
$ gcc --version
gcc (GCC) 12.2.1 20220819 (Red Hat 12.2.1-1)
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ make
$ cat fastk_source.fasta
>FastK input source
ACTG
$ ./FastK -v fastk_source.fasta
Partitioning 1 .fasta file into 4 parts
File is so small will use only 1 thread
Determining minimizer scheme & partition for fastk_source
FastK: Warming Sequences are on average smaller than 1.5x k-mer size!
Estimate -0.036K 40-mers
Handling data in a single block
Using 5-minimizers with 1024 core prefixes
Phase 1: Partitioning K-mers into 1 Super-mer Files
There are 1 reads totalling 4 bps
Part: 40-mer super-mers ave. length
0: 0 0 NA
Sum: 0 0 NA
Range 0 - 0
Resources for phase: 0.000u 0.001s 0.001w 127.8%
Phase 2: Sorting & Counting K-mers in 1 blocks
Processing block 1: Sorting super-mers Segmentation fault (core dumped)
Hi, @thegenemyers
Can fastk use the pigz
for parallel gzipping large fastq?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.