torognes / vsearch Goto Github PK

View Code? Open in Web Editor NEW

631.0 631.0 122.0 4.99 MB

Versatile open-source tool for microbiome analysis

License: Other

Shell 0.01% C++ 87.00% C 12.48% Makefile 0.26% M4 0.23% Dockerfile 0.03%

amplicon bioinformatics chimera clustering fasta fastq metabarcoding metagenomics microbiome search sequence-alignment

vsearch's People

Contributors

Stargazers

Watchers

Forkers

hobrien wasade prehensilecode fw1121 pschloss heianxing wwood mothur radiantgenomics ryneches liupfskygre ajwije astanabe giriartes dnieuw xtmgah txje tarah28 anitaannamale sdwfrost baiyuanxiang gmdzy2010 zhaolab2019 genomewalker joe-wan tong-chen sunhuaibo loupxing rajaldebnath colinbrislawn ishwarvh f-heeger hugocode alenzhao inambioinfo kilaza vpylro longbow0 healthvivo n8clark jvkersch plantarum mys721tx tw7649116 jianshu93 jonahventures bioslad lynngoodnight z-shiyi 672234993 aerijman davidealbanese peterjc ravinpoudel ditag blankenberg pipoca16 drmoca liuyangying shenweima jlljlong angrybee qrb prakashraaz stjordanis jungleblack007 ghc4 yaohuatj flywind2 joey10086 phymucs zhaoxia413 robmaz bellaqiao minghao2016 deyi-wangg liaochenlanruo jpetittowpi dengdi0929 joshualiuxu jhuanglabtools nielsonj nvt-1009 sreekanth370 besykes jaimeortiz-david ujjwalsh liuchao-nj ray-chun jianghexiliu snamper wangdi2014 feigeliudan01 fanchen26222 ayanlj martin-hartmann 25280841 dreamerjylx 00mjk schaudge

vsearch's Issues

API

Add a (Cython) API for the core functions of vsearch.

Fix values of some userfields

qlo, qhi, tlo, thi, raw should have meaningful values.

Fix number of columns in blast6out output for non-matching queries

Performance comparison to usearch

Compare sensitivity, precision, speed and memory usage

maxuniquesize option for dereplication

There is a minuniquesize option:

--minuniquesize positive integer
discard sequences with an abundance value smaller than integer.

Maybe the opposite option (maxuniquesize) would be useful?

Hello, very excited to see this project. If you are taking requests for usearch features, I would personally find the allpairs_global option useful. See http://www.drive5.com/usearch/manual/allpairs_global.html Thanks!

do not wrap alignments when rowlen=0

Consistency with output_no_hits and uc_allhits

Make output to alnout, blast6out, userout, uc etc consistent for non-matching queries with output_no_hits and uc_allhits.

Documentation

A man page and other documentation should be created. As for now, the online USEARCH documentation is useful, but VSEARCH should not be dependent on it.

Automatic regression testing

Convert array of top kmer hits into a min heap priority queue

The topscore_insert function in search.cc should be modified to improve performance. This may be significant at least with long lists.

Look at this blog for inspiration.

Improve selection of unique kmers in query

Should T and U be considered as identical during dereplication?

When aligning sequences during searching and clustering, T and U are considered identical, regardless of
their case. During dereplication, T and U are treated as different. Should we harmonize that behavior? T and U treated as identical in every situation?

Sequence masking

When searching a database against itself, sequence labels are not truncated correctly

On a database of sequences with labels structured as such: >sequenceID taxonomy, it seems that the query receives the full fasta header while the target is truncated correctly when using --self. That leads to incorrect results where a query is aligned with itself. The option --selfid yields the correct result (no alignment of a query against itself).

${VSEARCH} \
    --vsearch_global "${SUBJECTS}" \
    --db "${SUBJECTS}" \
    --self \
    --id 0.5 \
    --blast6out results.blast6

head -n 3 results.blast6

chout1 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Suessiales|Suessiaceae|Gymnodinium|simplex|simplex_X   chout1  100.0   162 0   0   1   162 1   162 *   *
chout2 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Suessiales|Suessiaceae|Polarella|glacialis|glacialis_X chout2  100.0   234 0   0   1   234 1   234 *   *
chout20 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Peridiniales|Peridiniales_X|Heterocapsa|sp|sp_X   chout20 100.0   229 0   0   1   229 1   229 *   *


