kevlar-dev / kevlar Goto Github PK
View Code? Open in Web Editor NEWReference-free variant discovery in large eukaryotic genomes
Home Page: https://kevlar.readthedocs.io
License: MIT License
Reference-free variant discovery in large eukaryotic genomes
Home Page: https://kevlar.readthedocs.io
License: MIT License
Currently, interesting k-mer abundances are recomputed by "kevlar filter", but these are not updated in the output. This should be addressed.
On Davis HPC, we've previously had (what we suspect are) issues with the system buffering output from Python scripts. I thought the problem would be restricted to stdout, but it seems to be happening even with handles to actual files on disk.
Flushing the file handle after each write takes more time, so we don't want it to be the default behavior. But we can definitely add it as an option when memory consumption is exorbitant.
The kevlar filter
command already recomputes k-mer abundances and discards k-mers whose correct abundance doesn't satisfy the threshold. Any read with no remaining "interesting" k-mers is discarded.
I was surprised today to find that running the same filter a second time will discard even more reads (no changes for third and subsequent attempts). At first I thought this might be a kind of sweep effect, in which some k-mers only reached the threshold in the first place because they were in reads that had a k-mer with a highly inflated abundance. But that doesn't make sense--if an interesting k-mer is discarded but the read has another, the read should not be discarded.
I need to investigate this further, and make sure the filter is doing the right thing. If it is, kevlar filter
should probably do both rounds of filtering with a single invocation of the command.
Currently kevlar localize
relies on the reference genome having already been processed by bwa index
. It shouldn't be too hard to autodetect whether the index files exist, and if they do not, to invoke bwa index
programmatically from Python.
Here is an example.
>contig4;cc=1
AAGGTGAAAGCAGAAAGTCAGTTGTCCATATGGCTTGGGGAGATAAAGAAGGCCCAGGAAGGCCTCCAGGAAAAGGCTGCCATGTCAGGCAGGACACAGAGGGCAATTGAGGAAAAGTGATTCTTACAAGATGGTGAAGGTGCCATTGTGGGTGTT
CAGGCAGGACACAGAGGGCAATTGAGGAAAA 18 0 0#
GGCAGGACACAGAGGGCAATTGAGGAAAAGT 18 0 0#
CACAGAGGGCAATTGAGGAAAAGTGATTCTT 18 0 0#
ACAGAGGGCAATTGAGGAAAAGTGATTCTTA 17 0 0#
CAGAGGGCAATTGAGGAAAAGTGATTCTTAC 17 0 0#
AGAGGGCAATTGAGGAAAAGTGATTCTTACA 17 0 0#
AGGGCAATTGAGGAAAAGTGATTCTTACAAG 16 0 0#
GCAATTGAGGAAAAGTGATTCTTACAAGATG 17 0 0#
Suggestion per @ctb: see #36 (comment).
The kevlar partition
code constructs a graph of shared interesting k-mers, the connected components of which are distinct "partitions" which are handled separately in subsequent processing. The data structure used to sort the graph and the function used to return the connected components doesn't return CCs in a deterministic order, so running the partition code multiple times on the same input gives different output (well, presumably the same output, but in different orders).
Although it's not scientifically or even technically important, for the sake of evaluation and rapid iterative development it might help to have the partitions reported in a deterministic order.
I've always been thinking in terms of "more memory is always better because lower FPR". This is true, but there's evidence, at least for some computations, that runtime scales with memory consumption. So there's an optimization to be made here: what is the smallest amount of memory I can request that won't lead to big problems with FPR?
If this observation holds and if @ctb and I aren't completely of our rockers, we should actually get a speedup when we parallelize (as in #15). If we run all the parallel processes at once, then the total memory consumption will still be the same, but we would have the option to split it up.
pip install .
vs PYTHONPATH
)Another comment from @drtamermansour: if the samples have different sequencing depths, you might expect to see some k-mers in high abundance in some samples and absent in other samples simply due to chance. Are there circumstances where we need to normalize by depth of sequencing?
Also, if library complexity is lower in one sample than another, the effective depth will be higher even if the same number of reads are generated. Tamer has seen this in RNA-seq samples, is this a problem with WGS libraries? If so, how do you measure it?
Suggestion from @drtamermansour: in a single pass, write count tables to N files (one for each band) in a single pass. Then running kevlar find
in N bands would not require N passes over the entire data set, just loading the count tables from disk N times.
I just wanted to capture this suggestion, I have some concerns and I'm not sure it would yield much benefit.
kevlar find
. Loading from count table files rather than the Fastq files again probably won't make a huge difference in overall runtime.And in any case, this is all optimization: there's still work to do to get reliable results first!
[1] There are ways we could investigate to do this in a streaming fashion, but for now I'm happy with saying we have to do a second pass over the reads. :-)
In addition to picking appropriate memory sizes, there's also the question of picking an appropriate value of k. For the human stuff we're currently limited by khmer's maximum of 32, but that should be solved very soon.
I'm wasting a lot of development time loading and re-loading Fastq data into counting data structures. Having the convenience of supporting Fasta/Fastq files as input is a priority, but I desperately need to take advantage of pre-computed countgraphs/counttables when they are available.
I propose the following behavior, which I think is 1) fairly reasonable from a user's perspective and 2) fairly straightforward to implement.
.cg
or .countgraph
, load the file with khmer.load_countgraph()
..ct
or .counttable
, load the file with khmer.load_counttable()
.ct
and load with ct.consume_seqfile()
.In case one needs to adjust match / mismatch scores or gap open / extend penalties.
docs/
directory of master
branch (see this)In the current assembly code, edges in the "shared interesting k-mers" graph are processed in order of overlap: the largest overlaps first. I should investigate processing instead by node degree: high-degree nodes are less likely to lead to a premature truncation of the contig due to a sequencing error.
cc @ctb
Will migrate once open PR(s) are merged.
The find-novel-kmers.py
script (soon to be kevlar find
) currently supports three types of output.
#
)I propose that enable creation of 3 additional output files, all tabular.
All of the information is already there, but this will make it much more searchable with standard shell commands.
We're seeing some strange contamination in the 12175 SSC trio. This morning we discussed a read filtering strategy. While checking each read for "interesting" k-mers, if any k-mer has an abundance below some threshold \tau
then the entire read should be discarded. This should be fairly easy to implement as an option, and hopefully shouldn't slow down performance too much. If anything, it might improve performance by discarding reads before all k-mers have been checked.
There is a lot of room for improving code documentation with docstrings.
At the moment kevlar is painfully slow. We've been philosophizing for a while now whether this was more likely due to the poor cache locality of khmer's Count-Min Sketch implementation or to slow Fastq parsing/handling, or some other factor or combination of factors. It's time to evaluate this empirically.
I profiled both kevlar count
and kevlar novel
on the new simulated data set I created in #83. It's small enough for kevlar to process in several minutes, but large enough to observe meaningful numbers in the profile. For kevlar novel
, I computed k-mer counts from scratch, not using precomputed counts.
Here are the top 25 functions, sorted by time spent, for each.
# kevlar count
ncalls tottime percall cumtime percall filename:lineno(function)
51921279 61.111 0.000 61.111 0.000 {method 'get' of '_khmer.KHashtable ' objects}
36053659 46.578 0.000 46.578 0.000 {method 'add' of '_khmer.KHashtable ' objects}
3 38.588 12.863 187.651 62.550 counting.py:22(load_sample_seqfile)
887492 29.514 0.000 29.514 0.000 {method 'get_kmers' of '_khmer.KHashtable ' objects}
887495 5.401 0.000 5.401 0.000 __init__.py:97(multi_file_iter_khmer) 887492 2.150 0.000 2.150 0.000 {method 'split' of '_sre.SRE_Pattern' objects}
1774984 1.703 0.000 6.031 0.000 __init__.py:103(clean_subseqs)
887819 1.123 0.000 1.157 0.000 re.py:278(_compile)
887492 0.905 0.000 4.177 0.000 re.py:195(split)
3 0.323 0.108 0.323 0.108 {method 'save' of '_khmer.KHashtable ' objects}
82/54 0.154 0.002 0.209 0.004 {built-in method _imp.create_dynamic}
896021/895669 0.153 0.000 0.153 0.000 {built-in method builtins.len}
257 0.101 0.000 0.101 0.000 {built-in method __new__ of type object at 0x10ffa3080}
522 0.050 0.000 0.050 0.000 {built-in method marshal.loads}
3 0.027 0.009 0.027 0.009 {method 'read' of '_io.BufferedReader' objects}
803/1 0.021 0.000 188.297 188.297 {built-in method builtins.exec}
2761 0.019 0.000 0.019 0.000 {built-in method posix.stat}
613/609 0.016 0.000 0.028 0.000 {built-in method builtins.__build_class__}
243 0.015 0.000 0.015 0.000 {built-in method builtins.compile}
1190 0.015 0.000 0.064 0.000 <frozen importlib._bootstrap_external>:1215(find_spec)
522 0.012 0.000 0.018 0.000 <frozen importlib._bootstrap_external>:816(get_data)
2 0.008 0.004 0.021 0.010 machar.py:116(_do_init)
6226 0.007 0.000 0.019 0.000 <frozen importlib._bootstrap_external>:50(_path_join)
6226 0.007 0.000 0.010 0.000 <frozen importlib._bootstrap_external>:52(<listcomp>)
704 0.007 0.000 0.087 0.000 <frozen importlib._bootstrap>:879(_find_spec)
# kevlar novel
ncalls tottime percall cumtime percall filename:lineno(function)
67446100 80.216 0.000 80.216 0.000 {method 'add' of '_khmer.KHashtable ' objects}
1183526 36.393 0.000 36.393 0.000 {method 'get_kmers' of '_khmer.KHashtable ' objects}
29633082 34.668 0.000 34.668 0.000 {method 'get' of '_khmer.KHashtable ' objects}
3 25.096 8.365 142.916 47.639 ./kevlar/counting.py:22(load_sample_seqfile)
22496672 18.973 0.000 54.691 0.000 ./kevlar/novel.py:23(kmer_is_interesting)
1 14.966 14.966 234.252 234.252 ./kevlar/novel.py:42(main)
887495 5.090 0.000 5.090 0.000 ./kevlar/__init__.py:97(multi_file_iter_khmer)
296035 3.894 0.000 10.359 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/site-packages/screed/fastq.py:14(fastq_iter)
887492 1.877 0.000 1.877 0.000 {method 'split' of '_sre.SRE_Pattern' objects}
1774984 1.387 0.000 4.978 0.000 ./kevlar/__init__.py:103(clean_subseqs)
1184137 1.336 0.000 1.980 0.000 {method 'readline' of '_io.BufferedReader' objects}
7784325 1.175 0.000 1.175 0.000 {method 'append' of 'list' objects}
1183855 1.153 0.000 1.187 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/re.py:278(_compile)
1184137 0.998 0.000 3.905 0.000 /usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/lib/python3.5/gzip.py:370(readline)
887492 0.753 0.000 3.462 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/re.py:195(split)
1184138 0.635 0.000 0.927 0.000 /usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/lib/python3.5/_compression.py:12(_check_not_closed)
296037 0.544 0.000 0.544 0.000 {method 'search' of '_sre.SRE_Pattern' objects}
1184137 0.544 0.000 0.968 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/site-packages/screed/utils.py:4(to_str)
296034 0.506 0.000 0.625 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/site-packages/screed/screedRecord.py:23(__init__)
1184179 0.424 0.000 0.424 0.000 {method 'decode' of 'bytes' objects}
1186054 0.366 0.000 0.366 0.000 {method 'startswith' of 'str' objects}
2429318/2428966 0.359 0.000 0.359 0.000 {built-in method builtins.len}
8852 0.317 0.000 0.317 0.000 {method 'decompress' of 'zlib.Decompress' objects}
1184140 0.292 0.000 0.292 0.000 /usr/local/Cellar/python3/3.5.2_3/Frameworks/Python.framework/Versions/3.5/lib/python3.5/gzip.py:296(closed)
296037 0.289 0.000 1.154 0.000 /Users/standage/Projects/kevlar/env/lib/python3.5/re.py:170(search)
In both cases, the add
and get
functions for incrementing and querying k-mer counts from the Count-Min sketch dominate the runtime, with the get_kmers
and re.split
functions as heavy hitters as well. The latter two have to do with the fact that the khmer's bulk loading functions don't support the kind of preprocessing we need for kevlar, so I'm doing it in Python. This incurs overhead in sending the data from C++ to Python objects, and then doing the processing in Python.
As far as priorities, I'm not sure there's much we can do about the Count-Min sketch implementation. We could try to implement buffering (collect N≈1e4 add operations before actually incrementing the tables) but honestly once the CQF-based counter is integrated we may already get much better performance from that. For the sequence loading, I'll continue to lobby for better support for multiple pre-processing strategies with khmer bulk Fastq loading code.
A while ago I replaced linear path assembly with junction count assembly in the kevlar collect
step. This eliminated problems we previously had with cyclical structure in the assembly graph, but it doesn't seem to have done much about truncation at sequencing errors. I've been playing around with khmer's simple labeled assembler, but this will take some effort to get working.
@ctb suggested investigating abundance-based trimming / spur removal. Probably a medium-term priority.
In the mean time, I assembled the reads described in dib-lab/khmer#1657 with velvet (using vanilla parameters) and got almost the entire consensus sequence of the reads' multiple sequence alignment. If we can get partitioning working well, one preliminary approach could be to simply delegate assembly of each distinct CC to a third-party assembler.
Our current workflow and scripts support the idea of splitting data by chromosome, which is intended to:
The problem is that we want this to be a reference-free method, but we're relying increasingly on the reference and read mappings against the reference. The current criteria in kevlar dump
are:
The first two points are fine, and when mapping data happens to be available we can leverage it. But we keep on running into problems with the second two points.
Today we discussed the following alternative approach.
A couple of considerations.
The latter two can eventually be refactored into create_countgraph and create_nodegraph.
Currently, assembled contigs are aligned against a small cutout of the reference genome for variant calling. We need to be able to convert the coordinates relative to the reference back to the global coordinates relative to the entire chromosome.
See tests/test_rcgrep.py
for an example.
The kevlar find
command prints a summary of the results at the end of the program execution, like so.
processed 629000000 reads...
processed 630000000 reads...
processed 631000000 reads...
[kevlar::find] Found 2374 novel kmers in 40279 reads, 2366 linear paths
However, there's a discrepancy between the number of reads reported and the number of reads found in the output file. The output file is the reads in Fastq format, adorned with additional lines containing the novel k-mers (and their abundances) with a trailing #
character.
@ERR894723.92068452
GAATAGAATGGAATGGAAAGAATTGGAATGGAATGGAATCGAATGGAATGGAATGGAATGGAATGGAATGGAATCAACCCGAGTGCAGGGGAATGTAATGGATCGGAATGCAGTGGAATGGAATCA
+
70<0070B0<<0<00<<BB<0<0BB<<B0<B<70B<BBB7<0B<B<<B0B<<BBBBBBBBBB<<BBBBBBBBBBBBB<B7<B<BBBBBBBBBBB<BBBBBBBB7BBBBBBBB<BBBBBB<<B<<<<
CAGGGGAATGTAATGGATCGGAATGCAGTGG 26 0 0#
AGGGGAATGTAATGGATCGGAATGCAGTGGA 25 0 0#
GGGGAATGTAATGGATCGGAATGCAGTGGAA 26 0 0#
GGGAATGTAATGGATCGGAATGCAGTGGAAT 27 0 0#
GTAATGGATCGGAATGCAGTGGAATGGAATC 27 0 0#
TAATGGATCGGAATGCAGTGGAATGGAATCA 26 0 0#
@ERR899709.54045623
TGGAATGGAAAGAATTGGAATGGAATGGAATCGAATGGAATGGAATGGAATGGAATGGAATGGAATCAACCCGAGTGCAGGGGAATGTAATGGATCGGAATGCAGTGGAATGGAATCATCCGGAAT
+
<B<BBB<BBBBBBB<B<BBBBBBBBBBBBBB<7BB<BBBBBBB7B<<BB<<BBBB<BBBB0<<<BBB<B<B<7BB<B<BBBB<<B0B0<0BBBB0<7BBBBBBB0BBBB0BBB<B00000000<<<
CAGGGGAATGTAATGGATCGGAATGCAGTGG 26 0 0#
AGGGGAATGTAATGGATCGGAATGCAGTGGA 25 0 0#
Removing the extra metadata to get strict Fastq format is simple: grep -v '#$' file.txt > file.fq
. We can then determine the number of reads in the Fastq file by computing the line count and dividing by 4.
$ wc -l NA19240.novel.fq
161688 NA19240.novel.fq
$ python -c 'print(161688 / 4)'
40422
The kevlar dump
command implements a data reduction step that discards any reads that match the reference genome perfectly. After this step, however, kevlar makes no attempt to determine whether a k-mer is in the reference genome. As a result, we are seeing some false positives reported in cases where certain k-mers are "not present" in the parents (because all overlapping reads match the genome perfectly and are thus discarded before counting) but "high abundance" in the proband (a sufficient number of reads overlapping the k-mer have errors elsewhere in the read and are thus NOT discarded).
For novel variant discovery, we're not really interested in k-mers that are in the reference genome, and we need to update our counting strategy to reflect this. There are potential implications for memory consumption to consider and investigate.
When case and controls are loaded into counttables for kevlar find
, the script should calculate FPR. This code already exists in khmer: in the short term, it's probably quickest and easiest to duplicate the code, but in the future it is probably worth refactoring the khmer code to make it reusable.
I had all but condemned kevlar collect
, which uses khmer's linear path assembly foo, to deprecation and phasing out. But it turns out that it seems to work just fine on kevlar-partitioned and abundance-trimmed data. I need to revisit this and perhaps integrated it with the existing kevlar assemble
command.
Currently, the kevlar filter
command has the option to output reads partitioned by shared novel k-mers. The kevlar assemble
command is stricter: it builds its graph not only with shared novel k-mers but also enforces that the entire overlaps must be identical. Two proposals.
kevlar partition
commandCurrently kevlar find
reports:
In an ideal world, each variant would be associated with a single contig/path that is long enough to determine its position in the genome uniquely. Practically, that's not happening.
Consider the following example, newly introduced in #16. The case has a 5bp deletion relative to the controls. The kevlar find
script found 11 "interesting" k-mers (k=13) in 16 reads, and these 11 k-mers pulled out 4 unique linear paths from the assembly graph.
GGGTCCCTATTGACCTCTTTACCA
GGGTCCCTATTGACC
CTATTGACCTCTTTACCA
TGGGCTCGGGGTCCCTATTGACC
A couple of things to note.
kevlar find
shows only two linear paths).This second bullet point is what perplexes me now. All along when I've observed this in the human data, I assumed it was the result of sequencing errors introducing structure in the assembly graph and truncating linear paths. But for this example, I simulated reads with no errors or mutations. Why do either of these cases occur?
Currently, kevlar novel
loads each sample into a dedicated Bloom filter, each of size M
, while kevlar count
loads the first sample into a BF of size M
and all subsequent samples into BFs of size M/x
for some factor x
. This latter behavior should be an option in kevlar count
.
Nothing major, but most of the commands output a list of read or contig sequences. The idea would be to convert the existing "main" methods into a generator that yields records instead of printing them, and then use a minimal wrapper around that as the "true" main method to actually print the output.
This would make it trivial to programmatically string together subcommands.
Currently, kevlar find
only supports pre-computed count graphs. It would be nice if there was also an option to compute counts on the fly from Fastq input. This will clutter up the CLI a bit, especially if there's a reason to support different table sizes for different input files (although I don't think there is).
The kevlar find
procedure already has support for an arbitrary number of controls. It should be trivial to introduce support for an arbitrary number of cases. For situations where we have many cases and controls and memory is a limiting factor, we'll have to use k-mer banding and multiple passes over the data.
The LinearPathAssembler kevlar has been using thus far to retrieve contigs from the assembly graph doesn’t handle non-branching cycles in the graph very well. As a stopgap solution, I added an option to kevlar collect
to specify a list of k-mers to ignore when retrieving contigs from the graph. To find out which k-mers were problematic, I had to run kevlar collect
a few dozen times with debugging print statements to find out which k-mers were causing it to choke. Luckily kevlar collect
is pretty quick and this didn’t take a huge amount of time, even though it was tedious. For my latest run (results from 32 kevlar find
runs with k-mer banding), there were 12 problematic k-mers. Adding these to the kevlar collect --ignore
flag leads to a successful run.
Since all of these problematic k-mers were AT-rich and low-complexity, my first thought for a more robust solution was to filter out low-complexity novel k-mers. I poked around with this a bit, but I’m concerned that it will lack specificity, and it turned out complexity scores are not as trivial to compute as I initially thought.
And then of course I thought of the fact that we’ve been talking about using a less naive assembly approach for a long time and tried substituting the LinearPathAssembler with the JunctionCountAssembler. This was a straightforward substitution, and handled all of the previously problematic k-mers with finesse. I still need to evaluate the contigs assembled by the two different assemblers, but this seems like the best long-term solution, especially if it gives us longer contigs.
The kevlar commands don't use pairing information, so I haven't really worried to much about it before. But some parts of kevlar do rely on reads having unique IDs, and kevlar dump
currently prints the same ID for reads belonging to the same pair. This is a straightforward fix and should be done.
Related to #52, making kevlar
an entry point rather than a script and reorganizing the argument parsers a bit will not only make them easier to test (increasing coverage), but will also make it easier to replace StringIO with capsys for capturing test output.
Since pre-processing input BAM files prior to analysis isn't an option, we should be able to handle this on a partition-by-partition basis. We don't have access to both ends of a read, so duplicate sequences MIGHT not be PCR duplicates, but we can deduplicate anyway (make it an option to disable deduplication).
It's heavy and kludgy for CI (which means it's probably a pain on a lot of HPC as well). Need Python-side support for BAM parsing, but something in Pure python would be preferred.
The cythonized C file is distributed with the code base, so the user shouldn't need to have Cython installed to compile and install and run kevlar.
When using the banding approach, there is a discrepancy between the reported abundance and actual abundance in some of the k-mers reported as novel by kevlar find
. This is particularly puzzling because the overall false positive rate is fairly low (0.14), but the k-mers in question have a reported abundance of >10 but an actual abundance of 1 (in all of the cases I’ve had time to fully investigate). Even more puzzling is that when I load the data into memory later (with an interactive python session) I cannot reproduce the error.
Currently the variantset.collapse
method takes a very naive approach to collapsing contigs that are subsequences of other contigs. The contigs are first sorted by length, and then added to a new set of unique contigs only if they are not a subsequence of any contig already in the unique set.
unique_contigs = set()
for contig in sorted(contigs, key=len, reverse=True):
contigrc = kevlar.revcom(contig)
keep = True
for uc in unique_configs:
if contig in uc or contigrc in uc:
keep = False
break
if keep:
unique_contigs.add(contig)
I’ve tried both a set (fast membership queries, slow iteration) and list (fast iteration, slow membership queries) as the container variable unique_contigs
, but the naive O(N^2)
algorithm washes out any small benefit of changing the data structure.
I explored an alternative approach using a suffix tree and got a huge speedup (≈15min to ≈20 seconds).
unique_contigs = set()
suffix_tree = trie()
for contig in sorted(contigs, key=len, reverse=True):
contigrc = kevlar.revcom(contig)
query = suffix_tree.with_prefix(contig) + \
suffix_tree.with_prefix(contigrc)
if query == []:
unique_contigs.add(contig)
for i in range(len(contig)):
suffix_tree[contig[i:]] = 1
Two concerns with this approach.
trie
implementation). I couldn’t find any Pure python suffix tree implementations: probably a lost cause with respect to memory overhead.From @ctb:
In update/filter_below_abund branch of khmer, I’ve fixed a few sandbox scripts and implemented dedup.py. I then ran the attached commands.
Briefly, they —
de-duplicate your reads based on name
load them into a counting table
slice out reads between coverage 2 and 120 (into .selected)
trim off low-abundance k-mers
partition
extract partitions
build a cDBG for group 12 of the partitions, which contains 6 partitions.
python dedup.py all-fq > all-fq.dedup
../scripts/abundance-dist.py all-fq.dedup all-fq.dedup.hist -k 31
../scripts/abundance-dist-single.py -k 31 -M 1e8 all-fq.dedup all-fq.dedup.hist
../scripts/load-into-counting.py -k 31 -M 1e8 dedup.kh all-fq.dedup
../sandbox/calc-median-distribution.py dedup.kh all-fq.dedup all-fq.dedup.medhist
../sandbox/slice-reads-by-coverage.py -m 2 -M 120 dedup.kh all-fq.dedup all-fq.dedup.selected
../scripts/trim-low-abund.py -k 21 -M 1e8 all-fq.dedup.selected
../scripts/abundance-dist-single.py -k 21 all-fq.dedup.selected.abundtrim all-fq.dedup.selected.hist -M 1e8 -s
../scripts/do-partition.py -k 31 -M 1e9 y2 all-fq.dedup.selected.abundtrim
../scripts/extract-partitions.py -X 5000 foo all-fq.dedup.selected.abundtrim.part
../sandbox/extract-compact-dbg.py -x 1e9 -k 31 foo.group0012.fq -o foo.group0012.fq.gml
grep ^@ foo.group0012.fq | cut -f2 | sort | uniq
When kevlar find
is run in banding mode, some reads will show up multiple times across the various outputs ("interesting" k-mers from each read will be distributed across different bands). Currently, kevlar collect
handles:
The first 3 filtering steps should be decoupled from the graph / assembly steps in kevlar. Maybe I'll call it kevlar recalibrate
, because NGS.
...for reproducibility's sake, yo.
There are a variety of parameters in the variant discovery problem that we should explore at some point.
In #16 I started down this road, and began to integrate this into the automated test suite. That was a step in the wrong direction. The test suite should guarantee correctness of the basic behavior and test for bug regressions. The type of exploratory analysis needed here should be handled separately.
Already supported for khmer consume_
functions, but not yet for augmented Fastq files.
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.