using vsearch with qiime?

[Vanessa V.]

Hi, 

I was wondering if someone could let me know how I can use the vsearch de novo clustering with QIIME?  I do not see this as an option in the pick_otus.py?

Thanks for your help in advance, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

Qiime 1 does not officially use vsearch, but this is still possible one of two ways. 

1. use vsearch directly to pick de novo OTUs outside of qiime. Then return to qiime for taxonomy assignment, tree building etc.

2. use vsearch to replace the default OTU picking program (which is uclust), and use the qiime script directly. 

Interested in trying one of these?

Colin

[Vanessa V.]

Hi Colin, 

Yes! So would I pass -m vesearch to use your second option?  For the pick_otus.py step?

[Colin Brislawn]

Hello Vanessa,

Vsearch is an open-source, 64-bit implementation of usearch, so option 2 is a sort of hack that uses vsearch to replace '-m usearch61'. 

Here is how it works. You tell the script pick_otus.py to use usearch by passing '-m usearch61.' But then we use an alias to have 'usearch61' actually run the program 'vsearch' instead. You can set up this alias by running this command:

alias usearch61='/absolute/path/to/vsearch'

To test that it's worked, try running

usearch61

which should then output a vsearch help file, showing that commands sent to 'usearch61' will actually reach vsearch. 

Colin

[Vanessa V.]

Ah, OK.  I will give this a shot and let you know how it goes!  

Thanks!!!

[Vanessa V.]

Hi Colin, 

I tried the vsearch with pick_open_reference_otus.py on our HPCC as SBATCH script.  It ran but how do I verify it actually ran with vsearch?  In the log file it shows in the pick_otus.py step with -m as usearch61_ref.  Making me think it did not run it as vsearch?

Thanks for your help in advance!

Vanessa

#SBATCH --Nodes=1

#SBATCH --ntasks=40

alias usearch61='~/../lab/tools/vsearch/1.11.1/bin/vsearch'

pick_open_reference_otus.py -i /gpfs0/home/resgoodman/vav002/170106_amplicons_qiime/seqs_chimeras_pos_filtered.fna -r /reference/greengenes/gg_13_8_otus/rep_set/97_otus.fasta -o /gpfs0/home/resgoodman/vav002/170106_amplicons_qiime/pos_OTUs_gg -p -m usearch61 -f -a -O 40

[Colin Brislawn]

Hello Vanessa,

Great question! To track this one down, look for the log files left over in the OTU picking folder. Depending on the settings, some of the vsearch / usearch log files will list the settings at the top of the file (before the results).

You could also add this line into you sbatch script, just before your pick_open_reference_otus script

usearch61 -v

This should actually reach the vsearch program and return something like

vsearch v1.11.1_linux_x86_64, 16.0GB RAM, 8 cores

https://github.com/torognes/vsearch

Keep in touch,

Colin

[Vanessa V.]

Hi Colin, 

It gives me this error in the log file.  Says its invalid.

Wed Feb 15 18:00:12 EST 2017

+ module load python-2.7.10

++ /usr/bin/modulecmd bash load python-2.7.10

+ eval LOADEDMODULES=python-2.7.10 ';export' 'LOADEDMODULES;MANPATH=/gpfs0/export/opt/anaconda-2.3.0/share:/usr/man:/usr/share/man:/usr/local/man:/usr/local/share/man:/usr/X11R6/man:/$

++ LOADEDMODULES=python-2.7.10

++ export LOADEDMODULES

++ MANPATH=/gpfs0/export/opt/anaconda-2.3.0/share:/usr/man:/usr/share/man:/usr/local/man:/usr/local/share/man:/usr/X11R6/man:/act/man

++ export MANPATH

++ PATH=/gpfs0/export/opt/anaconda-2.3.0/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/sbin:/act/bin:/act/bin:/gpfs0/home/resgoodman/vav002/bin

++ export PATH

++ PYTHONPATH=/gpfs0/export/opt/anaconda-2.3.0/lib/python2.7/site-packages/

++ export PYTHONPATH

++ _LMFILES_=/act/modulefiles/python-2.7.10

++ export _LMFILES_

+ module load qiime/1.9.1

