cazy_extract.pl downloading same sequence multiple times

[cmorganl]

Hi Dallas,

I've noticed that for the larger CAZy families (say >900 sequences) cazy_extract fails to download all of the unique fasta records from NCBI, but rather downloads multiple copies of the each fasta record until approximately the correct number of sequences have been downloaded. Here's an example: GH31 has 7637 representative sequences on CAZy. 7247 sequences are in the resulting fasta file from cazy_extract.pl, but only 982 of them are unique.

Since my perl sucks I don't know at which exact step (scraping the cazy pages or enrez e-utilities) its going wrong so any help would be appreciated!

Thanks, Connor

[DallasThomas]

Hello Connor,

So I have started looking into GH31. The cazy extract script does pull the the accession number given for each protein name/organism in each family and group. These accession numbers are used to pull the fasta information from NCBI.

It does download all of the fasta information given by Cazy. The one thing to note is not all of the fasta information is unique for each protein name/Organism. For example, lets take a look at the following four:

SAMN05428937_2635 Achromobacter sp. MFA1 R4 SIT23620.1 CVS48_02335 Achromobacter spanius DSM 23806 AUA54971.1   CLM73_03885 Achromobacter spanius MYb73 AVJ26317.1 AXYL_00816 Achromobacter xylosoxidans A8 ADP14166.1

These are 4 of the 7637 unique listed in GH31 All. When you follow each of their Accession numbers in GenBank - you have 4 unique GenBank entries - each with identical fasta sequences.

Hence the output from Cazy Extract from these four looks like this:

SIT23620.1 Alpha-glucosidase, glycosyl hydrolase family GH31 [Achromobacter sp. MFA1 R4] MPPKFDFAHTTQLDSVELLTSRPSRIDFDAGDGLHLVVEAHAPGVFRLRCGEANHFNDDKPSARARAVAE MLLARQEAVGEAALAPREGGDGWRITQGEVALEILTNPVRVALYRHEDLVFESEVSAHAPAFGHNALDQA EDAVWTACFTLAAHERVFGLGETPGDLNRREESVVSDDPEHRALPLAWSPRGWGLYVNTMRRVEHSVGVA PADSAYVLTVDDAVLDLFLFAGEPTEILNQYSALTGRAGQPVLWAMGAWLQQAAGEMPTQTVALVARMRE NQIPVDAVTLAQPAAWGFQADKTVFEWDAQRFPDAKQMLALFHKHNVHVAASGFPGVIKTSPLFEELEDR GWLLTRDDGSAQVFDGNAATSGQPYGLLDLTYRDIYNLWVERHRQLVEDGLDAPACDAQIAIPDGVTARN GETGPALRTMYPLLVRRALFDAVAGLKVPPEGVVPSADLFPAAQRLPWQVGPQVDNTWEGLRHSLRTSLS IGASGVPVQVHNIGSAARDRSGMTSELYLRWLTAGIFSANFAFQGVDGLMPWDFGDDVLAHAKTWMQWRY RLVPYVLGAIEDSARTGLPVQRSMALSFPNDPHAHEWDLQYLLGPALLVAPITEPGNQVRVYLPKGEAWW DLNTGHRYEGGTTWTIDCELGQFPVFGREGHMLCLGPAAQHTGEFNSARILDEVWMFGMPVHNPVVMRNK IRVMQMQGSSYIKGLEGLRISPSEGLEVKRRGAEVRISRAR AUA54971.1 glycosyl hydrolase [Achromobacter spanius] MPPKFDFAHTTQLDTIELLTSRPSRFDFDAGDGLHLVVEAHAPGVFRLRCGEANHLNDDKPGARARAVAE MLLARQEAVGEASIAPREGGDGWRITQGDVALEIQSNPVRVAIYRNEERVFESEVSAHTPAFGHNALDQE DDAVWTAGFTLLADERIFGLGETPGDLNRREESVVSDDPEHRALPLAWSPRGWGVYVNTMRRVEHSVGVA PADSAYVLTVDDAVLDVFLFAGEPAEILNQYTALTGRAGQPVLWAMGAWLQQATGEMPTQTVALLARMRE NQIPVDAVTLAQPAAWGFQSDKAVFEWDAQRFPDAKQMLALFHKHNVHVAASGFPGVLKTSPLFEELEDR GWLLAKEDGAAQVFDGIAATSGQPFGLLDLTYRDIYNMWVERHRQLIEDGLDAPACDAQLAIPDGVTARN GETGPALRSMYPLLVRRALFDAVAGLKVPPEGVVPSSDLFPAAQRLPWQVGPQAENSWDGLRHTLRTALS IGASGVPVQVHGIGSASRDRAGMTSELYLRWLTVGVFSANFSFQGVEGLMPWDFGDETLAQAKTWMQWRY RLVPYVLGAIEDSARTGLPVQRSMALSFPNDPHAHEWDLQYLLGPALLVAPVTEPGKQVRVYLPKGEAWW DLNTGHRYEGGTTWTIDCELGQFPVFGREGHMLCLGPAAQHTGEFNSARILDEVWMFGMPVHNPVVMRNK IRVMQMQGSSYIKGLEGLRISPSEGLEVKRRGAEVRISRAR AVJ26317.1 glycosyl hydrolase [Achromobacter spanius] MPPKFDFAHTTQLDNVELLTSRPSRFDFDAGDGLHLVVEAHAPGVFRLRCGEANHFNDDKPGARARAVAE MLLARQEAVGEGAIAPREGGDGWRITQGDVALEIQTSPVRVALYRNEDRVFESEVSAHTPAFGYNALDQA DDAVWTAGFTLAADERVFGLGETPGDLNRREESVVSDDPEHRALPLAWSPRGWGVYANSMRRVEHSVGVA PADSAYVLTVDDAVLDLFLFAGEPAEILNQYTALTGRAGQPVLWAMGAWLQQAAGEMPTQTVALVARMRE NQIPVDAVTLAQPIAWGFQADKTVFEWDAQRFPDAKQMLALFHKHNVHVAASGFPGVIKTSPLFEELEDR GWLLTKDDGAAQVFDGNAATSGQPYGLLDLTYRDIYNLWVERHRQLVEDGLDAPACDAQIAIPDGVTARN GETGAALRTMYPLLVRRALFDAVAGLKVPPEGVVPSADLFPAAQRLPWQVGPQVDNTWEGLRHTLRTSLS IGSSGVPVQVHNIGSAARDRSGMTSELYLRWLTAGVFSANFAFQGVEGLMPWDFGDETLSHAKTWMQWRY RLVPYVLGAIEDSARTGLPVQRSMALSFPNDPHAHEWDLQYLLGPALLVAPITEPGNQVRVYLPKGEAWW DLNTGHRYEGGTTWTIDCEPGQFPVFGREGHMLCLGPAAQHTGEFNSARILDEVWMFGMPVHNPVVMRNK IRVMQMQGSSYIKGLEGLRISPSEGLEVKRRGAEVRISRAR ADP14166.1 glycosyl hydrolase, family 31 protein [Achromobacter xylosoxidans A8] MPPKFDFAHTTQLDNVELLTSRPSRIDFDAGDGLHLVVEAHAPGVFRLRCGEANHLNDDKPGARARAVAE MLLARQEAVGEAVIAPRDDGEGWRITQGDVSLEILSNPVRVAIYRNDERVFESEVSAHTPAFGHNALDQA DDAVWTAGFTLQSDERVYGLGETPGDLNRREESVVSDDPEHRALPLAWSPRGWGVYVNTMRRVEHSVGVA PADSAYVLTVDDVVLDVFLFAGEPGEILNQYTALTGRAGQPSLWAMGAWLQQAPGEMPTDTVALVARMRE NQIPVDAVTLAQPAAWGFQSDKTVFEWDAQRFPDAKQMLALFHKHHVHVAASGFPGVLKTSPLFEELEDR GWLLAKDDGSAQVFDGNAATSGQPYGLLDLTYRDIYNLWVERHRQLVEDGLDAPACDAQIAIPDGVTARN GETGPALRSMYPLLVRRALFDAVAGLKVPPEGVVPSSDLFPAAQRLPWQVGPQAENTWEGLRHTLRTALS IGASGLPVQVHGIGSAARDRASMTPELYLRWLTAGVFSSNFAFQGVDGLMPWDFGEETLAHAKTWMQWRY RLVPYVLGAIEDSARTGLPVQRSMAMSFPNDPVSHDWDLQYLLGPALLVAPITEPGKQVRVYLPKGEAWW DLNTGHRYEGGTTWTIDCELGQFPVFGREGHMLCLGPAAQHTGEFNSARILDEVWMFGMPVHNPVVMRNK IRVMQMQGSSYIKGLEGLRISPSEGLEVKRRGAEVRISRAR

This is very common and expected as this is how ancestral relatedness and variance is discovered.

[DallasThomas]

I also have it default that fragments are removed, which will also cause a decrease in number of output sequences. You can disable this with the Fragments flag when calling either Saccharis or cazy_extract

[cmorganl]

Sorry for the confusion, I was referring to the fasta headers and not just the sequences. I agree, unique sequences would be expected. I don't think this problem arises from nature, though. Here are the commands (and results in comments) I've used to test this:

perl cazy_extract.pl --family GH31 --group all --Fragment false
grep -c "^>" GH31_cazy.fasta # Number of headers was 7992
grep "^>" GH31_cazy.fasta | sort | uniq | wc -l # Number of unique headers was 999
perl cazy_extract.pl --family GH31 --group all # Just to make sure this wasn't due to the fragment flag
grep -c "^>" GH31_cazy.fasta # 7896
grep "^>" GH31_cazy.fasta | sort | uniq | wc -l # 982

Can you reproduce this?

[DallasThomas]

Was able to reproduce and similar issue creeping up in a new way, arrrgghhhh.

Updated - I do sincerely apologize. Please give this one a shot.

[cmorganl]

No worries, Dallas! I really appreciate the software - definitely beats the alternative of endlessly clicking!

I think it works now! Feel free to close the issue.

Thanks again!