${VSEARCH} \
    --vsearch_global "${SUBJECTS}" \
    --db "${SUBJECTS}" \
    --selfid \
    --id 0.5 \
    --blast6out results.blast6

head -n 3 results.blast6

chout1 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Suessiales|Suessiaceae|Gymnodinium|simplex|simplex_X   chvA18  50.7    148 70  1   1   162 1   148 *   *
chout2 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Suessiales|Suessiaceae|Polarella|glacialis|glacialis_X chout21 63.9    241 58  8   1   234 1   235 *   *
chout20 Eukaryota|Alveolata|Dinophyta|Dinophyceae|Peridiniales|Peridiniales_X|Heterocapsa|sp|sp_X   chout21 89.8    235 17  1   1   229 1   235 *   *

The problem is the same when using -userout filename --userfields query+id+target. The problem disappears when using --notrunclabels.

When sorting, the --minsize and --maxsize options don't work (no change in the output)

Minor bug in Makefile

After a make, if I want to clean the directory, the make clean command does not remove the object files:

$ make
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o vsearch.o vsearch.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o query.o query.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o db.o db.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o dbindex.o dbindex.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o vsearch_global.o vsearch_global.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o align.o align.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o showalign.o showalign.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o util.o util.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o unique.o unique.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o userfields.o userfields.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o results.o results.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o sortbysize.o sortbysize.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o derep.o derep.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o maps.o maps.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o arch.o arch.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o sortbylength.o sortbylength.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o shuffle.o shuffle.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o align_simd.o align_simd.cc
align_simd.cc: In function ‘void search16(s16info_s*, unsigned int, unsigned int*, CELL*, short unsigned int*, short unsigned int*, short unsigned int*, short unsigned int*, char**)’:
align_simd.cc:730:11: warning: variable ‘QR_query_left’ set but not used [-Wunused-but-set-variable]
   __m128i QR_query_left, R_query_left;
           ^
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o mask.o mask.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o minheap.o minheap.cc
minheap.cc: In function ‘void minheap_test()’:
minheap.cc:210:26: warning: narrowing conversion of ‘random()’ from ‘long int’ to ‘unsigned int’ inside { } is ill-formed in C++11 [-Wnarrowing]
       elem_t x = {random(),0,1};
                          ^
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o bitmap.o bitmap.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o cluster.o cluster.cc
g++ -O3 -msse4.1 -mtune=core2 -Icityhash -Wall    -c -o searchcore.o searchcore.cc
g++  -o vsearch vsearch.o query.o db.o dbindex.o vsearch_global.o align.o showalign.o util.o unique.o userfields.o results.o sortbysize.o derep.o maps.o arch.o sortbylength.o shuffle.o align_simd.o mask.o minheap.o bitmap.o cluster.o searchcore.o cityhash/city.o -lpthread
$ make clean
rm -f *~ *.{o,S,gcda,gcno,gcov} vsearch gmon.out cityhash/*.{o,S,gcda,gcno,gcov}
$ ls *.o
align.o       arch.o    cluster.o  db.o     maps.o  minheap.o  results.o     showalign.o  sortbylength.o  unique.o      util.o            vsearch.o
align_simd.o  bitmap.o  dbindex.o  derep.o  mask.o  query.o    searchcore.o  shuffle.o    sortbysize.o    userfields.o  vsearch_global.o

I tried different modifications of the Makefile, but nothing worked. If I launch the rm command manually, the objects files are deleted as expected. I am using GNU Make 4.0.

That's weird.

Parallelisation with pthreads

Prioritized features/options

Here are some often used options that should be implemented (from Frédéric):

Chimera checking

--usersort
--chimeras FILENAME
--abskew 2.0
--uchime FILENAME
--log FILENAME
--nonchimeras FILENAME

Clustering

--cluster FILENAME
--usersort
--id 0.97
--slots 16769023
--sizeout
--uc FILENAME
--w 64
--seedsout FILENAME
--maxrejects 500
--sizein
--log FILENAME

Add --quiet option

Implement --cons_truncate clustering output option

Bug in sortbysize

Sortbysize does not work. Will be fixed in next version.

Clustering similar to usearch

Segfault in dereplication

Check if illegal options are specified for the different commands

Add --log option to uchime_ref (and other commands?)

Reference Chimera checking.

Hi,
I'm testing out vsearch where I was using USEARCH. I have used vsearch uchime_denovo with these files and it worked great. I am running into a few problems with the uchime_ref command. My uchime command looked like this:

usearch -uchime_ref derep.fasta  -db custom_nt_database.fasta -chimeras chimeras.fasta -nonchimeras non_chimeras.fasta -uchimeout derep.uchime -log uchime_ref.log --abskew 2

Which would run fine. I removed the log command as it is not supported in vsearch and ran:

vsearch --uchime_ref derep.fasta  --db custom_nt_database.fasta --chimeras chimeras.fasta --nonchimeras non_chimeras.fasta --uchimeout derep.uchime --abskew 2

However with vsearch I get an error that I don't understand when I run that command.

Reading file custom_nt_database.fasta 100%  
684514 nt in 1036 seqs, min 204, max 5800, avg 661
WARNING: 47 sequences longer than 50000 nucleotides discarded.
Indexing sequences 100%  
Counting unique k-mers 100%  
Creating index of unique k-mers 100%  
Detecting chimeras 0%*** glibc detected *** /home/labop/vsearch/vsearch-1.0.5-linux-x86_64: realloc(): invalid next size: 0x00007fc300025400 ***
*** glibc detected *** /home/labop/vsearch/vsearch-1.0.5-linux-x86_64: double free or corruption (out): 0x00007fc31002bae0 ***

And it just hangs there.

I've used the vsearch derep_fulllength to get derep.fasta.

head of my derep.fasta gives:

>Lib_pool22_M00366:44:000000000-A92F0:1:1117:6872:7515_1:N:0:4;size=4033;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAGAAGGGTGCCTATTATACACCAATAACAGCACTAGATTTCGAA
GGTCTATATCCATCAATTATGATGGCACATAATTTATGTTATTCATCAATGGTTATGGATTCTAAATACGAAAATATACC
TGGTATAACATACGAAACGTTTGGGTTTTATAAGTTTGCACAAGATGTTCCAAGTCTTTTACCAAGTATTCTTCTAGAAC
TAAAACAGTTTCGTAAACAAGCTAAAAAGGATATGGCACAATCATCCGGTGCACTAAA
>Lib_pool21_M00366:44:000000000-A92F0:1:1115:18924:7329_1:N:0:3;size=2837;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAAAAGGGGGCATATTATACACCGATTACAGCCCTAGATTTTGAA
GCACTATATCCATCAATCATGATGGCACACAATCTATGTTACTCTACACTCGTTATGGATGAGTATCGGTACGGTAATGT
TCCTGGTATAAACTATGAATCATTCAAAATTGGCGACAAAATATATAAATTTGCGCAAGGTGTACCCAGTCTTTTACCCG
labop@mesa:/Data/Julia/Thesis/Thesis-Overall/Miseq_Initial_Run_Processing/scripts$ head ../results/Divided_by_primers/All_libs_combined_OTUs_picked/Total_AVS_R1_annotated_trimmed_paired_derep.fasta 
>Lib_pool22_M00366:44:000000000-A92F0:1:1117:6872:7515_1:N:0:4;size=4033;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAGAAGGGTGCCTATTATACACCAATAACAGCACTAGATTTCGAA
GGTCTATATCCATCAATTATGATGGCACATAATTTATGTTATTCATCAATGGTTATGGATTCTAAATACGAAAATATACC
TGGTATAACATACGAAACGTTTGGGTTTTATAAGTTTGCACAAGATGTTCCAAGTCTTTTACCAAGTATTCTTCTAGAAC
TAAAACAGTTTCGTAAACAAGCTAAAAAGGATATGGCACAATCATCCGGTGCACTAAA
>Lib_pool21_M00366:44:000000000-A92F0:1:1115:18924:7329_1:N:0:3;size=2837;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAAAAGGGGGCATATTATACACCGATTACAGCCCTAGATTTTGAA
GCACTATATCCATCAATCATGATGGCACACAATCTATGTTACTCTACACTCGTTATGGATGAGTATCGGTACGGTAATGT
TCCTGGTATAAACTATGAATCATTCAAAATTGGCGACAAAATATATAAATTTGCGCAAGGTGTACCCAGTCTTTTACCCG
CAATTCTTCTTGAACTCAAACAATTTCGTAAAAAAGCTAAAAAGGATATGGCCGCTGCGAC

So I wondered if --uchime_ref is having a problem with the line wrapping generated from derep?

So I removed the wrapping from the derep.fasta so that it looked like this head:

>Lib_pool22_M00366:44:000000000-A92F0:1:1117:6872:7515_1:N:0:4;size=4033;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAGAAGGGTGCCTATTATACACCAATAACAGCACTAGATTTCGAAGGTCTATATCCATCAATTATGATGGCACATAATTTATGTTATTCATCAATGGTTATGGATTCTAAATACGAAAATATACCTGGTATAACATACGAAACGTTTGGGTTTTATAAGTTTGCACAAGATGTTCCAAGTCTTTTACCAAGTATTCTTCTAGAACTAAAACAGTTTCGTAAACAAGCTAAAAAGGATATGGCACAATCATCCGGTGCACTAAA
>Lib_pool21_M00366:44:000000000-A92F0:1:1115:18924:7329_1:N:0:3;size=2837;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGACGCACAAAAGGGGGCATATTATACACCGATTACAGCCCTAGATTTTGAAGCACTATATCCATCAATCATGATGGCACACAATCTATGTTACTCTACACTCGTTATGGATGAGTATCGGTACGGTAATGTTCCTGGTATAAACTATGAATCATTCAAAATTGGCGACAAAATATATAAATTTGCGCAAGGTGTACCCAGTCTTTTACCCGCAATTCTTCTTGAACTCAAACAATTTCGTAAAAAAGCTAAAAAGGATATGGCCGCTGCGAC
>Lib_pool57_M00366:44:000000000-A92F0:1:1101:6581:8826_1:N:0:39;size=2079;
AATTCGCCCTTGAGGGGGCGACGGTGCTGGATGCACAAAAGGGGGCATATTATACACCGATTACAGCCCTAGATTTTGAAGCACTATATCCATCAATCATGATGGCACACAATCTATGTTACTCTACACTCGTTATGGATGAGTATCGGTACGGTAATGTTCCTGGTATAAACTATGAATCATTCAAAATTGGCGACAAAATATATAAATTTGCGCAAGGTGTACCCAGTCTTTTACCCGCAATTCTTCTTGAACTCAAACAATTTCGTAAAAAAGCTAAAAAGGATATGGCCGCTGCGAC
>Lib_pool2_Julia.2_2886829_M01170:9:000000000-A6JT1:1:1108:16338:8526_1:N:0:0;size=1828;
GAGGGGGCGACGGTGTTGGACGCCATGCGTAAACAGTTATGGGTCAGCCAAGCTGCACCTAACCCCAAGCTGCGTATGTCCAACATTGGAAAGCCTTGTAGCAGGGCGCTCTGGTACGACATCAATGGTGACGAGCAAGCGGAGAGCTTTACGCCACAGACCAAGCTCAAGTTCATCGTAGGTGACATTGTTGAAGCAATGTTGATCTACTTGGCCAAGGAAGCTGGACATTCTGTGACTGAGATGCAATC
>Lib_pool1_Julia.1_8087229_M01170:9:000000000-A6JT1:1:2107:10781:13498_1:N:0:0;size=1809;

But when I run the vsearch uchime command with it I get a new error:

Reading file custom_nt_database.fasta 100%  
684514 nt in 1036 seqs, min 204, max 5800, avg 661
WARNING: 47 sequences longer than 50000 nucleotides discarded.
Indexing sequences 100%  
Counting unique k-mers 100%  
Creating index of unique k-mers 100%  
Detecting chimeras 0%

Fatal error: Illegal header line in query fasta file

Any ideas of things to try?

Write --help to stdout

Please consider making --help write to stdout so that it can be piped to a pager, grep, etc. Thanks!

Support for reading from / writing to pipes (input from stdin and output to stdout)

fix for maxaccepts=0 and maxrejects=0

Support for long (>15 nt) and gapped seeds

Requires changes to code that identified unique seeds in a sequence. Also the current direct index of kmers to a list of the database sequences where the kmer is present must be changed to a hash table.

Shuffling not random on unix with --seed 0

Ok on Mac

Add support for amino acid sequences

The README says that this "may be added in the future." Any updates on those plans? Would adding AA support be a large undertaking or something as simple as adding a new alphabet and scoring matrix?

Vectorization of global alignment - single query vs multiple targets

Difference between --blast6out and --userout for id field

in --blast6out, the third field (id) should be 0.0 (real) when there is no alignment, not 0 (integer). The behavior is correct in --userout.

Speed-up searching when using the --top_hits_only option

During a search, the target sequences are sorted by decreasing k-mer composition similarity with the query sequence. The target sequences that are the most likely to get a high similarity score with the query are the first to be treated (i.e. pairwise aligned).

In a taxonomic assignment context (with --maxaccepts 0 --maxrejects 0 --id 0.5), and when using the option --top_hits_only, vsearch ends up comparing all the targets to the query, which is inefficient. The search process could stop once the query-target similarities start to drop. In practice, the algorithm could count the number of blocks (a block contains 16 targets if I am correct) where observed similarities are all lower than the best obtained so far. After two blocks of decreasing similarity values, the search can safely stop and start working on the next query.

Search both strands

Stable sorting to avoid input-order dependencies

Sorting performed by cluster_fast and uchime_denovo should be stable to avoid input order dependencies. Sort by length and labels (and sequences?) as secondary keys.

Clustering

USEARCH clustering commands:

cluster_fast
cluster_smallmem
cluster_otus

The cluster_otus is a more complex clustering method that involves chimera detection.

Clustering options suggested:

--cluster FILENAME
--usersort
--id 0.97
--slots 16769023
--sizeout
--uc FILENAME
--w 64
--seedsout FILENAME
--maxrejects 500
--sizein
--log FILENAME

Support for --sizein and --sizeout when clustering

Read input sizes when --sizein is specified and write size of clusters to centroids file if --sizeout is specified.

Manuscript

A manuscript describing VSEARCH should be submitted to a preprint journal.

Bug in masking

With the latest vsearch version (0.2.0), when masking the vderep.fsa file, the output is empty when using --qmask none or --qmask soft. The output is ok when using --qmask dust:

${VSEARCH} --maskfasta vderep.fsa --output vderep.fsa_masked -qmask soft
wc -l < vderep.fsa_masked

The computation takes a few seconds with dust, but is immediate with none and soft.

Difference between --blast6out and --userout for evalue and bits fields

With the last version (0.2.2), the evalue and bits fields return different values. With --blast6out the fields are set to "*". With --userout evalue = -1 and bits = 0.

The documentation says it should be evalue = -1 and bits = 0.

Add option to define identity as including terminal gaps

Vectorization of global alignment - single query vs single target

SIMD-vectorization of global alignment of one query vs one target sequence. Should use the striped approach of Farrar (2007) or the vertical vector approach of Rognes (2000).

Add support for bzipped/gzipped fasta files

Chimera detection

Commands for chimera detection:

uchime_ref
uparse_ref

Suggested options for chimera checking

--usersort
--chimeras FILENAME
--abskew 2.0
--uchime FILENAME
--log FILENAME
--nonchimeras FILENAME

add usearch's fastq_mergepairs function and options

I had been using a function in usearch to detect overlap of paired-end reads. This used the following options which would be nice to have available in vsearch.

usearch_v.7.0.1090_i86linux32 -fastq_mergepairs read1.fq 
                              -reverse read2.fq
                              -fastq_maxdiffs 6
                              -fastq_truncqual 2 
                              -fastq_minlen 36 
                              -fastq_minmergelen 50 
                              -fastqout merge.fq 
                              -fastq_allowmergestagger

Illegal instruction fault on some cpus

vsearch 1.0.1 crashes with an illegal instruction fault on some cpus, notably AMD cpus. The problem is identified (it is due to lack of ssse3 and/or sse41 features) and a fix (that only requires sse2) is in the works.