++ /usr/bin/modulecmd bash load qiime/1.9.1

+ eval LOADEDMODULES=python-2.7.10:qiime/1.9.1 ';export' 'LOADEDMODULES;PATH=/opt/qiime/1.9.1/../usearch:/opt/qiime/1.9.1/../seqprep/bin:/opt/qiime/1.9.1/../ea-utils/bin:/opt/qiime/1.$

++ LOADEDMODULES=python-2.7.10:qiime/1.9.1

++ export LOADEDMODULES

++ PATH=/opt/qiime/1.9.1/../usearch:/opt/qiime/1.9.1/../seqprep/bin:/opt/qiime/1.9.1/../ea-utils/bin:/opt/qiime/1.9.1/bin:/gpfs0/export/opt/anaconda-2.3.0/bin:/usr/lib64/qt-3.3/bin:/u$

++ export PATH

++ PYTHONPATH=/opt/qiime/1.9.1/lib/python2.7/site-packages:/gpfs0/export/opt/anaconda-2.3.0/lib/python2.7/site-packages/

++ export PYTHONPATH

++ QIIME_CONFIG_FP=/opt/qiime/1.9.1/.qiime_config

++ export QIIME_CONFIG_FP

++ _LMFILES_=/act/modulefiles/python-2.7.10:/act/modulefiles/qiime/1.9.1

++ export _LMFILES_

+ cd /home/resgoodman/vav002/170601_amplicons_qiime

/var/spool/slurmd/job340501/slurm_script: line 34: cd: /home/resgoodman/vav002/170601_amplicons_qiime: No such file or directory

+ alias 'usearch61=~/../lab/tools/vsearch/1.11.1/bin/vsearch'

+ usearch61 -v

Invalid command line

Option --v is invalid

+ pick_open_reference_otus.py -i /gpfs0/home/resgoodman/vav002/170106_amplicons_qiime/seqs_chimeras_pos_filtered.fna -r /reference/greengenes/gg_13_8_otus/rep_set/97_otus.fasta -o /gp$

[Colin Brislawn]

Hello Vanessa,

Maybe this old version of vsearch does not have a version flag, or maybe it wants --v or --version.

Also, it looks like this directory had an issue too, so I would double check the paths throughout the document.

+ cd /home/resgoodman/vav002/170601_amplicons_qiime

/var/spool/slurmd/job340501/slurm_script: line 34: cd: /home/resgoodman/vav002/170601_amplicons_qiime: No such file or directory

Let me know what you find!

Colin

[Vanessa V.]

Hi Colin, 

I'm still having trouble with this.  It seems to work fine on the command line on the log in node.  I can do the alias and then check with usearch61 -v and it gives the output you listed on Feb. 15.  But once I use this with SBATCH via SLURM on our cluster, it gives me the invalid error once I use usearch61 -v (or the others you had me try).  It shows me in the log file that it was an invalid command.  I definitely don't want to run the pick otus step via the log in node.  Is there another way to check that it actually ran vsearch instead of usearch61? 

Thanks for your help, 

Vanessa

I attached the log file after running the script - seemingly with the vsearch alias.  Any clues it ran vsearch that I'm missing??

[Colin Brislawn]

Hello Vanessa,

Great detective work. Something about submitting the bash script (vs running it on the SLURM node) is causing this to fail.

What steps do you perform on the command line before you run qiime or vsearch? Do you activate an environment, or load a module, or something like that? How do you make sure your SBATCH script does these same things so that it can also use qiime?

We may have to do something special to make sure the SBATCH script knows to use the other version of vsearch.

Colin

[Vanessa V.]

Hi Colin, 

I have attached the .sh file that I use.  The first time I ran it, I did not have the usearch61 -v line in there.  It did not give me any errors, so I *think* alias is working.  Even with the invalid usearch -v in there, it says invalid option but it still runs the rest of the code.  I just don't want to assume its running vsearch when it's really not for some reason, I just wanted some way to double check!

I am also attaching the slurm log too (pick_otus).  Where it says invalid option.  

Thanks, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

Thanks for sending me both of those. I think we are really close, and I've got a bunch of small ideas which will help us close the gap.

I noticed that the alias uses your home directory. To be safe, you could the absolute path.

