Using CAZy Database with SAMSA2

[damien...@brown.edu]

Hello,

I am hoping to use the SAMSA2 pipeline to analyze my metatranscriptomic data using the CAZy database but I am having some trouble modifying the analysis counter scripts to parse the resulting m8 files to count and condense the output. Whenever I try running the command, I get the following error message: 

Starting database analysis now.

Traceback (most recent call last):

  File "/users/djcabral/scratch/samsa2/python_scripts/DIAMOND_analysis_counter.py", line 138, in <module>

    db_entry = db_entry[1][:-1]

IndexError: list index out of range

I downloaded the CAZy database in fasta format from the dbCAN website (http://csbl.bmb.uga.edu/dbCAN/download.php) and an example of the format can be seen below:

>ACI19890.1|GT2

MVNEKVAIVILNYNGWQDTIECLESVYQIDYPNYDVILVDNASSEDSIVKIKEYCDGKIKVESKFFEYNPNNKPIMVAEYEKNELEDCSIDPDFRNLPSNRKLILIKNDENYGFAGGNNVGVRFAMKYLKPEFILFLNNDTVVDEKFLMSLVGKIKDNSIIGSAQSLLLRKDCETIDSLGQELKMRGLQISDKAMGEKISNYKDIKDAEIFGPCCAAAIYRSEVLNKIGFFDEDFFYIYEDVDLSFRNRLAGFKSYLIVDSMVFHKRGISGGKDKTISMKSYFSLRNYLSILIRYYPSNLIIRNLPKILYIISRLFVSSLYFKKLNETIMFFYKSLRYRYSIEDKLKLRMIIKEWVR

>ASB58792.1|GT4

MTKKILFCATVDYHFKAFHLPYFKWFKQMGWEVHVAANGQTKLPYVDEKFSIPIRRSPFDPQNLAVYRQLKKVIDTYEYDIVHCHTPVGGVLARLAARQARRHGTKVLYTAHGFHFCKGAPMKNWLLYYPVEKWLSAYTDCLITINEEDYIRAKGLQRPGGRTQKIHGIGVNTERFRPVSPIEQQRLREKHGFREDDFILVYPAELNLNKNQKQLIEAAALLKEKIPSLRLVFAGEGAMEHTYQTLAEKLGASANVCFYGFCSDIHELIQLADVSVASSIREGLGMNVLEGMAAEKPAIATDNRGHREIIRDGENGFLVKIGDSAAFARRIEQLYHKPELCRKLGQEGRKTALRFSEARTVEEMADIYSAYMDMDTKEKSV

>ACZ83404.1|GT0

MSEQTPEEPKYQRLGPINWLPEERKVIERFEERVRDLERRLAIAEARADYAQWKLESTKVQRPYRVAEAITGAKGSGVVRLPGKLRMALRPRKGPEAPRPVAEVIEELDQRGPAVDVPQVKWPAGPINRPDLKVAVILDDFSRMAFKYEWDQIEFGLRDWPEIFAERRPELLFVESAWHGNQGRWRYQMTGSNAPKPELRALIDWCREEGIPTVFWNKEDPPNFDFFIDTAKLFDYVFTCDGDMVPRYREALGHDRIDVLQFAAQPRVHNPIQERRGRLHDVVFAGMYFRDKHPERREQMETVLAPIRELGLHIFARNGTVDEKYAWPEEYVPHIVGELPYDQMLAAYKMYKVFLNVNSVLDSPTMCARRVFELSACSTSVVSGWSRAVQETFGPLIPIAHDELESYNQVLHLINSPELRARQGHLAMREVFDKHLFTHRVDQILGFLGRPVEQRSRSVSVVLPTNRASQLEHAIASVAKQIHRPLQLVMVLHGLDIDPVVVADKARMAGISDVVVLPADPSLSLGACMNLGIAAAEGDLIAKMDDDNLYGEHYLSDLVRAFDYSDAELVGKGAHYAYFEGRNTTMLRMPGLEHRYSWLVQGGTFLGKADMFRHYGFEDVTRGEDTRLVRRLKEDRVKIYSADRFNFVYWRSGDPSMHTWQPDDHKLTRGAQFSFVGHPDAHVMI

Do you have any suggestions of how I can modify the DIAMOND_analysis_counter.py scripts to accommodate this database? Thanks in advance for all your help!

Best,

Damien

[Sam Westreich]

Hi Damien,

Yes, I know what's going wrong here - the DIAMOND_analysis_counter.py script looks for a space in order to separate the ID from function/organism annotation:

>XXX123.1  Function  [Organism]

The script is having an error here because there's no space in the CAZy database header lines.

The simplest modification here might be to just replace the pipe symbol (|) in the CAZy database with a couple of spaces.  This can be done with a one-line command:

sed 's/|/  /g' cazy_db.fasta > cazy_db_w_spaces.fasta

Since the CAZy database doesn't contain organism annotations, you'll want to run the DIAMOND_analysis_counter.py script with the -F flag, to get functional results.

I tested this on my machine and it works, so let me know if you have issues, or if you want to modify the DIAMOND_analysis_counter.py script to split the database on a pipe instead of on a space.

Best,

Sam

--
You received this message because you are subscribed to the Google Groups "SAMSA bioinformatics group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to samsa-bioinformatics-group+unsub...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/samsa-bioinformatics-group/147fa48d-12dd-40b7-a061-7baf3e017a56%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Sam Westreich

Microbiome Scientist, DNAnexus, 

http://www.mosaicbiome.com

[damien...@brown.edu]

Hi Sam,

Thanks for the quick response! I tried your suggestion and that seems to have helped, but now I am running into another error:

Now reading through the m8 results infile.

Analysis of /users/djcabral/scratch/Metatranscriptomics/SAMSA2/CAZy/A1T12.new_cazy.m8 complete.

Number of total lines: 104528

Number of unique sequences: 104528

Time elapsed: 0.21 seconds.

Starting database analysis now.

Success!

Time elapsed: 9.05 seconds.

Number of lines: 863368

Number of errors: 0

Dictionary database assembled.

Time elapsed: 9.05 seconds.

Number of errors: 0

Top ten function matches:

Traceback (most recent call last):

  File "/users/djcabral/scratch/samsa2/python_scripts/DIAMOND_analysis_counter.py", line 229, in <module>

    for k, v in sorted(condensed_RefSeq_hit_db.items(), key=lambda k,v: -v)[:10]:

TypeError: <lambda>() takes exactly 2 arguments (1 given)

Any suggestions of how to fix that?

-Damien

[Sam Westreich]

Hi Damien,

Hmm.  You can try pulling the latest version of the script from Github (https://github.com/transcript/samsa2/blob/master/python_scripts/DIAMOND_analysis_counter.py), but I'm not sure that's the issue here.  The updated script has brackets around the k, v given to lambda to better explicitly declare them.

If this doesn't work, you could send me the first couple thousand lines of the input query file, and I can see if I can replicate the error.

Best,

Sam

--

You received this message because you are subscribed to the Google Groups "SAMSA bioinformatics group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to samsa-bioinformatics-group+unsub...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/samsa-bioinformatics-group/88308b92-5af2-40ec-9795-48fa60c9e774%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.