AGAT
AGAT copied to clipboard
NBISweden
agat_sp_load_function_from_protein_align.pl - works only for Exonerate via MAKER.
I used pblat to align proteins to my draft genome, then used psl2gff to generate a gff3 file.
When I run agat_sp_load_function_from_protein_align.pl like so: agat_sp_load_function_from_protein_align.pl --annotation manage_functional_annotation/strict_filter_A_atri_tm2_pnm_braker_plus_EDTA.gff3 -pgff pblat/Combined_insect_proteins2A_atri_genome_m2.gff -pfasta /home/soppenheim/mendel-nas1/EvolAppl/Proteins/Align/Combined_insect_proteins.fasta --method add -w --pe 1 -v --value 80 --output load_function_from_protein_align
the script fails with numerous instances of: Use of uninitialized value in numeric ne (!=) at /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/AGAT/OmniscientTool.pm line 2600
And when I run it with the species flag, like so: agat_sp_load_function_from_protein_align.pl --annotation manage_functional_annotation/strict_filter_A_atri_tm2_pnm_braker_plus_EDTA.gff3 -pgff pblat/Combined_insect_proteins2A_atri_genome_m2.gff -pfasta /home/soppenheim/mendel-nas1/EvolAppl/Proteins/Align/Combined_insect_proteins.fasta --method add -w --pe 1 -v --value 80 --sp "Nasonia vitripennis, Trichogramma brassicae, Trichogramma sarcophagae, Trichogramma pretiosum" --output sp_load_function_from_protein_align
it fails with: Can't locate object method "new" via package "XML::Twig" at /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/Bio/DB/Taxonomy/entrez.pm line 676.
I did a fresh install of agat via conda in hopes of fixing, but the results did not change. Other scripts (most recently agat_sp_manage_functional_annotation.pl) are running fine.
I am grateful for any help! Sara
Additional errors: Problem ! protein ID not found into the protein fasta file! Use of uninitialized value $prot_tag in concatenation (.) or string at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 1000. We cannot check the real protein length, let's continue without this one: I cannot calculate the absolute match because the tag Gap is absent ! I cannot calculate the absolute match because the tag Gap is absent ! Use of uninitialized value $correct_id in pattern match (m//) at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 814. Use of uninitialized value $correct_id in pattern match (m//) at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 818. Use of uninitialized value $correct_id in pattern match (m//) at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 822. Use of uninitialized value $seq_id_correct in lc at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 791. Use of uninitialized value $seq_id_correct in concatenation (.) or string at /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl line 803.
The first issue make me think that there are error in the strand column of your gff file. Could you show a sample of your file?
The second error is a dependency missing. I will have to add it in conda. For now you can install it using conda install -c bioconda perl-xml-twig
For the third problem it is a bit more complex. Aligned protein must have a Name attribute in the gff, and a CIGAR string like in a Gap attribute (https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md#the-gap-attribute).
Exonerate produce such type of output, it sounds that it is not the case using pblat + psl2gff. I should probably update the Readme to be more precise about expected input.
Hi Jaques, I'm going to try an exonerate run and I've installed the missing dependency.
In the meantime, here are some lines from my pblat-based gff:
Aatri_ctg9 BLAT match_part 1933664 1933714 5.75 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid1;Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 55 105
Aatri_scaf2 BLAT match_part 30326210 30326231 3.64 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid2;Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 192 213
Aatri_scaf3 BLAT match_part 25427846 25427867 7.85 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid3; Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 39 60
Aatri_scaf4 BLAT match_part 28345344 28345359 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1023 1038
Aatri_scaf4 BLAT match_part 28346954 28346993 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1082 1121
Aatri_scaf4 BLAT match_part 28347104 28347129 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1132 1157
The native psl output is too messy to paste in a useful way, but here's a link to a description of the format: https://useast.ensembl.org/info/website/upload/psl.html
Thank you for your fast response! Sara
On Fri, Aug 6, 2021 at 4:07 AM Jacques Dainat @.***> wrote:
The first issue make me think that there are error in the strand column of your gff file. Could you show a sample of your file?
The second error is a dependency missing. I will have to add it in conda. For now you can install it using conda install -c bioconda perl-xml-twig
For the third problem it is a bit more complex. Aligned protein must have a Name attribute in the gff, and a CIGAR string like in a Gap attribute ( https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md#the-gap-attribute ). Exonerate produce such type of output, it sounds that it is not the case using pblat + psl2gff. I should probably update the Readme to be more precise about expected input.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-894084225, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJOAUBFY2LEGCZDJAYDT3OJ5PANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
Follow up: when I tried installing perl-xml-twig in my agat env, I got "All requested packages already installed."
On Fri, Aug 6, 2021 at 1:01 PM Sara Oppenheim @.***> wrote:
Hi Jaques, I'm going to try an exonerate run and I've installed the missing dependency.
In the meantime, here are some lines from my pblat-based gff: Aatri_ctg9 BLAT match_part 1933664 1933714 5.75 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid1;Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 55 105 Aatri_scaf2 BLAT match_part 30326210 30326231 3.64 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid2;Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 192 213 Aatri_scaf3 BLAT match_part 25427846 25427867 7.85 ++ . Parent=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI_mid3;Target=tr|A0A7M7M7L2|A0A7M7M7L2_NASVI 39 60 Aatri_scaf4 BLAT match_part 28345344 28345359 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1023 1038 Aatri_scaf4 BLAT match_part 28346954 28346993 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1082 1121 Aatri_scaf4 BLAT match_part 28347104 28347129 7.03 +- . Parent=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI_mid8;Target=tr|A0A7M7LNR4|A0A7M7LNR4_NASVI 1132 1157
The native psl output is too messy to paste in a useful way, but here's a link to a description of the format: https://useast.ensembl.org/info/website/upload/psl.html
Thank you for your fast response! Sara
On Fri, Aug 6, 2021 at 4:07 AM Jacques Dainat @.***> wrote:
The first issue make me think that there are error in the strand column of your gff file. Could you show a sample of your file?
The second error is a dependency missing. I will have to add it in conda. For now you can install it using conda install -c bioconda perl-xml-twig
For the third problem it is a bit more complex. Aligned protein must have a Name attribute in the gff, and a CIGAR string like in a Gap attribute ( https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md#the-gap-attribute ). Exonerate produce such type of output, it sounds that it is not the case using pblat + psl2gff. I should probably update the Readme to be more precise about expected input.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-894084225, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJOAUBFY2LEGCZDJAYDT3OJ5PANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
I now have some exonerate output but I'm still getting errors when I run the script: -------------------------------- parse features -------------------------------- => GFF version parser used: 2 [Fri Aug 6 17:34:01 2021] Parsing: 0%
------------- EXCEPTION: Bio::Root::Exception ------------- MSG: [vulgar: tr|Q17ER3|Q17ER3_AEDAE 0 124 . Aatri_scaf1 35696922 35696547 - 573 M 6 18 G 0 3 M 118 354] does not look like GFF2 to me STACK: Error::throw STACK: Bio::Root::Root::throw /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/Bio/Root/Root.pm:447 STACK: Bio::Tools::GFF::_from_gff2_string /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/Bio/Tools/GFF.pm:523 STACK: Bio::Tools::GFF::from_gff_string /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/Bio/Tools/GFF.pm:439 STACK: Bio::Tools::GFF::next_feature /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/Bio/Tools/GFF.pm:395 STACK: AGAT::OmniscientI::slurp_gff3_file_JD /mendel-nas1/soppenheim/miniconda3/envs/agat/lib/site_perl/5.26.2/AGAT/OmniscientI.pm:317 STACK: /home/soppenheim/mendel-nas1/miniconda3/envs/agat/bin/agat_sp_load_function_from_protein_align.pl:221
Here is the top of my exonerate output. Do I need to convert it somehow before using with agat? vulgar: tr|Q17ER3|Q17ER3_AEDAE 0 124 . Aatri_scaf1 35696922 35696547 - 573 M 6 18 G 0 3 M 118 354
--- START OF GFF DUMP ---
##gff-version 2 ##source-version exonerate:protein2genome:local 2.4.0 ##date 2021-08-06 ##type DNA
I got it to run by manually editing the exonerate output, but the output from agat_sp_load_function_from_protein_align.pl was incorrect. Agat reported "no protein aligned" for all genes, but in my exonerate and annotation gffs I can see that there are alignments. One example: From the annotation gff:
Aatri_scaf1 AUGUSTUS gene 35696545 35696922 . - . ID=AatriG00000012879;Name=his2a_26;gene_id=jg6900;locus=RLOC_00003748;makerName=jg6900;transcript_id=jg6900.t1
From the exonerate gff:
Aatri_scaf1 exonerate:protein2genome:local gene 35696548 35696922 573 - . gene_id 1 ; sequence tr|Q17ER3|Q17ER3_AEDAE ; gene_orientation . ; identity 95.16 ; similarity 97.58
From the agat_sp_load_function_from_protein_align report: No protein aligned to this gene model: aatrig00000012879
My exonerate command:
exonerate --model protein2genome --cores 24 --verbose \
--softmasktarget --showalignment no --showtargetgff --score 500 --percent 70 --bestn 1 \
--query /nas4/dbobo/sardeans/Aphelinus_genomes/Proteins/Combined_insect_proteins.fasta \
--target /nas4/dbobo/sardeans/Aphelinus_genomes/softMasked/A_atri_genome.edta_softMask2.fasta \
> Combined_insect_proteins2Aatri_m2_exonerate
Where am I going wrong? And, is there a recommended tool for converting the exonerate output into a format that agat likes?
Thanks, Sara
When developing this script I was using exonerate output From MAKER. Either you do not use the correct command, or MAKER was doing extra manipulation before outputing the exonerate gff.
The gff must contain a Name and/or product attribute(s) that will be lifted to the annotation and a Gap attribute (like Gap=M8 D3 M6 I1 M6 ) use to compute the similarity on the fly.
Could you share a few lines of a Maker-produced exonerate file so I can compare?
Sara
On Mon, Aug 9, 2021 at 6:51 AM Jacques Dainat @.***> wrote:
When developing this script I was using exonerate output From MAKER. Either you do not use the correct command, or MAKER was doing extra manipulation before outputing the exonerate gff.
The gff must contain a Name and/or product attribute(s) that will be lifted to the annotation and a Gap attribute (like Gap=M8 D3 M6 I1 M6 ) use to compute the similarity on the fly.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-895125824, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJKSYRJJR2PNA6RGSSDT36XKBANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
Yes sure:
unitig_6|quiver protein2genome protein_match 162828 163298 359 - . ID=unitig_6|quiver:hit:874:3.10.0.1;Name=B4G8R5.1
unitig_6|quiver protein2genome match_part 162828 163298 359 - . ID=unitig_6|quiver:hsp:1613:3.10.0.1;Parent=unitig_6|quiver:hit:874:3.10.0.1;Target=B4G8R5.1 130 287;Length=311;Gap=M35 I2 M51 I2 M31 D3 M37
unitig_6|quiver protein2genome protein_match 162828 163298 353 - . ID=unitig_6|quiver:hit:875:3.10.0.1;Name=B3MU92.1
unitig_6|quiver protein2genome match_part 162828 163298 353 - . ID=unitig_6|quiver:hsp:1614:3.10.0.1;Parent=unitig_6|quiver:hit:875:3.10.0.1;Target=B3MU92.1 129 286;Length=310;Gap=M35 I2 M52 I1 M6 I1 M24 D3 M37
So looking through Maker it is clear that a lot of post-processing is applied to the exonerate output, see /maker-3.01.03/lib/Widget/exonerate/ protein2genome.pm. Do you have any interest in modifying agat_sp_load_function_from_protein_align so that it can handle native exonerate output? If not I will try and hack a solution, but it is limiting if the agat script only works with maker-generated exonerate files. I'm using braker2 rather than maker and would prefer not to re-do everything in order to use agat_sp_load_function_from_protein_align.
Sara
On Mon, Aug 9, 2021 at 12:34 PM Jacques Dainat @.***> wrote:
Yes sure:
unitig_6|quiver protein2genome protein_match 162828 163298 359 - . ID=unitig_6|quiver:hit:874:3.10.0.1;Name=B4G8R5.1 unitig_6|quiver protein2genome match_part 162828 163298 359 - . ID=unitig_6|quiver:hsp:1613:3.10.0.1;Parent=unitig_6|quiver:hit:874:3.10.0.1;Target=B4G8R5.1 130 287;Length=311;Gap=M35 I2 M51 I2 M31 D3 M37 unitig_6|quiver protein2genome protein_match 162828 163298 353 - . ID=unitig_6|quiver:hit:875:3.10.0.1;Name=B3MU92.1 unitig_6|quiver protein2genome match_part 162828 163298 353 - . ID=unitig_6|quiver:hsp:1614:3.10.0.1;Parent=unitig_6|quiver:hit:875:3.10.0.1;Target=B3MU92.1 129 286;Length=310;Gap=M35 I2 M52 I1 M6 I1 M24 D3 M37
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-895368631, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJIFHP2532Y3U3GSCXDT377QDANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
Instead of using the Gap attribute to compute the identity on the fly I can directly take the identity or the similarity score. In the same way we can use the sequence instead of the Target attribute to find the reference protein file.
It should not be a big deal to implement something that. Could you send me data /file examples?
Here is my exonerate command: exonerate --model protein2genome --softmasktarget --showalignment no --showtargetgff yes --score 500 --percent 70 --bestn 10 --query testSet_insectProts.fasta --target scaffold2.fasta > insectProts2scaffold2.gff
The files are attached! insectProts2scaffold2.gff https://drive.google.com/file/d/14ajU0inRl2vfbn406Ei2qnpjgO0isahY/view?usp=drive_web scaffold2.fasta https://drive.google.com/file/d/1Ixb4xL1V0XiLqfav82v2QE_tcfh2s7KN/view?usp=drive_web testSet_insectProts.fasta https://drive.google.com/file/d/1M_nr_Vg1VjqbbpZkM71IplBDJ3SUsYgw/view?usp=drive_web
On Mon, Aug 9, 2021 at 4:22 PM Jacques Dainat @.***> wrote:
Instead of using the Gap attribute to compute the identity on the fly I can directly take the identity or the similarity score. In the same way we can use the sequence instead of the Target attribute to find the reference protein file. It should not be a big deal to implement something that. Could you send me data /file examples?
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-895515708, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJIB2ODLUYNPPWNNBPLT4A2ILANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
I have updated the example files and shared them here: https://drive.google.com/drive/folders/1QlgX2Zp06Ugk6HCoiddYP82xbvdAMz1Y?usp=sharing
My exonerate command:
exonerate --model protein2genome --showvulgar yes --showquerygff no --softmasktarget
--showalignment no --showtargetgff --score 300 --percent 50
--ryo ">%qi length=%ql alnlen=%qal\n>%ti length=%tl alnlen=%tal\n"
--target seq.fasta --query protein.fasta > exonerate.out
I used "process_exonerate_gff3.pl" (https://github.com/hyphaltip/genome-scripts/blob/master/data_format/process_exonerate_gff3.pl) to convert the exonerate output to gff3.
I need a sample of your annotatoin gff file on which you want to lift the protein information on
I added a sample to the Drive folder!
On Fri, Aug 13, 2021 at 6:43 AM Jacques Dainat @.***> wrote:
I need a sample of your annotatoin gff file on which you want to lift the protein information on
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-898367205, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJMDAACXOQV7ZWTNIFDT4TZOLANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
Actually the use of CIGAR is preferable than vulgar. We need --showcigar yes rather than showvulgar.
The GAP attribute in GFF3 is based on CIGAR. That might simplify a lot the problem
OK, I am re-running with --showcigar yes will you need a new example output file?
On Tue, Aug 17, 2021 at 4:53 AM Jacques Dainat @.***> wrote:
Actually the use of CIGAR is preferable than vulgar. We need --showcigar yes rather than showvulgar. The GAP attribute in GFF3 is based on CIGAR. That might simplify a lot the problem
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-900116511, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJJI2RTUE6LO5MFXAMDT5IPQZANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
Okay I think everything you need is in the shared folder!
On Wed, Aug 18, 2021 at 3:56 AM Jacques Dainat @.***> wrote:
yes please
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/NBISweden/AGAT/issues/158#issuecomment-900902137, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC5QDJPMZWRFAWFBDPVSMHTT5NRTDANCNFSM5BU7BKEQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .
I have not seen any easy way to fix that. I'm not able to reproduce the MAKER polishing allowing to have a neat information to judge the overlapping quality to infer or not the function. For now the script is usable only for mapping done using Exonerate via MAKER.