The slurm script lists 

module load python-2.7.10
module load qiime/1.9.1

Is there a module load vsearch? This might help vsearch get to the node.

If the program does not like -v, you could try --v or -version or --version or just usearch61 by itself. I'm hoping to find some setting which will give an output, and that output will let us make sure it's working. 

We are deep in the 'guess and check' phase of troubleshooting, so let me know any clues you find along the way.

Thanks Vanessa,

Colin

[Vanessa V.]

Hi Colin, 

So I actually did get the usearch61 --version to work in the SBATCH script to check if alias worked.  It did not run alias as vsearch but came back as usearch61 despite the alias command in the line above it (see below).  I will test some of the ideas you listed above to get vsearch to work (just as a note, the alias command and running usearch61 -v does work on the command line in terms of coming back as vsearch, it's just something in the sbatch/slurm code is not quite right).  

I am not sure what you mean by using the absolute path not the home directory for vsearch.  What I have listed below is exactly where vsearch is located on our cluster?  Can you explain that in more detail?

So should I add the below to where I have module load python, qiime, etc. in the sbatch?

module use ~/../lab/tools/share/modulefiles module load vsearch 

From slurm log file: 

alias usearch61='~/../lab/tools/vsearch/1.11.1/bin/vsearch'

usearch61 --version

usearch_i86linux32 v6.1.544

Thanks a ton for your help, I really appreciate it!

Vanessa

[Colin Brislawn]

Hello Vanessa,

Thanks for the detailed updates. I'm glad we were able to get that version and confirm that, yes, it's still a little funky on the SBATCH queue. 

I think now would be a good time to contact the folks who support the supercomputer. They should be able to close this final gap when it comes to making one of the worker nodes recognize this alias. This should be easy for them. 

As for 'absolute path' I meant replace this 

alias usearch61='~/../lab/tools/vsearch/1.11.1/bin/vsearch'

with something that avoids using the ~ like this

alias usearch61='/root/folder/folder/lab/tools/vsearch/1.11.1/bin/vsearch'

The team that build and updated that vsearch module should be able to solve this pretty quick, so I would just ask them.

Colin

[Vanessa V.]

Hi Colin, 

Indeed with the help from HPC staff, I think we got the alias to work in SBATCH.  She said: 

Alias is not expanded in a non-interactive shell script, in this case, the batch script. To make it work, you need to change the shell option. 

add a line above the alias definition: "shopt -s expand_aliases “  That should get you the right command vsearch. :)

It seems to work by doing this, I get the vsearch version to come up after I run usearch61 --version!

alias 'usearch61=~/../lab/tools/vsearch/1.11.1/bin/vsearch'

+ /gpfs0/home/resgoodman/vav002/../lab/tools/vsearch/1.11.1/bin/vsearch --version

vsearch v1.11.1_linux_x86_64, 125.9GB RAM, 20 cores

https://github.com/torognes/vsearch

However, I am still not 100% sure it's running vsearch during pick_open_reference_otus.py!!  I put in -m as usearch61 when I submit the SBATCH for this script, but if you see below (from the log file) it has -m as usearch61_ref during pick_otus.py.  Is it not going to use the alias because it's not usearch61, it's technically running usearch61_ref?? 

# Pick Reference OTUs command

pick_otus.py -i /gpfs0/home/resgoodman/vav002/170106_amplicons_qiime/seqs_chimeras_pos_filtered.fna -o /gpfs0/home/resgoodman/vav002/170106_amplicons_qiime/pos_OTUs_gg_vsearch/step1_otus -r /reference/greengenes/gg_13_8_otus/rep_set/97_otus.fasta -m usearch61_ref  --suppress_new_clusters

Thanks for your help!

Vanessa

[Colin Brislawn]

Hello Vanessa,

I find HPC support folks so supportive. I'm glad you have a friendly team over there.

However, I am still not 100% sure it's running vsearch during pick_open_reference_otus.py!!  I put in -m as usearch61 when I submit the SBATCH for this script, but if you see below (from the log file) it has -m as usearch61_ref during pick_otus.py.  Is it not going to use the alias because it's not usearch61, it's technically running usearch61_ref?? 

Good question. Remember, we are 'tricking' the qiime script to use vsearch. Everything that qiime does and says will be usearch related. Because qiime thinks it's running usearch, so we can't look for confirmation here. 

Instead we need to look for confirmation that our alias trick worked, by calling a 'usearch61' command and finding a vsearch result underneath. 

I get the vsearch version to come up after I run usearch61 --version!

Exactly right!

Based on the code you copied in, it looks like you ran  

/gpfs0/home/resgoodman/vav002/../lab/tools/vsearch/1.11.1/bin/vsearch --version

but this may just be how the output works. Our goal is to run

usearch61 --version

and get that a vsearch result. Running usearch61 is exactly what the qiime script will do, so this is how we make sure our alisas worked.

We are so close!

Colin

PS. Happy Friday.

[Vanessa V.]

Hi Colin, 

I did input usearch61 --version (as a line in my sbatch script) but in the log file, it writes it like below.  

alias 'usearch61=~/../lab/tools/vsearch/1.11.1/bin/vsearch'

+ /gpfs0/home/resgoodman/vav002/../lab/tools/vsearch/1.11.1/bin/vsearch --version

vsearch v1.11.1_linux_x86_64, 125.9GB RAM, 20 cores

https://github.com/torognes/vsearch

However, I did find this though when looking at pick_open_reference otus, 

http://qiime.org/scripts/pick_open_reference_otus.html

-m, --otu_picking_method

The OTU picking method to use for reference and de novo steps. Passing usearch61, for example, means that usearch61 will be used for the de novo steps and usearch61_ref will be used for reference steps. [default: uclust]

I think it might be running a vsearch and usearch61_ref combo??!!  Is there another alias that we can use?  Can you explain why vsearch is the best way to go? :)

Thanks for all of your help, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

Ah, OK. Looks like the alias is working just like we want, and any time qiime calls usearch61, it's getting out version of vsearch.

I think it might be running a vsearch and usearch61_ref combo??!!  Is there another alias that we can use?  Can you explain why vsearch is the best way to go? :)

Remember that qiime is going to say usearch, because we are sneakily replacing the usearch program with the vsearch program. The help guides and errors from qiime will all keep saying usearch. We will know that we have swapped in vsearch, and ignore them.

Part of the confusion here is that qiime uses the program usearch (which we are replacing) and several workflows for otu picking, called things like -m usearch61 and -m usearch61_ref. These all use the program usearch61 (which we have replaced!) to do OTU picking in different ways. 

The program usearch, vsearch, uclust all perform search and clustering in a similar way. They combine the speed of a k-mer filter with the precision of an alignment. Out of these options, I like vsearch because it's

• always free
• open source
• 64-bit (for large databases and clusters)
• uses exact, optimal alignments

When you requested, "I can use the vsearch de novo clustering with QIIME?"

Why do you want to use vsearch de novo :-)

Colin

[Vanessa V.]

Hi Colin, 

OK, that makes sense.  Although I will say that I get the same answer in terms of OTUs generated from usearch61 (when alias wasn't working previously) to now using vsearch with alias working.  Is that to be expected?  That's why I was thinking something still wasn't right here...

Yes, I kept seeing here on the forum and elsewhere that vsearch was the better to use over uclust for OTU picking so I wanted to try it out and compare.  :)

Thanks, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

Thanks for sticking with me while we figured this out. Under a month for a strange technical issue is pretty good. 

Because they are based on the same algorithm, we would expect similar results. Probably not identical through! While the number of OTUs may be the same, are all the exact same OTU centroids selected? Is the number of reads in each OTU in each sample the same (like, are the OTU tables identical)? That would make me think vsearch was not being used at all...

Let's keep comparing them and see!

Colin

[Vanessa V.]

Hi Colin, 

I am looking at a mock community positive control sample.  It is only one sample.  I was using this to vet the pipeline to use for my actual samples (I wanted to run through the pipelines/clustering algorithms quickly and see which gave me the best representation of the mock community).  I didn't look at centroid sequences within the OTUs but I did look at the taxa summaries and the OTUs and # sequences within OTUs.  The usearch and vsearch results are exactly identical on every front.  I also ran uclust and while the results are similar there are definitely differences, the results are not identical.  I have attached the biom table results in case you want to take a peek.  This is why I'm not sure that vsearch is actually being run - we wouldn't expect the results from usearch and vsearch to be exactly identical, right?  

Thanks, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

we wouldn't expect the results from usearch and vsearch to be exactly identical, right?  

Right!

Well... unless it was in a very simple sample, with a single 'right' answer. I mean, if you clustered 5 reads, different programs could give you the same answer. If you clustered 20x copies of the same 5 reads, you would still get identical results, just like if you clustered 50x or 1000x of those reads. Simple communities yield simple, consistent results. Still, we are trying to benchmark programs with real data so I would expect some changes.

I've cc'ed one of the qiime devs who may be able to help us more. Let's see where this takes us... 

Any advice on using vsearch through one of the qiime workflow scripts?

Colin

[TonyWalters]

Hello,

As long as vsearch is being seen in the $PATH environment as the "usearch61" executable, it should be using it anytime you call pick_otus.py or identify_chimeric_seqs.py with -m usearch61 (or usearch61_ref). You might convince yourself that it's calling that executable by renaming the vsearch "usearch61" executable file that's in the $PATH, and confirming that the code errors out when calling pick_otus with an error about not finding the software.

But for direct comparisons of OTU picking methods, I would go outside of QIIME itself, and call uclust, usearch (there are versions of this not wrapped in QIIME), and vsearch directly, although that will involve manually calling several of the steps that are wrapped in QIIME. See code for denovo OTU picking with usearch61/vsearch here: https://github.com/biocore/burrito-fillings/blob/master/bfillings/usearch.py#L1860 Direct comparisons would probably have to be at the individual steps involved, but I'd imagine most of the differences would be at the step of actual clustering (i.e. heuristics to minimize comparisons and percent identity calculations).

[Colin Brislawn]

Good morning Tony. Thanks for 'qiimeing' in. 

We have been playing around with path variables too, and I'm not feeling confident that vsearch is being called. Perhaps Vanessa's HPC support could meet with her to test this too.  

I really like your idea about going through OTU picking using the programs directly. This would let us identify the exact steps in which changes occur, and track these changes through to the final OTU table.  

I was using this to vet the pipeline to use for my actual samples (I wanted to run through the pipelines/clustering algorithms quickly and see which gave me the best representation of the mock community). 

Vanessa, I think running them directly is a great way to do that. Also, if you try any of the new non-OTU methods like Swarm or dada2, this will let you explore their conceptual and practical differences.

Let me know what you try!

Colin

[Vanessa V.]

Hi all, 

Thank you very much for your help.  I was hoping to be able to run these in QIIME directly as I'm trying to establish a pipeline for our lab group within QIIME.  I'm testing a few other programs as well.  If it's not easy to run Vsearch within the QIIME environment, I think it might not be the best choice to use as of now.  I will definitely look into some of the other clustering methods within QIIME that are available.  How do these cluster the data (Swarm or dada2) without OTUs?? Interesting!

I'll see if our HPC team has any ideas, too.  Do you think vsearch will eventually be implemented as an otu picking method in any of the newer versions of QIIME?

Thanks again, 

Vanessa

[Colin Brislawn]

Hello Vanessa,

Do you think vsearch will eventually be implemented as an otu picking method in any of the newer versions of QIIME?

I hope so! Qiime 2 will use a plugin based design so that developers can more easily integrate their software. Qiime 2 is still in development, so stay tuned.

How do these cluster the data (Swarm or dada2) without OTUs?? Interesting!

Instead of clustering, they perform denoising / error correction. So maybe you could say they produce Denoised Sequence Variants (DSVs)

See the computational biologists talk about what to call these post-OTUs, in this long thread.

https://github.com/benjjneb/dada2/issues/62 

Colin

[Jay T]

Colin hey we spoke earlier about Vsearch in a PM. I am able to successfully use the alias however it seems to disappear once I start a new terminal session. Is there a more permanent solution? Also using the actual alias in my workflow seems to produce an error. Would you mind looking over this? 

Thanks,
JT

MacQIIME Todds-Mac-Pro:Soil $ pick_open_reference_otus.py -m usearch61 -i soil_seqs_final/combined_seqs.fna -o vsearch_picked_otus -f -a -O 6

Traceback (most recent call last):

  File "/macqiime/anaconda/bin/pick_open_reference_otus.py", line 453, in <module>

    main()

  File "/macqiime/anaconda/bin/pick_open_reference_otus.py", line 432, in main

    minimum_failure_threshold=minimum_failure_threshold)

  File "/macqiime/anaconda/lib/python2.7/site-packages/qiime/workflow/pick_open_reference_otus.py", line 713, in pick_subsampled_open_reference_otus

    close_logger_on_success=False)

  File "/macqiime/anaconda/lib/python2.7/site-packages/qiime/workflow/util.py", line 122, in call_commands_serially

    raise WorkflowError(msg)

qiime.workflow.util.WorkflowError: 

*** ERROR RAISED DURING STEP: Pick Reference OTUs

Command run was:

 pick_otus.py -i soil_seqs_final/combined_seqs.fna -o vsearch_picked_otus/step1_otus -r /macqiime/anaconda/lib/python2.7/site-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta -m usearch61_ref  --suppress_new_clusters

Command returned exit status: 1

Stdout:

Stderr

Traceback (most recent call last):

  File "/macqiime/anaconda/bin/pick_otus.py", line 1004, in <module>

    main()

  File "/macqiime/anaconda/bin/pick_otus.py", line 897, in main

    otu_prefix=otu_prefix, HALT_EXEC=False)

  File "/macqiime/anaconda/lib/python2.7/site-packages/qiime/pick_otus.py", line 1800, in __call__

    HALT_EXEC=HALT_EXEC

  File "/macqiime/anaconda/lib/python2.7/site-packages/bfillings/usearch.py", line 1844, in usearch61_ref_cluster

    raise ApplicationError('Error running usearch61. Possible causes are '

burrito.util.ApplicationError: Error running usearch61. Possible causes are unsupported version (current supported version is usearch v6.1.544) is installed or improperly formatted input file was provided

[Colin Brislawn]

Hello JT,

Glad you got it working!

however it seems to disappear once I start a new terminal session

I'm not super familiar with alias, but this works for me; adding an alias to my .bashrc file. 

https://askubuntu.com/questions/17536/how-do-i-create-a-permanent-bash-alias 

When new terminals are opened, that alias you add will be automatically configured for you. 

Depending on how you run your workflow scripts, getting it to recognize an alias could be tricky. What software platform is used to run your workflow?

Colin

[Jay T]

Hi Colin

I am using a Mac so I'm within the MacQiime environment. 

Thanks, 
JT

[Colin Brislawn]

Hello JT,

Ah, OK. MacQiime should play nice with any of the alisas you may want to pass it.

It's possible that this is a compatability issue between usearch61 and vsearch. This issue was reported here, so hopefully it will be address soon. 

https://github.com/torognes/vsearch/issues/229 

If the incompatibility has to do with the version of vsearch and not the alias, we will just have to wait for a fix. :-/

Colin

[Jay T]

Do I need to have usearch61 installed? I'm working on that and it seems to be a pain in the you know what.

[Colin Brislawn]

Nope. Usearch is an optional dependency, which you do not need to perform analysis. 

Qiime includes the program uclust (and older version of usearch), which works pretty well.

Colin

[Jay T]

Colin - I know this is contrary to what you said but it seems after following the instructions below from the MacQiime web site --Vsearch now works. I had to install both usearch5 and usearch6.1. Then create an alias for usearch61.

USEARCH - If you don't have USEARCH installed, you can get it here. There are two required versions and you will need both in QIIME 1.9.0: USEARCH v5.2.236 and USEARCH 6.1, each required for DIFFERENT scripts, unfortunately. Name the 5.2.236 executable "usearch" and the 6.1 executable "usearch61" and make sure they're in your PATH.

I installed both and it finally ran. 

¯\_(ツ)_/¯

Thanks!

JT

[Colin Brislawn]

Oh yeah, some scripts do make use of usearch and you will need it to run those scripts. For other scripts, usearch is not needed. 

If it works, great!

Colin