Mercurial > repos > bgruening > tapscan
changeset 0:196795831b6a draft
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/tapscan commit 2e1cd301fb38af8a1e9267fc60fcb5ca3c576aeb
| author | bgruening |
|---|---|
| date | Wed, 14 Feb 2024 13:54:16 +0000 |
| parents | |
| children | c4f865bd101a |
| files | tapscan.xml tapscan_coverage_values_v10.txt tapscan_domains_v12.txt.gz tapscan_rules_v81.txt tapscan_script_v74.pl test-data/PUBLIC_Ectocarpus-sp7_proteins_head.fa test-data/output.1.tsv test-data/output.2.tsv test-data/output.3.tsv test-data/output.domtbl.tsv |
| diffstat | 10 files changed, 1713 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tapscan.xml Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,126 @@ +<tool id="tapscan_classify" name="TAPScan Classify" version="4.74+galaxy0" profile="23.0"> + <description>Detect Transcription Associated Proteins (TAPs)</description> + <edam_topics> + <edam_topic>topic_0121</edam_topic> + </edam_topics> + <requirements> + <requirement type="package" version="3.3.2">hmmer</requirement> + <requirement type="package" version="5.26">perl</requirement> + <requirement type="package" version="4.8">sed</requirement> + </requirements> + <required_files> + <include type="literal" path="tapscan_script_v74.pl"/> + <include type="literal" path="tapscan_domains_v12.txt"/> + <include type="literal" path="tapscan_rules_v81.txt"/> + <include type="literal" path="tapscan_coverage_values_v10.txt"/> + </required_files> + <command detect_errors="aggressive"><![CDATA[ + +hmmsearch + --domtblout domtblout.txt + --cut_ga + '${__tool_directory__}/tapscan_domains_v12.txt.gz' + '$protein_fasta_in' + +&& + +perl '${__tool_directory__}/tapscan_script_v74.pl' + domtblout.txt + '${__tool_directory__}/tapscan_rules_v81.txt' + '$taps_detected' + '$taps_family_counts' + '$taps_detected_extra' + '${__tool_directory__}/tapscan_coverage_values_v10.txt' + +&& + +## make the outputs tab-separated for Galaxy compatibility +sed -i -e 's/;/\t/' -e 's/;/\t/' -e 's/;/\t/' -e "1d" '$taps_detected' && +sed -i -e 's/;/\t/g' -e "1d" '$taps_family_counts' && +sed -i -e 's/;/\t/4' -e 's/;/\t/' -e 's/;/\t/' -e 's/;/\t/' -e "1d" '$taps_detected_extra' && + +## add header lines for clarity +sed -i '1s/^/sequence ID\tTAP family\tnumber of classifications\tdomains\n/' '$taps_detected' && +sed -i '1s/^/TAP family\tnumber of detected proteins\n/' '$taps_family_counts' && +sed -i '1s/^/sequence ID\tTAP family\tsubfamily\tnumber of classifications\tdomains\n/' '$taps_detected_extra' + + ]]></command> + <inputs> + <param name="protein_fasta_in" type="data" format="fasta" optional="false" label="Proteins in FASTA format" help=""/> + <param name="domtblout" type="boolean" checked="false" label="Output the HMMer domain hits table?"/> + </inputs> + <outputs> + <data name="taps_detected" format="tabular" label="${tool.name} on ${on_string}: Detected TAPs - domains and assigned TAP family +for each gene ID"> + <actions> + <action name="column_names" type="metadata" default="sequence ID,TAP family,number of classifications,domains" /> + </actions> + </data> + <data name="taps_family_counts" format="tabular" label="${tool.name} on ${on_string}: Count - Summary of the number of +members for each TAP family"> + <actions> + <action name="column_names" type="metadata" default="TAP family,number of detected proteins" /> + </actions> + </data> + <data name="taps_detected_extra" format="tabular" label="${tool.name} on ${on_string}: Detected TAPs Extra - with subfamiliy information"> + <actions> + <action name="column_names" type="metadata" default="sequence ID,TAP family,subfamily,number of classifications,domains" /> + </actions> + </data> + <data name="domtbl" format="tabular" from_work_dir="domtblout.txt" label="${tool.name} on ${on_string}: HMMer Domain Hits Table"> + <filter>domtblout</filter> + </data> + </outputs> + <tests> + <test expect_num_outputs="3"> + <param name="protein_fasta_in" value="PUBLIC_Ectocarpus-sp7_proteins_head.fa" ftype="fasta"/> + <output name="taps_detected" file="output.1.tsv" ftype="tabular" lines_diff="2" /> + <output name="taps_family_counts" file="output.2.tsv" ftype="tabular" lines_diff="2" /> + <output name="taps_detected_extra" file="output.3.tsv" ftype="tabular" lines_diff="2" /> + </test> + <test expect_num_outputs="4"><!-- test with domtblout --> + <param name="protein_fasta_in" value="PUBLIC_Ectocarpus-sp7_proteins_head.fa" ftype="fasta"/> + <param name="domtblout" value="true"/> + <output name="taps_detected" file="output.1.tsv" ftype="tabular" /> + <output name="taps_family_counts" file="output.2.tsv" ftype="tabular" /> + <output name="taps_detected_extra" file="output.3.tsv" ftype="tabular" /> + <output name="domtbl" file="output.domtbl.tsv" ftype="tabular" lines_diff="10"/> + </test> + </tests> + <help><![CDATA[ +**What it does** + +TAPscan is a comprehensive tool for annotating TAPs with a special focus on species +belonging to the Archaeplastida. In general, the detection of TAPs is based on the detection +of highly conserved protein domains. + +During the first step, each sequence out of a species protein set is scanned for protein domains +(stored as profile Hidden Markov Models) using hmmsearch. The domains list consists of 154 profile +HMMs and functions as the domain reference during the hmmsearch command. + +Afterwards, by running TAPscan, specialized rules are applied to finally +assign the protein sequences to TAP families based on the detected domains in the previous step. +With the latest TAPscan v4, a protein set can be scanned for 137 different TAP families with high +accuracy through applying GA-thresholds and coverage values. + +**Output Files** + +TAPscan provides the user with three different output files. Each output file is +tab-separated. + +- **Output 1: "Detected TAPs"** - contains the detected domains and finally assigned TAP family + for each gene ID. If domains are assigned to a sequence but not all rules are fulfilled, the + sequence is assigned to “0_no_family_found”. +- **Output 2: "Family Counts"** is a summary of the number of members for each TAP family. +- **Output 3: "Detected TAPs Extra"** - is similar to output 1 but contains additional + information about subfamilies. + + ]]></help> + <creator> + <organization name="Rensing Lab" url="https://plantco.de"/> + </creator> + <citations> + <!-- TODO: add citation for TAPscan v4 paper when published --> + <citation type="doi">10.3390/genes12071055</citation> + </citations> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tapscan_coverage_values_v10.txt Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,154 @@ +Acetyltransf_1 0.0931034483 +AP2 0.3214285714 +ARID 0.0516393443 +AUX_IAA 0.1077981651 +Auxin_resp 0.4308510638 +B3 0.196875 +BES1_N 0.3888888889 +BSD 0.485915493 +BTB 0.3105095541 +bZIP_1 0.64453125 +bZIP_2 0.5948275862 +C1_2 0.1612903226 +CAF1C_H4-bd 0.1764705882 +CBFB_NFYA 0.2019230769 +CCT 0.4952830189 +CG-1 0.5568181818 +CSD 0.5878378378 +DDT 0.2678571429 +DEAD 0.1115591398 +dsrm 0.2330097087 +DUF260 0.3638059701 +DUF296 0.155075188 +DUF547 0.1181672026 +DUF573 0.3515625 +DUF632 0.3302603037 +DUF702 0.2314126394 +E2F_TDP 0.2078313253 +EIN3 0.5874125874 +FHA 0.1777456647 +FLO_LFY 0.4569138277 +FYRC 0.3928571429 +FYRN 0.46875 +GAGA_bind 0.2265774379 +GATA 0.5406976744 +GRAS 0.3990825688 +Helicase_C 0.118852459 +HLH 0.3191489362 +HMG_box 0.6866197183 +Homeobox 0.5485074627 +HSF_DNA-bind 0.071641791 +IQ 0.6428571429 +JmjC 0.5058139535 +JmjN 0.4261363636 +K-box 0.438 +KNOX1 0.578125 +KNOX2 0.6346153846 +LIM 0.4726027397 +MBF1 0.3049450549 +Med26 0.1939655172 +Med31 0.451048951 +Med6 0.1577868852 +Med7 0.1818181818 +MEKHLA 0.4357541899 +mTERF 0.4850543478 +Myb_DNA-binding 0.3636363636 +NAM 0.0142857143 +O-FucT 0.0914866582 +Ovate 0.3137755102 +PAH 0.1194267516 +PAZ 0.1564569536 +PC4 0.3848684211 +PHD 0.37 +Piwi 0.4082278481 +PLATZ 0.30625 +PP2C 0.4110962567 +QLQ 0.61875 +RB_B 0.1663987138 +Rcd1 0.4131355932 +Response_reg 0.5017241379 +RFX_DNA_binding 0.5257731959 +RHD_DNA_bind 0.4581447964 +Ribonuclease_3 0.0935013263 +RRN3 0.2342427093 +Runt 0.5714285714 +RWP-RK 0.4542253521 +S1FA 0.7386363636 +SBP 0.4274193548 +SET 0.0392857143 +SH2 0.4197247706 +Sigma70_r2 0.1675531915 +Sigma70_r3 0.6346153846 +Sigma70_r4 0.6388888889 +SIR2 0.45703125 +SNF2_N 0.1150895141 +SRF-TF 0.406779661 +SSXT 0.5856164384 +START 0.5067567568 +STAT_bind 0.3781869688 +SWIB 0.3267326733 +SWIRM 0.2532467532 +TANGO2 0.125 +TCP 0.219665272 +TCR 0.4943181818 +TEA 0.4078282828 +TF_AP-2 0.5660377358 +Tfb2 0.1757668712 +tify 0.5955882353 +Tub 0.0872576177 +VEFS-Box 0.4918831169 +WD40 0.0563909774 +WHIM1 0.5263157895 +Whirly 0.5648148148 +WRC 0.6223404255 +WRKY 0.2196969697 +WSD 0.1768867925 +YABBY 0.2804621849 +zf-AN1 0.4027777778 +zf-B_box 0.2746478873 +zf-C2H2 0.5080645161 +zf-C5HC2 0.253125 +zf-CCCH 0.5689655172 +zf-Dof 0.5955882353 +ZF-HD_dimer 0.5294117647 +zf-MIZ 0.6805555556 +zf-TAZ 0.1566455696 +zf-ZPR1 0.1608910891 +Zn_clus 0.5480769231 +Alfin-like 0.75 +BEL 0.75 +DNC 0.75 +FIE_clipped_for_HMM 0.75 +G2-like_Domain 0.75 +HRT 0.75 +KNOXC 0.75 +LUFS_Domain 0.75 +NF-YB 0.75 +NF-YC 0.75 +NOZZLE 0.75 +PINTOX 0.75 +STER_AP 0.75 +trihelix 0.75 +ULT_Domain 0.75 +VARL 0.75 +VOZ_Domain 0.75 +WOX_HD 0.75 +CXC 0.75 +bZIP_AUREO 0.75 +bZIP_CDD 0.50 +ALOG 0.75 +C2H2-IDD 0.75 +zf-MYST 0.75 +CBP 0.75 +DUF3591 0.75 +LOB2 0.75 +zz-ADA2 0.75 +NLP 0.75 +CRF 0.75 +GIY_YIG 0.75 +ZPR 0.75 +LD 0.75 +NDX 0.75 +SAWADEE 0.75 +C1HDZ 0.75 +C2HDZ 0.75
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tapscan_rules_v81.txt Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,298 @@ +ABI3/VP1;AP2;should not +ABI3/VP1;Auxin_resp;should not +ABI3/VP1;B3;should +ABI3/VP1;WRKY;should not +Alfin-like;Alfin-like;should +Alfin-like;Homeobox;should not +Alfin-like;zf-TAZ;should not +Alfin-like;PHD;should not +AP2;AP2;should +AP2;CRF;should not +ARF;Auxin_resp;should +Argonaute;Piwi;should +Argonaute;PAZ;should +ARID;ARID;should +Aux/IAA;AUX_IAA;should +Aux/IAA;Auxin_resp;should not +Aux/IAA;B3;should not +BBR/BPC;GAGA_bind;should +BES1;BES1_N;should +bHLH;HLH;should +bHLH;TCP;should not +bHLH_TCP;TCP;should +bHSH;TF_AP-2;should +BSD domain containing;BSD;should +bZIP1;bZIP_1;should +bZIP1;HLH;should not +bZIP1;Homeobox;should not +bZIP2;bZIP_2;should +bZIP2;HLH;should not +bZIP2;Homeobox;should not +bZIPAUREO;bZIP_AUREO;should +bZIPAUREO;HLH;should not +bZIPAUREO;Homeobox;should not +bZIPCDD;bZIP_CDD;should +bZIPCDD;HLH;should not +bZIPCDD;Homeobox;should not +C2C2_CO-like;CCT;should +C2C2_CO-like;GATA;should not +C2C2_CO-like;tify;should not +C2C2_CO-like;PLATZ;should not +C2C2_CO-like;zf-B_box;should +C2C2_Dof;zf-Dof;should +C2C2_Dof;GATA;should not +C2C2_GATA;GATA;should +C2C2_GATA;tify;should not +C2C2_GATA;zf-Dof;should not +C2C2_YABBY;YABBY;should +C2H2;zf-C2H2;should +C2H2;zf-MIZ;should not +C3H;AP2;should not +C3H;SRF-TF;should not +C3H;MYB-2R;should not +C3H;MYB-3R;should not +C3H;MYB-4R;should not +C3H;zf-C2H2;should not +C3H;zf-CCCH;should +CAMTA;CG-1;should +CAMTA;IQ;should +NF-YA;bZIP_1;should not +NF-YA;bZIP_2;should not +NF-YA;CBFB_NFYA;should +NF-YB;NF-YB;should +NF-YB;NF-YC;should not +NF-YC;NF-YB;should not +NF-YC;NF-YC;should +NF-YC;HMG_box;should not +Coactivator p15;PC4;should +CPP;TCR;should +CSD;CSD;should +CudA;STAT_bind;should +CudA;SH2;should +DBP;DNC;should +DBP;PP2C;should +DDT;DDT;should +DDT;Homeobox;should not +DDT;Alfin-like;should not +Dicer;Piwi;should not +Dicer;DEAD;should +Dicer;Helicase_C;should +Dicer;Ribonuclease_3;should +Dicer;dsrm;should +DUF246 domain containing/O-FucT;O-FucT;should +DUF296 domain containing;DUF296;should +DUF547 domain containing;DUF547;should +DUF632 domain containing;DUF632;should +DUF833 domain containing/TANGO2;TANGO2;should +E2F/DP;E2F_TDP;should +EIL;EIN3;should +FHA;FHA;should +GARP_ARR-B_G2;CCT;should not +GARP_ARR-B_G2;G2-like_Domain;should +GARP_ARR-B_G2;Response_reg;should +GARP_ARR-B_Myb;CCT;should not +GARP_ARR-B_Myb;Response_reg;should +GARP_ARR-B_Myb;Myb_DNA-binding;should +GARP_G2-like;G2-like_Domain;should +GARP_G2-like;Response_reg;should not +GARP_G2-like;Myb_DNA-binding;should not +GeBP;DUF573;should +GIF;SSXT;should +GNAT;Acetyltransf_1;should +GNAT;PHD;should not +GRAS;GRAS;should +GRF;QLQ;should +GRF;WRC;should +C3HDZ;Homeobox;should +C3HDZ;START;should +C3HDZ;MEKHLA;should +C4HDZ;Homeobox;should +C4HDZ;START;should +C4HDZ;MEKHLA;should not +HD_PLINC;ZF-HD_dimer;should +HD_WOX;WOX_HD;should +HD_DDT;Homeobox;should +HD_DDT;DDT;should +HD_DDT;WHIM1;should +HD_DDT;WSD;should +HD_PHD;PHD;should +HD_PHD;Homeobox;should +HD_PINTOX;Homeobox;should +HD_PINTOX;PINTOX;should +HD_BEL;Homeobox;should +HD_BEL;BEL;should +HD_KNOX1;Homeobox;should +HD_KNOX1;KNOX1;should +HD_KNOX1;KNOX2;should +HD_KNOX1;KNOXC;should not +HD_KNOX2;Homeobox;should +HD_KNOX2;KNOX1;should +HD_KNOX2;KNOX2;should +HD_KNOX2;KNOXC;should +HD-other;EIN3;should not +HD-other;Homeobox;should +HD-other;bZIP_1;should not +HD-other;WOX_HD;should not +HD-other;PINTOX;should not +HD-other;PHD;should not +HD-other;BEL;should not +HMG;ARID;should not +HMG;HMG_box;should +HMG;YABBY;should not +HRT;HRT;should +HSF;HSF_DNA-bind;should +IWS1;Med26;should +Jumonji_PKDM7;JmjC;should +Jumonji_PKDM7;JmjN;should +Jumonji_PKDM7;zf-C5HC2;should +Jumonji_PKDM7;FYRN;should +Jumonji_PKDM7;FYRC;should +Jumonji_Other;JmjC;should +LFY;FLO_LFY;should +LIM;two_or_more_LIM;should +LUG;LUFS_Domain;should +MADS;SRF-TF;should +MADS;K-box;should not +MADS_MIKC;SRF-TF;should +MADS_MIKC;K-box;should +MBF1;MBF1;should +Med6;Med6;should +Med7;Med7;should +mTERF;mTERF;should +MYB-2R;G2-like_Domain;should not +MYB-2R;Response_reg;should not +MYB-2R;trihelix;should not +MYB-2R;MYB-2R;should +MYB-3R;G2-like_Domain;should not +MYB-3R;Response_reg;should not +MYB-3R;trihelix;should not +MYB-3R;MYB-3R;should +MYB-4R;G2-like_Domain;should not +MYB-4R;Response_reg;should not +MYB-4R;trihelix;should not +MYB-4R;MYB-4R;should +MYB-related;ARID;should not +MYB-related;G2-like_Domain;should not +MYB-related;Myb_DNA-binding;should +MYB-related;Response_reg;should not +MYB-related;trihelix;should not +MYB-related;MYB-2R;should not +MYB-related;MYB-3R;should not +MYB-related;MYB-4R;should not +NAC;NAM;should +NZZ;NOZZLE;should +OFP;Ovate;should +PcG_EZ;CXC;should +PcG_EZ;SET;should +PcG_FIE;FIE_clipped_for_HMM;should +PcG_FIE;WD40;should +PcG_VEFS;VEFS-Box;should +PcG_VEFS;zf-C2H2;should not +PcG_MSI;WD40;should +PcG_MSI;CAF1C_H4-bd;should +PcG_MSI;FIE_clipped_for_HMM;should not +PHD;Myb_DNA-binding;should not +PHD;Alfin-like;should not +PHD;ARID;should not +PHD;DDT;should not +PHD;Homeobox;should not +PHD;JmjC;should not +PHD;JmjN;should not +PHD;PHD;should +PHD;SWIB;should not +PHD;zf-TAZ;should not +PHD;zf-MIZ;should not +PHD;zf-CCCH;should not +PHD;HMG_box;should not +PLATZ;PLATZ;should +Pseudo ARR-B;CCT;should +Pseudo ARR-B;Response_reg;should +Pseudo ARR-B;tify;should not +RB;RB_B;should +Rcd1-like;Rcd1;should +Rel;RHD_DNA_bind;should +RF-X;RFX_DNA_binding;should +RRN3;RRN3;should +Runt;Runt;should +S1Fa-like;S1FA;should +SAP;STER_AP;should +SBP;SBP;should +SET;zf-C2H2;should not +SET;TCR;should not +SET;CXC;should not +SET;PHD;should not +SET;Myb_DNA-binding;should not +SET;SET;should +Sigma70-like;Sigma70_r2;should +Sigma70-like;Sigma70_r3;should +Sigma70-like;Sigma70_r4;should +Sin3;PAH;should +Sin3;WRKY;should not +Sir2;SIR2;should +SOH1;Med31;should +SRS;DUF702;should +SWI/SNF_BAF60b;SWIB;should +SWI/SNF_SNF2;AP2;should not +SWI/SNF_SNF2;PHD;should not +SWI/SNF_SNF2;SNF2_N;should +SWI/SNF_SNF2;zf-CCCH;should not +SWI/SNF_SNF2;Myb_DNA-binding;should not +SWI/SNF_SNF2;HMG_box;should not +SWI/SNF_SWI3;SWIRM;should +SWI/SNF_SWI3;Myb_DNA-binding;should +TEA;TEA;should +TFb2;Tfb2;should +tify;tify;should +TRAF;BTB;should +TRAF;zf-TAZ;should not +Trihelix;trihelix;should +TUB;Tub;should +ULT;ULT_Domain;should +VARL;VARL;should +VOZ;VOZ_Domain;should +Whirly;Whirly;should +WRKY;WRKY;should +Zinc finger, AN1 and A20 type;zf-AN1;should +Zinc finger, AN1 and A20 type;zf-C2H2;should not +Zinc finger, MIZ type;zf-MIZ;should +Zinc finger, MIZ type;zf-C2H2;should not +Zinc finger, ZPR1;zf-ZPR1;should +Zn_clus;Zn_clus;should +ALOG;ALOG;should +C2H2;C2H2-IDD;should not +C2H2_IDD;C2H2-IDD;should +C2H2_IDD;zf-C2H2;should +MYST;zf-MYST;should +CBP;CBP;should +CBP;zf-TAZ;should +CBP;BTB;should not +TAFII250;DUF3591;should +LOB1;bZIP_1;should not +LOB1;bZIP_2;should not +LOB1;DUF260;should +LOB1;HLH;should not +LOB1;Homeobox;should not +LOB2;LOB2;should +LDL/FLD;SWIRM;should +LDL/FLD;Myb_DNA-binding;should not +ADA2;zz-ADA2;should +ADA2;Myb_DNA-binding;should +RKD;RWP-RK;should +RKD;NLP;should not +NLP;RWP-RK;should +NLP;NLP;should +CRF;CRF;should +CRF;AP2;should +GIY_YIG;GIY_YIG;should +ZPR;ZPR;should +HD-LD;LD;should +HD-NDX;NDX;should +HD-SAWADEE;SAWADEE;should +C1HDZ;C1HDZ;should +C1HDZ;Homeobox;should +C1HDZ;START;should not +C1HDZ;MEKHLA;should not +C2HDZ;C2HDZ;should +C2HDZ;Homeobox;should +C2HDZ;START;should not +C2HDZ;MEKHLA;should not
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tapscan_script_v74.pl Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,916 @@ +#!/usr/bin/perl +use strict; +use warnings; + +# Written by Gerrit Timmerhaus (gerrit.timmerhaus@biologie.uni-freiburg.de). +# Changes included by Kristian Ullrich, Per Wilhelmsson and Romy Petroll. + +# Script to extract all detected domains out of a hmmsearch results file and classify the families of all used proteins based on these domains. +# The classification depends on a table which contains all known classification rules for the protein families of interest and on specific coverage values defined for every domain. +# The script provides three outputs, namely output.1, output.2 and output.3. The output files are tables in ";"-delimited format. +# The structure of output.1 is: "sequence ID ; TAP family ; number of classifications ; domains". +# Output.3 shares in principle the same structure as output.1, except that subfamilies are considered. ("sequence ID ; TAP family ; Subfamily ; number of classifications ; domains") +# The superior TAP family is specified first, followed by the subfamily. If a TAP family has no subfamily, the TAP family is specified first and then a "-". +# The structure of output.2 is: "TAP family";"number of detected proteins". +# More than one entry for a protein is possible because the classification rules may allow more than one classification. +# +# The script must be startet with the arguments <hmmsearch output file> <classification rules> <output classifications file> <output family statistics file> <output subfamily classifications file> <"filter" if desired> + +if (!@ARGV or ($ARGV [0] eq "-h") or ($ARGV [0] eq "-help")) { + print "Usage: extract.and.classify.pl <hmmsearch output file> <classification rules> <output classifications file> <output family statistics file> <output subfamily classifications file> <\"filter\" (if desired)>\n\n"; + exit; +} + +# hmmsearch_output: domtblout file +my $hmmsearch_output = $ARGV [0]; +# decision_table: rules file +my $decision_table = $ARGV [1]; +# family_classifications: output.1 +my $family_classifications = $ARGV [2]; +# family_statistics: output.2 +my $family_statistics = $ARGV [3]; +# subfamily_classifications: output.3 +my $subfamily_classifications = $ARGV [4]; +# domspec_cuts: coverage values file +my $domspec_cuts = $ARGV [5]; +# gene_model_filte: filter for ARATH and ORYSA +my $gene_model_filter = $ARGV [6]; + +if ($family_statistics eq "") { + print "Usage: extract.and.classify.pl <hmmsearch output file> <classification rules> <output classifications file> <output family statistics file> <output subfamil classifications file> <\"filter\" (if desired)>\n\n"; + exit; +} + +if ($gene_model_filter and $gene_model_filter eq "filter") { + print "\nGene model filter is activated. It only works for TAIR (Arabidopsis) and TIGR (Rice) proteins up to now\n"; +} + +# Array where the $hmmsearch-output/domtblout file will be stored +my @output = (); +# Array with domain-specific coverage values +my @cuts = (); +# Array with rules +my @dec_table = (); +# Counter for the number of detected domains in the hmmsearch output file +my $entry_counter = 0; +# Containes the actual result for a query sequence +my $akt_entry = ""; +# Used to define query entry to ignore similar domains +my $whole_entry = ""; +# Includes the final entries after ignoring similar domains +my @results_of_extraction = (); +# Used to define query entry to ignore similar domains +my $extracted_domain = ""; +# Used to define query entry to ignore similar domains +my $present = ""; +# Used to define query entry to ignore similar domains +my $protein = ""; + +my $lek = ""; + +############################################ +### 1. Read in the hmmsearch output file ### +############################################ + +print "\n*** reading in $hmmsearch_output ***\n\n"; + +@output = get_file_data("$hmmsearch_output"); + +print "*** Parsing $hmmsearch_output ***\n\n"; + +# If wrong format exit the program, ninth row from the end +if ($output [-9] !~ /^# Program: hmmsearch.*/) { + print "wrong file format. $hmmsearch_output has to be a hmmsearch output file\n\n"; + exit; +} + +### 1.1 Trim input file (hmmsearch) to fit script. Erase lines starting with #. +@output = grep !/^#.*/, @output; + +# Read domain-specific coverage values file +# domname cut-off +# EIN3 0.5874125874 + +open(FILENAME,$domspec_cuts); +my %cuts = map { chomp; split /\t/ } <FILENAME>; + +### 1.2 Use first (prot name), third (domain name), eight (dom evalue), fourteenth (# overlapping) and fifteenth (# domain) column of data input(array element). +foreach my $output (@output) { + $output =~ /^(\S+)\s+\S+\s+\S+\s+(\S+)\s+\S+\s+(\S+)\s+\S+\s+\S+\s+\S+\s+\S+\s+\S+\s+(\S+)\s+\S+\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+.*/; + $lek = (($6-$5)+1)/$3; + $output = $1."\t".$2."\t".$4."\t".$lek; + } + +# Stucture: +# ARATHwo_AT2G03430.1 Ank 1.8e-08 0.93 +# ARATHwo_AT2G03430.1 Ank 7.8e-10 0.95 +# ARATHwo_AT2G03430.1 Ank 1.2e-08 0.93 +# ARATHwo_AT2G03430.1 Ank 1.2e-08 0.63 + +################################ +### 2. Modify the input data ### +################################ + +### 2.1 Run through length_cut_off loop (throws out entries below coverage values) +my @cutarray; +foreach my $output (@output) { + $output =~ /^(\S+)\t+(\S+)\t(\S+)\t(\S+)/; + if ($4 > $cuts{$2}) { + push (@cutarray, $output); + } +} + +# End up with (prot name) (domain name) (dom evalue) (overlapping) + +# Structure: +# ARATHwo_AT2G03430.1 Ank 1.8e-08 0.93 +# ARATHwo_AT2G03430.1 Ank 7.8e-10 0.95 +# ARATHwo_AT2G03430.1 Ank 1.2e-08 0.93 + + +### 2.2 Retain hit with the lowest evalue + +# String for the previous protein name +my $old_prot = ""; +# Coordinates for array filling +my $x = -1; +my $y = 1; +my $scoreeval; +my @newarray; +foreach my $cutarray (@cutarray) { + $cutarray =~ /^(\S+)\t+(\S+)\t(\S+)\t(\S+)/; + $akt_entry = $1.$2; + # If the protein contains more than one of the same domain + if ($old_prot eq $akt_entry) { + if ($3>$scoreeval) { + $y++; + $newarray[$x] = $1."\t".$2."\t".$3."\t".$y; + } + else { + $y++; + $newarray[$x] = $1."\t".$2."\t".$scoreeval."\t".$y; + } + } + # If the protein is new reset $j and push the entry in the next array + else { + $y = 1; + $x++; + $newarray[$x] = $1."\t".$2."\t".$3."\t".$y; + } + $old_prot = $akt_entry; + $scoreeval = $3 + +} +# Final structure: +# ARATHwo_AT2G03430.1 Ank 7.8e-10 3 + +### 2.3 Then sort @output primarily on prot name($1) and then secondarily on eval($3) to fit the "similar"-procedure. +@newarray = sort { (split '\t', $a)[0] cmp (split '\t', $b)[0] || (split '\t', $a)[2] <=> (split '\t', $b)[2] } @newarray; + +# Flags for merged domains +my $ARR_B_counter = 0; +my $bZIP_counter = 0; +my $ET_counter = 0; + +# Flags for similar domain omission +my $WD40domain = 0; # modified for FIE rule +my $FIEdomain = 0; # --------'''''-------- +my $MYBdomain = 0; +my $G2domain = 0; +my $YBdomain = 0; +my $YCdomain = 0; +my $PHDdomain = 0; +my $Alfindomain = 0; +my $Dr1domain = 0; +my $GATAdomain = 0; +my $Dofdomain = 0; +my $C1domain = 0; + +foreach my $line (@newarray) { + +# Reset the flags for domain omission + if ($line !~ /^$protein.*$/ ) { + $WD40domain = 0; + $FIEdomain = 0; + $MYBdomain = 0; + $G2domain = 0; + $YBdomain = 0; + $YCdomain = 0; + $PHDdomain = 0; + $Alfindomain = 0; + $Dr1domain = 0; + $GATAdomain = 0; + $Dofdomain = 0; + $C1domain = 0; + + } + + # Get the Query entry + $line =~ /^(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/; + $protein = "$1"; + $extracted_domain = $2; + $present = "$4"; # FIX! number for many Myb + $whole_entry = $1.";".$2.";".$3.";".$4; + + + # Ignore FIE if it is not better scored than WD40 + if ($extracted_domain eq "WD40" and $FIEdomain == 0) {$WD40domain = 1;} + + # Ignore similar domains part 1 + + if ($extracted_domain eq "Myb_DNA-binding" and $G2domain == 0) {$MYBdomain = 1;} + if ($extracted_domain eq "G2-like_Domain" and $MYBdomain == 0) {$G2domain = 1;} + + if ($extracted_domain eq "PHD" and $Alfindomain == 0 and $C1domain == 0) {$PHDdomain = 1;} + if ($extracted_domain eq "Alfin-like" and $PHDdomain == 0 and $C1domain == 0) {$Alfindomain = 1;} + if ($extracted_domain eq "C1_2" and $Alfindomain == 0 and $PHDdomain == 0) {$C1domain = 1;} + + if ($extracted_domain eq "GATA" and $Dofdomain == 0) {$GATAdomain = 1;} + if ($extracted_domain eq "zf-Dof" and $GATAdomain == 0) {$Dofdomain = 1;} + + # Ignore similar domains part 2 + + if ($extracted_domain eq "FIE_clipped_for_HMM" and $WD40domain == 1) {next;} + if ($extracted_domain eq "Myb_DNA-binding" and $G2domain == 1) {next;} + if ($extracted_domain eq "G2-like_Domain" and $MYBdomain == 1) {next;} + + if ($extracted_domain eq "PHD" and ($Alfindomain == 1 or $C1domain == 1)) {next;} + if ($extracted_domain eq "Alfin-like" and ($PHDdomain == 1 or $C1domain ==1)) {next;} + if ($extracted_domain eq "C1_2" and ($PHDdomain == 1 or $Alfindomain == 1)) {next;} + + if ($extracted_domain eq "GATA" and $Dofdomain == 1) {next;} + if ($extracted_domain eq "zf-Dof" and $GATAdomain == 1) {next;} + + # Save the entry + push @results_of_extraction, $whole_entry; + $entry_counter++; +} +print "$entry_counter domain matches were found in $hmmsearch_output\n\n"; + +# Sort the entries: +my @sorted_results_of_extraction = @results_of_extraction; +print "*** Preparing hmmsearch results for classification ***\n\n"; + +# This array will be filled with the sorted proteins +my $array_of_arrays = []; + +# String for the previous protein name +my $old_entry = ""; + +# Coordinates for array filling +my $i = -1; +my $j = 0; + +@dec_table = get_file_data("$decision_table"); + +if ($dec_table [0] !~ /^[^;]+;[^;]+;[^;]/) { + print "wrong file format. $decision_table has to be in the format family;domain;type\n\n"; + exit; +} + +# Formate the entries +foreach my $domain_entry (@sorted_results_of_extraction) { + + $domain_entry =~ /^([^;]+);/; + $akt_entry = $1; + # If the protein has more than one domain push it behind the old entry + if ($old_entry eq $akt_entry) { + $j++; + $array_of_arrays->[$i][$j] = $domain_entry; + } + + # If the protein is new reset $j and push the entry in the next array + else { + $j = 0; + $i++; + $array_of_arrays->[$i][$j] = $domain_entry; + } + $old_entry = $akt_entry; + +} +$old_entry = ""; + +# Remember the number of entries to use them later +my $number_of_proteins = $i+1; +print "$number_of_proteins proteins were found in $hmmsearch_output\n"; +$number_of_proteins--; + +######################################### +### 3. Get the classification entries ### +######################################### + +# Make a list of all families in the classification rules file +my @liste_alle_familien = ('0_no_family_found'); + +my $ARR_B_in_list = 0; +my $bZIP_in_list = 0; +my $ET_in_list = 0; + +my $array_of_classifications = []; +$i = -1; +$j = 0; + +my @classifications_input = get_file_data("$decision_table"); + +# Sort the classifications +my @classifications = sort @classifications_input; + + +foreach my $classification (@classifications) { + $classification =~ /^([^;]+);/; + $akt_entry = $1; + + # If the family has more than one domain entry push it behind the old entry + if ($old_entry eq $akt_entry) { + + $j++; + $array_of_classifications->[$i][$j] = $classification; + } + + # If the family is new push it in the next array and reset $j + else { + $j = 0; + $i++; + $array_of_classifications->[$i][$j] = $classification; + # Add family to list of all families in the classification rules + push(@liste_alle_familien,$akt_entry); + + # Add hardcoded "OR" defined families to list of families if subfamily is present + # Do this just once for each family + if (($akt_entry eq "GARP_ARR-B_Myb") or ($akt_entry eq "GARP_ARR-B_G2")) { + if ($ARR_B_in_list == 0) { + push(@liste_alle_familien,"GARP_ARR-B"); + $ARR_B_in_list++; + }} + if (($akt_entry eq "bZIP1") or ($akt_entry eq "bZIP2") or ($akt_entry eq "bZIPAUREO") or ($akt_entry eq "bZIPCDD")) { + if ($bZIP_in_list == 0) { + push(@liste_alle_familien,"bZIP"); + $bZIP_in_list++; + }} + if (($akt_entry eq "HRT") or ($akt_entry eq "GIY_YIG")) { + if ($ET_in_list == 0) { + push(@liste_alle_familien,"ET"); + $ET_in_list++; + }} + } + $old_entry = $akt_entry; +} + +my $number_of_classifications = $i+1; +print "$number_of_classifications different families were found in $decision_table\n\n"; +$number_of_classifications--; + +print "*** begin classification ***\n\n"; + +# Create the output files (output.1 and output.3) +my $outputfile = "$family_classifications"; +my $subfamilyoutput = "$subfamily_classifications"; + +unless (open(FAMILY_CLASSIFICATIONS, ">$outputfile")) { + print "Cannot open file \"$outputfile\" to write to!!\n\n"; + exit; +} + +unless (open(SUBFAMILY_CLASSIFICATIONS, ">$subfamilyoutput")) { + print "Cannot open file \"$subfamilyoutput\" to write to!!\n\n"; + exit; +} + +print FAMILY_CLASSIFICATIONS "family classifications for $hmmsearch_output\n"; +print SUBFAMILY_CLASSIFICATIONS "family classifications for $hmmsearch_output\n"; + +# Counter for classified families +my $classified_families = 0; +# Counter for unclassified families +my $unclassified_families = 0; +# Counter for the different amount of entries in classification possibilities +my @entry_counter = []; +# Counter for the possible classifications for the current protein +my $possibilities = 0; +# For later famliy statistics +my @family_list = []; +# This array will be filled with the result. Later written in output.1 and output.3 +my @family_classifications_file = []; +my @subfamily_classifications_file = []; + +# One time for every protein in array of arrays +for (my $ii = 0; $ii <= $number_of_proteins; $ii++) { + my $jj = 0; + $possibilities = 0; + + # Flag to mark that any family was found for the current protein + my $member_found = 0; + + # In this string all the domains of one protein will be stored + # The ; signs are important to mark a single entry (for example to differ AP2 and TF_AP2) + my $domains = ";"; + + # Get the domains in all elements of the current protein array: + while ($array_of_arrays->[$ii][$jj]) { + + $array_of_arrays->[$ii][$jj] =~ /([^;]+);[^;]+;(\d+)$/; + # For discrimination between MYB and MYB-related. If more than 1 MYB domains was found replaced "Myb_DNA-binding" with "two_or_more_Myb_DNA-binding" + # Same for LIM + if ($1 eq "Myb_DNA-binding" and $2 eq 2) { + $domains .= "MYB-2R;Myb_DNA-binding;"; + } + elsif ($1 eq "Myb_DNA-binding" and $2 eq 3) { + $domains .= "MYB-3R;Myb_DNA-binding;"; + } + elsif ($1 eq "Myb_DNA-binding" and $2 eq 4) { + $domains .= "MYB-4R;Myb_DNA-binding;"; + } + elsif ($1 eq "LIM" and $2 > 1) { + $domains .= "two_or_more_LIM;LIM;"; + } + else { + $domains .= "$1;"; + } + $jj++; + } + +######################### +### 4. Classification ### +######################### + +# Get the name of the actual protein for the output after classification: +$array_of_arrays->[$ii][0] =~ /^([^;]+);/; +$akt_entry = $1; + +# Variables for double classification decisions: + +# Here the shoulds for the actual family will be stored +my $shoulds_act_family = ";"; +# Here the shoulds of the previous CLASSIFIED family will be stored +my $shoulds_prev_family = ""; + +# This loop will be made for every family entry in output.1 and output.3 +for ($i = 0;$i <= $number_of_classifications; $i++) { + $j = 0; + + # Flag to mark if the protein is a possible member of the current family + my $possible_member = 0; + + # Reset the domain entries for the current family + $shoulds_act_family = ";"; + + # This loop will be made for every domain classification entry for the current family + while ($array_of_classifications->[$i][$j]) { + + # For "should" entries + if ($array_of_classifications->[$i][$j] =~ /;should$/) { + + # Get the domain from the classification entry + $array_of_classifications->[$i][$j] =~ /^[^;]+;([^;]+);[^;]+/; + + my $domain_classification_entry = $1; + # Check if the protein might be a member of the current family. If yes check the next entry + if ($domains =~ /;$domain_classification_entry;/) { + $possible_member = 1; + $j++; + # All matched domains for the current family are stored here: + $shoulds_act_family .= "$domain_classification_entry;"; + next; + } + $possible_member = 0; + $j++; + last; + } + + # For "should not" entries + if ($array_of_classifications->[$i][$j] =~ /;should not$/) { + + # Get the domain from the classification entry + $array_of_classifications->[$i][$j] =~ /^[^;]+;([^;]+);[^;]+/; + + # Check if the protein has the domain. If yes go to the next family + my $domain_check = $1; + if ($domains =~ /;$domain_check;/) { + $possible_member = 0; + last; + } + # If not check the next entry + else { + $j++; + next; + } + } + + # This happens when there is no "should" or "should not" entry at the right position + print "error in classification entry $array_of_classifications->[$i][$j]\n" + } + + # If the check was successful print the found classification + if ($possible_member) { + + $array_of_classifications->[$i][0] =~ /^([^;]+);/; + my $currentFamily = $1; + $possibilities++; + # The double classified proteins solution: + if ($possibilities > 1) { + + # Fill an array with the current domains in $domains + my @domains_array = []; + my $number_of_domains = ($domains =~ tr/;/;/)-1; + my $tempdomains = $domains; + for (my $domain_counter = 0;$domain_counter != $number_of_domains;$domain_counter++) { + + $tempdomains =~ /^;([^;]+);/; + $domains_array[$domain_counter] = $1; + # Remove the inserted domain from domains + $tempdomains =~ s/^;[^;]+;/;/; + } + # Compare the two possible families (actual and previous family) + my $array_counter = 0; + while ($array_counter < $number_of_domains) { + if ($shoulds_prev_family =~ /$domains_array[$array_counter]/) { + if ($shoulds_act_family =~ /$domains_array[$array_counter]/) { + + # If actual domain is in previous and actual family go to next domain + $array_counter++; + next; + + } + # Protein was right classified in the previous classification. The actual entry will be ignored + + $possibilities--; + last; + } + + # If actual domain is in the actual family but not in the previous. The actual classification is right + if ($shoulds_act_family =~ /$domains_array[$array_counter]/) { + + pop @family_classifications_file; + pop @subfamily_classifications_file; + + # Merge "GARBP_ARR-B_Myb" and "GARP_ARR-B_G2" to "GARP_ARR-B" + if ( ($currentFamily eq "GARP_ARR-B_Myb") or ($currentFamily eq "GARP_ARR-B_G2")){ + + push @family_classifications_file, "$akt_entry;GARP_ARR-B;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;GARP_ARR-B;-;$possibilities$domains\n"; + if ($ARR_B_counter == 0) { + push @liste_alle_familien, 'GARP_ARR-B'; + $ARR_B_counter++; + } + } + elsif ( ($currentFamily eq "bZIP1") or ($currentFamily eq "bZIP2") or ($currentFamily eq "bZIPAUREO") or ($currentFamily eq "bZIPCDD")) { + + push @family_classifications_file, "$akt_entry;bZIP;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;bZIP;-;$possibilities$domains\n"; + if ($bZIP_counter == 0) { + push @liste_alle_familien, 'bZIP'; + $bZIP_counter++; + } + } + elsif ( ($currentFamily eq "HRT") or ($currentFamily eq "GIY_YIG")) { + push @subfamily_classifications_file, "$akt_entry;ET;-;$possibilities$domains\n"; + push @family_classifications_file, "$akt_entry;ET;$possibilities$domains\n"; + if ($ET_counter == 0) { + push @liste_alle_familien, 'ET'; + $ET_counter++; + } + } + else { + + push @family_classifications_file, "$akt_entry;$currentFamily;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;$currentFamily;-;$possibilities$domains\n"; + } + $possibilities--; + last; + + } + + $array_counter++; + } + + } + # Normal enty: if the family is classified for one family write it in the array + else { + + # Store the matched domains of the actual family for respectively later doubly classification decision + $shoulds_prev_family = $shoulds_act_family; + # Store the entry in the output array + # Merge "GARBP_ARR-B_Myb" and "GARP_ARR-B_G2" to "GARP_ARR-B" + # Add GARP_ARR-B to @liste_alle_familien ONCE + + if ( ($currentFamily eq "GARP_ARR-B_Myb") or ($currentFamily eq "GARP_ARR-B_G2")){ + push @family_classifications_file, "$akt_entry;GARP_ARR-B;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;GARP_ARR-B;-;$possibilities$domains\n"; + if ($ARR_B_counter == 0) { + push @liste_alle_familien, 'GARP_ARR-B'; + $ARR_B_counter++; + } + } + elsif (($currentFamily eq "bZIP1") or ($currentFamily eq "bZIP2") or ($currentFamily eq "bZIPAUREO") or ($currentFamily eq "bZIPCDD")){ + push @family_classifications_file, "$akt_entry;bZIP;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;bZIP;-;$possibilities$domains\n"; + if ($bZIP_counter == 0) { + push @liste_alle_familien, 'bZIP'; + $bZIP_counter++; + } + } + elsif (($currentFamily eq "HRT") or ($currentFamily eq "GIY_YIG")){ + push @family_classifications_file, "$akt_entry;ET;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;ET;-;$possibilities$domains\n"; + if ($ET_counter == 0) { + push @liste_alle_familien, 'ET'; + $ET_counter++; + } + } + else { + + push @family_classifications_file, "$akt_entry;$currentFamily;$possibilities$domains\n"; + push @subfamily_classifications_file, "$akt_entry;$currentFamily;-;$possibilities$domains\n"; + } + + $classified_families++; + } + + # Flag to mark that any family was found for the current protein + $member_found = 1; + } + + } + + if ($possibilities == 0) { + + # If no family was found for the current protein give this out + push @family_classifications_file, "$akt_entry;0_no_family_found;0$domains\n"; + push @subfamily_classifications_file, "$akt_entry;0_no_family_found;-;0$domains\n"; + $unclassified_families++; + } + + # Increment the number of possible members in the array at the right position + else { + + my $temp = $entry_counter[$possibilities]; + $temp++; + $entry_counter[$possibilities] = $temp; + } + +} +# Define for which TAP families also subfamilies may be present: + +foreach my $subfamily_classifications_file (@subfamily_classifications_file) { + + # If "C2H2" was assigned to a sequence, write "C2H2;C2H2" to output.3 as TAP family and subfamily, respectively. + if ($subfamily_classifications_file =~ /C2H2/ ) { + $subfamily_classifications_file =~ s/C2H2;-/C2H2;C2H2/ + } + # If "C2H2_IDD" was assigned to a sequence, write "C2H2;C2H2_IDD" to output.3 as TAP family and subfamily, respectively. + if ($subfamily_classifications_file =~ /C2H2_IDD/ ) { + $subfamily_classifications_file =~ s/C2H2_IDD;-/C2H2;C2H2_IDD/ + } + if ($subfamily_classifications_file =~ /bHLH/ ) { + $subfamily_classifications_file =~ s/bHLH;-/bHLH;bHLH/ + } + if ($subfamily_classifications_file =~ /bHLH_TCP/ ) { + $subfamily_classifications_file =~ s/bHLH_TCP;-/bHLH;bHLH_TCP/ + } + if ($subfamily_classifications_file =~ /NF-YA/ ) { + $subfamily_classifications_file =~ s/NF-YA;-/NFY;NF-YA/ + } + if ($subfamily_classifications_file =~ /NF-YB/ ) { + $subfamily_classifications_file =~ s/NF-YB;-/NFY;NF-YB/ + } + if ($subfamily_classifications_file =~ /NF-YC/ ) { + $subfamily_classifications_file =~ s/NF-YC;-/NFY;NF-YC/ + } + if ($subfamily_classifications_file =~ /C1HDZ/ ) { + $subfamily_classifications_file =~ s/C1HDZ;-/HDZ;C1HDZ/ + } + if ($subfamily_classifications_file =~ /C2HDZ/ ) { + $subfamily_classifications_file =~ s/C2HDZ;-/HDZ;C2HDZ/ + } + if ($subfamily_classifications_file =~ /C3HDZ/ ) { + $subfamily_classifications_file =~ s/C3HDZ;-/HDZ;C3HDZ/ + } + if ($subfamily_classifications_file =~ /C4HDZ/ ) { + $subfamily_classifications_file =~ s/C4HDZ;-/HDZ;C4HDZ/ + } + if ($subfamily_classifications_file =~ /MYST/ ) { + $subfamily_classifications_file =~ s/MYST;-/HAT;MYST/ + } + if ($subfamily_classifications_file =~ /CBP/ ) { + $subfamily_classifications_file =~ s/CBP;-/HAT;CBP/ + } + if ($subfamily_classifications_file =~ /TAFII250/ ) { + $subfamily_classifications_file =~ s/TAFII250;-/HAT;TAFII250/ + } + if ($subfamily_classifications_file =~ /GNAT/ ) { + $subfamily_classifications_file =~ s/GNAT;-/HAT;GNAT/ + } + if ($subfamily_classifications_file =~ /LOB1/ ) { + $subfamily_classifications_file =~ s/LOB1;-/LBD;LOB1/ + } + if ($subfamily_classifications_file =~ /LOB2/ ) { + $subfamily_classifications_file =~ s/LOB2;-/LBD;LOB2/ + } + if ($subfamily_classifications_file =~ /MYB-related/ ) { + $subfamily_classifications_file =~ s/MYB-related;-/MYB;MYB-related/ + } + if ($subfamily_classifications_file =~ /SWI\/SNF_SWI3/ ) { + $subfamily_classifications_file =~ s/SWI\/SNF_SWI3;-/MYB-related;SWI\/SNF_SWI3/ + } + if ($subfamily_classifications_file =~ /MYB-2R/ ) { + $subfamily_classifications_file =~ s/MYB-2R;-/MYB;MYB-2R/ + } + if ($subfamily_classifications_file =~ /MYB-3R/ ) { + $subfamily_classifications_file =~ s/MYB-3R;-/MYB;MYB-3R/ + } + if ($subfamily_classifications_file =~ /MYB-4R/ ) { + $subfamily_classifications_file =~ s/MYB-4R;-/MYB;MYB-4R/ + } + if ($subfamily_classifications_file =~ /RKD/ ) { + $subfamily_classifications_file =~ s/RKD;-/RWP-RK;RKD/ + } + if ($subfamily_classifications_file =~ /NLP/ ) { + $subfamily_classifications_file =~ s/NLP;-/RWP-RK;NLP/ + } + if ($subfamily_classifications_file =~ /AP2/ ) { + $subfamily_classifications_file =~ s/AP2;-/AP2;AP2/ + } + if ($subfamily_classifications_file =~ /CRF/ ) { + $subfamily_classifications_file =~ s/CRF;-/AP2;CRF/ + } +} + +### 4.1 Filter for gene models + +# Filter for gene models. All alternative splice variants will be removed from Arabidopsis and rice entries +shift @family_classifications_file; +shift @subfamily_classifications_file; +# Sort the array for filter and a sorted output +@family_classifications_file = sort @family_classifications_file; +@subfamily_classifications_file = sort @subfamily_classifications_file; +if ($gene_model_filter and $gene_model_filter eq "filter") { + print "*** filtering gene models: removing splice variants ***\n"; + + # temporary array for the filtered entries + my @filtered_entries = []; + my $models_removed_counter = 0; + + my $prev_entry = ""; + foreach my $entry_line (@family_classifications_file) { + if ($entry_line !~ /^([^.]+).\d+/) { + print "Bad format for filtering splice variants. Use \"filter\" only for TAIR (Aradopsis) or TIGR (Rice) files.\n\n"; + exit; + } + if ($1 eq $prev_entry) { + $models_removed_counter++; + next; + } + $prev_entry = $1; + push @filtered_entries, $entry_line; + } + print " -> $models_removed_counter gene models were removed\n\n"; + @family_classifications_file = @filtered_entries; + shift @family_classifications_file; +} + +# Print the results in the array to output.1 and output.3 and push the names of the TAP families and Subfamilies into $family_list to create output.2 +foreach my $fcf_line (@family_classifications_file) { + $fcf_line =~ /^[^;]+;([^;]+)/; + #push @family_list, "$1"; + print FAMILY_CLASSIFICATIONS "$fcf_line"; +} +close (FAMILY_CLASSIFICATIONS); + +foreach my $fcf_line (@subfamily_classifications_file) { + $fcf_line =~ /^[^;]+;([^;]+);([^;]+)/; + if ($2 eq "-") { + push @family_list, "$1"; + } + else { + push @family_list, "$2"; + } + print SUBFAMILY_CLASSIFICATIONS "$fcf_line"; +} + +close (SUBFAMILY_CLASSIFICATIONS); + +print "*** calculating the family statistics and write it in $family_statistics ***\n\n"; + +################################## +### 5. Create the output files ### +################################## + +my $statistics_outputfile = "$family_statistics"; + +unless (open(FAMILY_STATISTICS, ">$statistics_outputfile")) { + print "Cannot open file \"$statistics_outputfile\" to write to!!\n\n"; + exit; +} + +# Count the family entries +my @output_family_statistics = (); +my @gefundene_familien = (); +my $family_counter = 1; + +shift @family_list; +@family_list = sort @family_list; + + +my $old_family = ""; +push @family_list, 'BAD FIX'; # Makes to loop go through every fam and stops at non fam(BAD FIX). + +foreach my $line (@family_list) { + if ($line eq $old_family) { + $family_counter++; + } + elsif ($old_family ne "") { + push (@output_family_statistics,"$old_family;$family_counter\n"); + # Add all found families to a list of found families + push (@gefundene_familien,"$old_family"); + $family_counter=1; + } + $old_family = $line; +} + +my %hash = (); + +# Put all families from the classifictaion ruled and all found families in a hash +foreach my $element (@gefundene_familien,@liste_alle_familien) {$hash{$element}++;} + +# Remove merged families +delete $hash{'GARP_ARR-B_Myb'}; +delete $hash{'GARP_ARR-B_G2'}; +delete $hash{'bZIP1'}; +delete $hash{'bZIP2'}; +delete $hash{'bZIPAUREO'}; +delete $hash{'bZIPCDD'}; +delete $hash{'HRT'}; +delete $hash{'GIY_YIG'}; + +# Add all not found families from the classification rules to @output_family_statistics +# With zero as number of families found + +foreach my $element (keys %hash) { + if (($hash{$element} == 1) and ( $element eq "0_no_family_found")) { + push (@output_family_statistics,"$element;$unclassified_families\n"); + } + if (($hash{$element} == 1) and ($element ne "0_no_family_found")) { + push (@output_family_statistics,"$element;0\n"); + } +} + +# Sort @output_family_statistics caseinsensitive-alphabetically +my @sortierte_statistik = sort {lc $a cmp lc $b} @output_family_statistics; + +# Print headline to @sortierte_statistik +unshift (@sortierte_statistik,"family statistics for $hmmsearch_output\n"); + +print FAMILY_STATISTICS @sortierte_statistik; + +# Print FAMILY_STATISTICS "$old_family;$family_counter\n"; + +close (FAMILY_STATISTICS); + +################################################# +### 6. Give out some statistical informations ### +################################################# + +$entry_counter [0] = 0; +my $sum = 0; +foreach my $entry (@entry_counter) { + $sum += $entry; +} +print "$classified_families classifications were found for $sum proteins.\n"; +print "This classifications are divided in:\n"; +my $count = 0; +foreach my $element (@entry_counter) { + if ($count != 0) { + print "$element proteins were classified for $count"; + if ($count == 1) {print " family\n";} + else {print " different families\n";} + } + $count++; +} +print "\n$unclassified_families proteins could not be classified\n\n"; + +print "*** The results were written in $family_classifications and $subfamily_classifications ***\n"; +print "*** done ***\n\n"; + +exit; + +sub get_file_data { + + my ($filename) = @_; + + use strict; + use warnings; + + my @filedata = (); + + unless( open(GET_FILE_DATA, $filename)) { + print STDERR "Cannot open file \"$filename\"n\n"; + exit; + } + + @filedata = <GET_FILE_DATA>; + + close GET_FILE_DATA; + + return @filedata; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/PUBLIC_Ectocarpus-sp7_proteins_head.fa Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,32 @@ +>Ec-00_001160.1 Forkhead-associated (FHA) domain (809) ;mRNA; f:1452084-1459465 +MEPPPPVPAPIIVAGTRAKETGVDGNSSAATAGSTAAVSPPLDKAQQSALDPAKDGSASLPPPPPKTAGSLNGDGVKAVASGGYKPPSWGLTEAPGASGLSLTVLKGGVEVGSISLDNRTHVLLGRQQGVVDVLLEHPSISRKHAILQHGQNGALFLFDNGSTHGCSVNKKKIPPKEFHRLHVGDVIKFGESTRLYALEGPEELRPAEYESDNLRNLRLDAGRKQLAAKLAKIKAGGAGEGGGKGGGDSGEYGISWGFDEDAVAEEEDEDGDGAERDEDEVELPDYLKTEAQKRRRRDTKIGLTEDNVHKRDAKLFEKLQLKLTKIEEIEETIRSKNKARERGKEGGEGGSGGEEGGRKAGRGTEDGEEDDDYYDRTAPVVPTSSGTSKSDLASKKAEIKARRFGARDKTKKLSAVTPPADDTAAAERKRGVGPGQNEAAAQSLEALTRRGEAVVEDLERTQAGLAELEAEEAGEAALAEDGGVAADPLDMFMTENRRKERQQAIVRLTAKREALREEQALLKVMVEAARPSMPTLKKSPAPAAATASVAATKEETVPVPETTTRGSGSSSSDQAAEPDDRKDTRAAAVDREHMPGGGGYGEAMPGSTGSGSHEKYGDEVEKTAPSRAVLAPEAASTLGSMPSPVAPPCRSNAIPEARKSPAGTVASAGTRERGVEKGQAETRHPGVKEGKEASGTKKRGTPVGTSMLPPPPSKRQQRRAENSKAGLAGENDDDPVEPKAKRTVKGPAMPPPLGKPSKVGGGVSVPTAVARVGEKVAGKEALEGGDVDWVPPKDALEKMAALNRKFGY* +>Ec-00_001310.1 SET domain protein (668) ;mRNA; r:1563326-1571582 +MAIPTSKDEDLLDDEQAPAVAAAQDQELVAGSDAGTVAPSKKKKKKKKKKKSPQQNNQLAQVQVYETIKDLDTFKVTEDSVSGRCVIASRDLKAGELVLREPPFVKVVRRDCASRQCAYCCQQVTERGKIEADVPFAVYCSRACQAREDALRAAEASALGKLAGISAARDVDIDLLRMLLRLLITRAKALGLREPSGDSDSVSRGVDEEGEDGTMGEGLFLRQQWENLYALMHHREAMAPDWISVVREAGEDLLQLLPEWVRFDVEEVVQLACRVNVNAHGLRDDSGANLVIGVGMFPLTAMINHACRPNCTFVYFGGNLEVRTLEPVSAGAELSVYYIDLLQSTAARRQELLTSKHFLCKCSRCENPSSMDDYLDGVCCTDCGERGCLTPTPPPSAEDILAAQLAQLGEGLADESAANGSGKGMGKKKGGGGGSTTSRRETSVNPSAALGAKDGSGNGSRAGGVRGGSESVVQKVYCSACGREYPGTAVEESVARAKALWDAAMAVVRAKSFSLARKSLEKWLQDYDAGAVLPPTATKKKSVKKDRLKLHPANVMVVQTLVPLSNCCTFEEDHAASARHLRRAVSAMEAVYPANFPELGDFHAALADANDALLQKRGQTLPKKSRSQAVSERKQALERAATIRSVCLGKDHPATREAARALDRVTG* +>Ec-00_001360.1 aureochrome 2 (442) ;mRNA; r:1593684-1598236 +MPASVKPPVFTSMVHRKVQHHTSWQDADFQPDDLGLDLTDLSTMTGFLMNEVPDNTGHFYPPWANELSPLVKDEPSAFMMPDPAAPRQPQPRQQQQQQDQRLPAPEGLPAAPVADPALDVIMGGATGSSRPGSTTSSSSSGASSMLRAPKAAAAAGAAALGGVGSTFRKTPAGGVTRRRSSSKEEQAKKRRERNRVLARRTRLRKKFFFQSLQQQVARLQRENERLKGIVTTRCPDSVGEILMSCRSKMPSMVADCAGQATAVLDQSGFLLVKALQSSQPSFCVTDPQMPDNPIVYASDTFIELTGYDRAQVLGRNCRFLQGPDTDPDAVAKIRKGIEEGSDTSVYLRQYKADGTVFWNHVFVAALRNSEHKIINYVGIQHPLDKEPSPEVVACINNGKEQEIMSVQEEDRPAGWGGQWPEDVNGDLATLDHLMAGGWGTD* +>Ec-00_001730.1 WD40 repeat (1089) ;mRNA; r:1948051-1960039 +MSGYSRAGVGVGGGGGGGGGGGGGGGGGGGGANSNRGASASRGPGGALQGAHMAGTVATVAAAGGQQSKSVAVFRSLELCDLLKLEIGNITSEMGQHLEEREEYEKKFRQQLAEMDRIQHSLKQLQEAHMVMKQQYEDEIVRLRHQLESSHPLKSDPDSRGGGAGSSVQGVSGGMSMGSTGPSQSHSSGPPHKGLGGMPVGGPHGLGSPRKQPGQLEPGAASLRGPMLTAPGAQLVGGGGMGAGSGVVLPGMAMRGRGHGEDDYGGGRGQGPNGQTLEPLSKRPRLADLQGPPGPLPPHVLHPPHGFSQSQHQQHPQQHQQHPPPHHKGGFFGGRPGAPPAAVAGGGGGGGGGGDNGGRRPGGPEGQYRWPGGSGPGGGGGGGDSGGAGGGGGAGVPAMGPGKQRQLGPGGGGAPPRELSSSTGKGGGAGGGGGGNGNGMVEIPDAIPAELSYQVRYEAEEGEGGRAAAAKDGLAVELAKSQDLRSVVCCVRFSTDGTKIAAGSHSCVKVFDVNSFKQLYVCRKQVQEEQGAPQTADGGDPYVRAVCFSPDGLSIVAGMEKNSAKVLVLEEEGGRQGAITLSGHESEVYSLDWVSDMIASGSGDGRIRLWDSVTGACKASLGDMGGPQDGVTSVVLRQDTTMVAAASIDRVVHVWSTQTHKILHRLDGHSESVYAITFSADGNRLVSGGLDKTIKVWDLGPGSEGRLSPQAKTLPGGHKDYVLSVCFSQDGKYIISGSKDRSVTMWDARLMKRVATITGFKNSVIGVSASPFNSMFATGSGDNLVCVWNYGDRDDYRRQSESAVAPASSGRSSSAEKPRPPSRSSRSASPPISNASSLRKDDRASSPSPPRASPGRGGRGGDGDGGSSRGRGGVGLVEKGRGGSAGAGAGATNGKGRTSPSSDEREGRTAAGSEERKRSGSKESSRTGPRGGGGGHDSDRKQQRHESKHRQGRAGSPSDGRRSSSNSKERSSGRGGGVEPMDEDDDTEEEKQGEEAEGGGGKRAGNTDETPLGRKRGSGGVGNGSGGRKALNRSSSSSSSSSGSGGGGGGGGGGGGGGGRGPSPEKRGKAALPPKKELQRRDSPSSAE* +>Ec-00_010160.1 Phytochrome-like protein (978) ;mRNA; f:17223525-17226527 +MDPHTHTNERIHMRSSGGCDGDMRTASTHIQSCGCAFAIEETADDMYPSGLRILGVSQNAVEAPWAYASSVSDLLGKDLGHLLRVECVRTVRSLVTRYAQACKPSHEDDDHISPKANRITADACPSPRIRGELRPGAASGAEGDIASFTVTGSNPGVYLVDVERHGSDCARVEHTPGLLLLGDLLESIPVGSNPVESTAALCDALAKSMPAYDRVMVYRFAPDGSGQDDGSGEVVHESVRAGADIGSSYLNLRFPALDIPPIARKLFKLVGVRFIADTSAPAVPMITLHDQASSPLDLFRSALRAPAECHLRYLRNMGVKASLVVSIAVDGGTWGLFSFHSYTRTVHPSCEERLSVEMAASVVSSLISRYQREEIAATALSLSRTLGNLGNYTRVNDFLSADHHSLLGILDVDAVILCEHLRSVTLYGKKDITLSLEECQELRNGDGDESSEMAISFRTLGARGVAFFWVRSFFVAFLRGSIANSVKWAGNPDAPVNKDEVMTPRASFELFMRTSGARCKAWSPLTVDLLNMVRQGFSSQLYAEALPADLQETFARVSHELRTPFHGVMGALEILEAGNGIMGAEEQLDVIRSALRCGGSMMSTLNDILEIAKDRNNTEVVRGRFSASGPIVLAVAAMGLFAAAESVELTVEIGPPDDVLEVTGDMRRIKAIVQNLVNNAIKFTPSGGKVRISLVVFDSLQEVTDWWAKETGRFGAQTWMASSGGESAPGTGSSQAIKWHVYCVEDSGIGVLPADLPHLVEAYRQISHGASRSHAGTGLGLHICRTHMEAMCGSFGIASTFSEKDTSGGTLFAVVLPLDSEEPGAAVNTQEPLEAEVNLRLLDHKIRRFFKDNGADVEVMSATDGFIALEMHEAARRNQSHGSVLAGMFIDFHMPDLDGIECTKRIRLLEADNGWSRIMICGCTADPTQAIRRVFQNAGGDEVISKPWCPGQVESICNAMVANVLNAEQKSGGDGGA* +>Ec-00_001700.1 NIN-like transcription factor (579) ;mRNA; r:1923852-1928158 +MHTPLLNSHHRDTGLLETLEKPDKEPSTSRQQDTVMSLSTHSGTRGRPPANRASCGNSTIPDITVRSPADGLPPRLAIPHEPPKVPRRTPSKGNKRETLADIAKRIPVGLMRHYFNYPLRAAAEAMDISVTTLKRLCRRHGVKRWPHRQICGINRTLNDLETQHDTAKGDEVDSVADQLRQLYRRRDVIIELAFESDDESISNGSDDAPPAGRKLSKRGSFTSSDGGGSSASSSSFPSPPSSPSPSRSSSFVSPPASPTPTDERTHTAAGVTWLNNGAPGTGLPSLAPVLSVTVSSGGGGGGGGSSSSRSSKCPAVRSRGSSRGAKISSGGVGSSGSSGSYGKTSSGSMRPHRPSGTSGSGRSKGPRKPRSSAMFGSIPGLGKVPPPVVTVATGAAVSTTLSLPERTSASTSTSTSPVRADASTETSSRVTIPVPSRPLSASIPANSSSTHDGGGTSTVARLSLIGSGASWISDNDAPSPSTTSVDILGLGPHGSQKMDDLAILGDLLFGPDDAAVAAAAAAGVRAQSSAATYPPPSSTTSPPSSSYERGFSAGLASASNWGGGLCWDANGLSHLGPL* +>Ec-01_005970.1 Acetyltransferase (GNAT) domain protein (608) ;mRNA; f:5180258-5197307 +MERSVSLAKRRNLVIGWRKMARKKEIQVVLLAALISTLGIGPVAARIVVTDHLAFVASTVRTRTSSSSSPSTFQLQRDVRKTAAVPPRRYEQQQQTSRLRPRCGWGLERPSLSPDPVPDAFSAPLRRRRRRPRQQQQQHLVVAASSGLQEPLERVASIPVQEAAAGRGDRVSYGRGGTVVEAAAPEGVDAADGGGHERKTAGSVYGDEVGVLRGAVIWKGEENDADKDDGGGGGEDDDFEGDDIDIYLPGGMGAGLRPGAKPRRKRPKPPVMWSKPVRGSKRLAVVLAGREDLDEAGALCIKVFFGQPDSPWKAAQLRQLLHEQRQDLESRCSRRESVMFKAIDTRRKDGMVGFVEVSETSGSKYGMGAGITLADTRPVVSNLAVDPRVRRCGVGSALMEACEDLVKTWSFDEIILQVEEANEAALAFYGGKGFKKLFVDKAARRYDTSGFLLQNVRTSKLTMRKALGSSTKRSKSGDSLADHVSQLFGFLSKPFVVRPPPSKRLPARSTASVETLFGRPGAVGSGSGVRRPRSGGSGRRVTSSGGRVGSSSGGIGRSSGRSEAAAGGHDTSFSSSSSSSSSSSSSSSSAGEISTSRTRLPRTIRRR* +>Ec-01_004960.2 Ankyrin repeat-containing domain (797) ;mRNA; f:4270310-4278760 +MPAEGRRANPIELNIYIEVVLYSSTTRLAAASFFFCFSFFSSKTTVSGEPTPRGGDFADPPPIFMIARSTASPSLIDITKLSSAMMLGNMSSSVLSQLDRLKESFLWAAAKGGRVEECESLLEMGTDINWVSPEGDTPLLAACRNGHLQTALCLLSHGASANQVDKDGRTALHVSCRYGKEAVAEALILRGADVSARDHGGATAFESQGPHVPDGMLQRLDTIAQRSGRRTHHHLGISTAAINTSTFTTSSGRSAQDGGGSMSSRRSNTVLDSDTLSARATGGGVPENNIGDGRSGWRWRSNTGRIEHHVPLTDDAHDGAAGVGPAAISTGDGGRDDAQREERGGVEGRTTFPSISRKAGSTSADSRPSSRVHSAGRAASSSGDGVRSSDDEHDSAAIAAAIASASAPAPATTVTAAGSTIYPMMASNTLASPAPGPTSLRSPRRPLHFHGTPTASWDGAGADGGAGHGITVSNSRDSGSGGNRWNGSGSASGSSSSSSSSSSSNSGRGGGGERQSRPAASATNANSPATATGTTATTTTVTPTCGLYGSLSTNAVAERINSRSITAHWLSPPAAAEVRGPLGDHAMHWGGRTNGDLVSGDAPYESIAPAGAGGRGVLGGGGHGGDPAQGDGGSSDHVGLIRALWLQEKESRLASEAKADAVRDQSARLWNELVRSEKQVLSLEEGAAALREERDKLKATLNGEGFKGRSLQELEHLEKDLRKALEGVCGERDRIVQQQLAKEEQRLCVVCQENERSVLLLPCRHLCVCRGCSERQELTLCPLCRDHITESLVVYS* +>Ec-01_005360.1 Ferrochelatase-2, putative chloroplast precursor (583) ;mRNA; f:4601493-4607037 +MRREQKMRAAATLLCTLSVASAWFVAPPTSAASSHRHFDSTTSRASSPPFSSTAQRRRRNAAGGDAADASTTSGNIRGRRRAAAARTSHRMVASSRNWEWGRAKSSGDAPTPAAQAFGSDQELKLGVLLLNLGGPERPEDVQPFLFNLFADPDIIRLPKLVQWLQNPIAAVLAARRAPQSKSAYESIGGGSPIVSWTNAQAKGIASQLEAKGLSGTKCYVGMRYWHPFTEAALEAVEDDEINALVILPLYPQFSISTSGSSLRILNEEFTRRPEQWGHKNVVHTVVPSYHDRPGYVNAMASLIAREVAEYTPEQRMQGVQVLFSAHGVPKSYIDAGDPYKAQIESCVKLISEKVDGINAEGGPGAKPGSSGAAAGGVTYHLSYQSRVGPVEWLQPYTDAKIHELADNGCKNLVVVPVSFVSEHIETLEEIDMEYREVAEEAGITNWRRVPALNTDPAFIEDMADMVVEALALPTLTVSEAFTRNNCDRKEAEGFLEKALDGMYGMPKTGGKPPKSGKVGGAGAAGANASSGGGGGADGEGSGDKRKEAARVLSTLSGAAFAADGVGREIAGLFTATSDGIFF* +>Ec-01_005370.1 Hypothetical protein (68) ;mRNA; r:4610968-4611171 +MAHGVLLSSTTGWWRNLQVWCPCCCIFQFEPVGFYLHAKRGPSRLQTISLVGQRIMIFADPKQLESA* +>Ec-01_005380.1 Aminoalcohol phosphotransferase (396) ;mRNA; r:4612890-4621264 +MSDFRGGSVVTPRAVKYLRRYQYHGSDRSLLYKYVLSPLAETCLVFLPSWMAPNLVTTIGLGLTTASYLLLYLSMPGLVSNESTPWWVFPAAAAGLIVYQTLDNMDGKQARRTGSSSPLGLIFDHGCDAINCCFGVVFVSCILDAGSSLPLLAAIVLNQLVPFFFTTWEHYYTHELILPIVNGPSEGVVLGAVSACLRGVYGPEFFSAPREGLRGWPLGEVLMACSLLGVALTVFKQIVLVARARRLAGRGMVNPIRDASWFVVLMVLGGSWASVRPELFLTKPYTMVLLFGLLHVDMAVHLMVCHVCNMVCRSFRPILIPFMLVAANSFFPGGPLLGERTLVLGFTVLTFIYEALYLYLVVTETSLALDIYVFKLGKRSDATRNGGGTTTSKKD* +>Ec-01_005390.1 WD40 repeat (407) ;mRNA; f:4621750-4629304 +MSESKDGGGSSGGGGRAPYLDVSEDKMYDSRLRGDYELHDYAHSPSNQRLKEQRESKGGGSIADRVLGARLHPGLTISENSAEVFVTRFAPEGNLLAAACGDGTIRIFHVSTGRLAYNLQSSSSQSLPTTSLRFRPAIAQSKTRNVLVSCNAGGEIQHWHITSGKCLSTIKEEDQFYALDYRKDGSVFAATGKNHTVHIYDETTKHEISLLQGGSGYGSSSAPGHSNRVFAVKFHPEDPEVLLTAGWDNTVQFWDMRVGHSVRSIFGPHISGDSLDICGNEILTGSWRPNDPLEIWDYGTAELKETIPWNRSSAQGVQPTLLYTAQFSNTADGGRYIAAGGSGANEAKIFDHTANNQLVGTVTGLPRGVFTVDFSPSPDAKKVAVAGGDASIRIIDIVEEKTADVY* +>Ec-01_005400.1 hypothetical protein (1653) ;mRNA; r:4630582-4637601 +MSYVERLREELMGSLGFEVESQSCLSGNRGSAGERDQKARRRRSGGGAIAAVQVGVPYPHALRPRQIQQQQQQQQQQQGRRAADDSSCMVGGGCLRPQRADLSRVPRSRRRPASASSVRYLANKSNMNDGYSGANKDMDVLLKICRDQGRANSLNPADFDSSISAIANLPVKQFMTAAKDGWGPYFRNLHSNMRKETTARSTTTAKRMNRPPPERRRRSSQQQHPPLPLPRRASSGATKDEKKNNRRYSKVANDLHHMVRSADVASSDSNNRRRRRQSEDDNSDGLAASRARAAKRQAEEMVIRRMGVRLQALWQELKIPDPDRAYVTAAYLDAGGGGGGFSGQRARGGGGIEEPPVGSGGGGGGNSGPTSENVNRELTRQIRLLLEHRAATVKVLRCVNARESRLCEVEQALQAFQWHRGLDTDALGVVVSLAGLRDASLDVVRAVEEWRTNLWQPRAFCWRGVNYLEKMWQDTIFLRSPVGQSLLASVGLEAKDMTFVLYPSGNGGGAPSAIPGDSSTCGSELSDGERDTKGDVLRPPPAIVRLTPKAELVSAYRARYFEQSPVVGAGGRPGGESGDVGGERGRLKSTAAAATAAMEAVVVQEAALQRRVEAERVRLATKGCFVPLLRWLPGDSRNPAPKGKRRWSSSSPPSVPPAARNPPPGAESETHAAKDVVDSFAHNDGRSGVPAVRDVSVSAAIPPDPLTGGELPLQLQQQQQQQQQQQQPATVTNTSTSAGASLFRGSVRDSAQQERLTVSGEGDGHEEEDREDEDGDEGEDREQDEKGQASLLHTEGVRSQKQPERVEEGGHDQTQSQSEGPQAGQGSTAEVEGRSAFHENGVGEQTKMTNPLAPSFRLETSLHPPGAEEESRSETDIARRSACTVGAEARAGENGCLKVVATSECGVGSAGIPVGGVSSGDDGDGDTDKPSIVPGELRQEGIVDEERSVAPSVCYDDDFEEDAEGYSLSDSETSGEGGSADFGECLNGSISSDDAKDARRDQGVAHTTEEEGAAISCGEGRQNPVVGGNNVLTEHQRQQQLVEERETVAATAALRIQRSARAWLCRGGRLRNTKPDEKEAPQESAKEGKQEWLDAWRARRHSQVRQAQRPEEIPEEAAIRIQTTVRGGLASRQVEELRRELHARKMAATLKIQAFVRRRRARRNTTAINRREISENVDLRQQQEAEPDPELAGVGVCRGTARQGEEEILAVSSEFDSAEAATRSDRSARTIQKAVRRRLSSGNTTSKADGLPAMEIKGQPKTLEEEAQGQPDAGKRGDGPGAVLVNRCGSGDGSFIPSQQAPVVAAESSSPVVRAGDELDSPSSPAAVDGTPSSTPASGADGLDARLSDADGKNNRPDVAESAFVEAGLPSTSLPSEDSAESETADDAAGGDLSQTLQHPREGSNFQPEGASSWSPLPSPTTPRPEAAERGQQENDHQPRQVEGGGEAAVSLPPSSATACPPTLGPTVDVFVGLLPVDTDTSADAMMTTVGGTHAGGGAESLGFSAVGLSPVPERQALVAVAASASEYGSEDPDFDALSSSPGGSSRGGVSRGEEENAGVAARVATRSSVVTTGEKVEEPTESRGVERVAAPVAAAVAGAATEKNLAGDIEPVSLAAAEVATTKEVERGGTQCTEDGSDNSSIVSLSVSSGS* +>Ec-01_005410.1 Hypothetical protein (154) ;mRNA; f:4637945-4638913 +MGSAWSRKQATIHSSFQAAPLAVASWDLQHVKILHDRFRDDGYEFGIDRVSLRDLISSALPLARGVTGDLWEIYVEHEQEMLYPLELFAAVALQCHGLSRCGVSTIYSCPVFLCKNCKNPSYFHAHFGPCRSEGRRLDWCISSKASPGIRVHF* +>Ec-01_005410.2 Hypothetical protein (200) ;mRNA; f:4641160-4641759 +MATGDSKTEFSIPTPIPVNPCNGTCSLSSPVDASCKIVVTQHLYAVHAAIVLRGNSKQLKTRDLVRICARFDGKPVATLPPQSFHAHSSPAHCVVRLKKAGALHNRTAVPPSCRWPLATPLVYLFDQGSNHRHLRRQQHRSKRLLVCIAPTTIKPARPTTVVPSLRQPPPRLCPAPRPAPLVPATHHPRAIHQQRGRRC* +>Ec-01_005410.3 Hypothetical protein (102) ;mRNA; f:4639028-4639623 +MRGALLKLGVYRWLSVQRKRENSTKTRPQPTDSVCAVEVYYRVGVSFRAPRLLRSQPTEKRRGLCSTCLTLKADKRFRTTSWPSASRPSWEHFRRSLDPAV*
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.1.tsv Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,10 @@ +sequence ID TAP family number of classifications domains +Ec-00_001160.1 FHA 1 FHA; +Ec-00_001310.1 SET 1 SET; +Ec-00_001360.1 bZIP 2 bZIP_AUREO;bZIP_2;bZIP_CDD; +Ec-00_001700.1 RKD 1 RWP-RK; +Ec-00_001730.1 0_no_family_found 0 WD40; +Ec-00_010160.1 0_no_family_found 0 Response_reg; +Ec-01_004960.2 0_no_family_found 0 K-box; +Ec-01_005390.1 0_no_family_found 0 WD40; +Ec-01_005970.1 GNAT 1 Acetyltransf_1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.2.tsv Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,139 @@ +TAP family number of detected proteins +0_no_family_found 4 +ABI3/VP1 0 +ADA2 0 +Alfin-like 0 +ALOG 0 +AP2 0 +ARF 0 +Argonaute 0 +ARID 0 +Aux/IAA 0 +BBR/BPC 0 +BES1 0 +bHLH 0 +bHLH_TCP 0 +bHSH 0 +BSD domain containing 0 +bZIP 1 +C1HDZ 0 +C2C2_CO-like 0 +C2C2_Dof 0 +C2C2_GATA 0 +C2C2_YABBY 0 +C2H2 0 +C2H2_IDD 0 +C2HDZ 0 +C3H 0 +C3HDZ 0 +C4HDZ 0 +CAMTA 0 +CBP 0 +Coactivator p15 0 +CPP 0 +CRF 0 +CSD 0 +CudA 0 +DBP 0 +DDT 0 +Dicer 0 +DUF246 domain containing/O-FucT 0 +DUF296 domain containing 0 +DUF547 domain containing 0 +DUF632 domain containing 0 +DUF833 domain containing/TANGO2 0 +E2F/DP 0 +EIL 0 +ET 0 +FHA 1 +GARP_ARR-B 0 +GARP_G2-like 0 +GeBP 0 +GIF 0 +GNAT 1 +GRAS 0 +GRF 0 +HD-LD 0 +HD-NDX 0 +HD-other 0 +HD-SAWADEE 0 +HD_BEL 0 +HD_DDT 0 +HD_KNOX1 0 +HD_KNOX2 0 +HD_PHD 0 +HD_PINTOX 0 +HD_PLINC 0 +HD_WOX 0 +HMG 0 +HSF 0 +IWS1 0 +Jumonji_Other 0 +Jumonji_PKDM7 0 +LDL/FLD 0 +LFY 0 +LIM 0 +LOB1 0 +LOB2 0 +LUG 0 +MADS 0 +MADS_MIKC 0 +MBF1 0 +Med6 0 +Med7 0 +mTERF 0 +MYB-2R 0 +MYB-3R 0 +MYB-4R 0 +MYB-related 0 +MYST 0 +NAC 0 +NF-YA 0 +NF-YB 0 +NF-YC 0 +NLP 0 +NZZ 0 +OFP 0 +PcG_EZ 0 +PcG_FIE 0 +PcG_MSI 0 +PcG_VEFS 0 +PHD 0 +PLATZ 0 +Pseudo ARR-B 0 +RB 0 +Rcd1-like 0 +Rel 0 +RF-X 0 +RKD 1 +RRN3 0 +Runt 0 +S1Fa-like 0 +SAP 0 +SBP 0 +SET 1 +Sigma70-like 0 +Sin3 0 +Sir2 0 +SOH1 0 +SRS 0 +SWI/SNF_BAF60b 0 +SWI/SNF_SNF2 0 +SWI/SNF_SWI3 0 +TAFII250 0 +TEA 0 +TFb2 0 +tify 0 +TRAF 0 +Trihelix 0 +TUB 0 +ULT 0 +VARL 0 +VOZ 0 +Whirly 0 +WRKY 0 +Zinc finger, AN1 and A20 type 0 +Zinc finger, MIZ type 0 +Zinc finger, ZPR1 0 +Zn_clus 0 +ZPR 0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.3.tsv Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,10 @@ +sequence ID TAP family subfamily number of classifications domains +Ec-00_001160.1 FHA - 1 FHA; +Ec-00_001310.1 SET - 1 SET; +Ec-00_001360.1 bZIP - 2 bZIP_AUREO;bZIP_2;bZIP_CDD; +Ec-00_001700.1 RWP-RK RKD 1 RWP-RK; +Ec-00_001730.1 0_no_family_found - 0 WD40; +Ec-00_010160.1 0_no_family_found - 0 Response_reg; +Ec-01_004960.2 0_no_family_found - 0 K-box; +Ec-01_005390.1 0_no_family_found - 0 WD40; +Ec-01_005970.1 HAT GNAT 1 Acetyltransf_1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.domtbl.tsv Wed Feb 14 13:54:16 2024 +0000 @@ -0,0 +1,28 @@ +# --- full sequence --- -------------- this domain ------------- hmm coord ali coord env coord +# target name accession tlen query name accession qlen E-value score bias # of c-Evalue i-Evalue score bias from to from to from to acc description of target +#------------------- ---------- ----- -------------------- ---------- ----- --------- ------ ----- --- --- --------- --------- ------ ----- ----- ----- ----- ----- ----- ----- ---- --------------------- +Ec-01_005970.1 - 608 Acetyltransf_1 PF00583.26 117 1.6e-14 44.3 0.0 1 1 1.8e-15 2.8e-14 43.5 0.0 31 117 348 434 305 434 0.66 Acetyltransferase (GNAT) domain protein (608) ;mRNA; f:5180258-5197307 +Ec-00_001160.1 - 809 FHA PF00498.27 69 2.3e-19 59.7 0.2 1 1 2.4e-20 3.8e-19 59.0 0.2 2 69 123 190 122 190 0.96 Forkhead-associated (FHA) domain (809) ;mRNA; f:1452084-1459465 +Ec-01_004960.2 - 797 K-box PF01486.18 93 3.1e-06 17.5 4.8 1 1 3.7e-07 6e-06 16.6 4.8 11 85 683 754 674 756 0.81 Ankyrin repeat-containing domain (797) ;mRNA; f:4270310-4278760 +Ec-00_010160.1 - 978 Response_reg PF00072.25 112 1.2e-10 31.7 0.0 1 1 1.6e-11 2.6e-10 30.7 0.0 13 106 845 952 840 958 0.74 Phytochrome-like protein (978) ;mRNA; f:17223525-17226527 +Ec-00_001700.1 - 579 RWP-RK PF02042.16 49 2.1e-19 59.5 0.1 1 1 2.3e-20 3.7e-19 58.7 0.1 7 46 111 150 106 153 0.95 NIN-like transcription factor (579) ;mRNA; r:1923852-1928158 +Ec-00_001310.1 - 668 SET PF00856.29 169 4.1e-18 56.4 0.0 1 1 1.2e-18 2e-17 54.2 0.0 1 168 83 338 83 339 0.72 SET domain protein (668) ;mRNA; r:1563326-1571582 +Ec-00_001730.1 - 1089 WD40 PF00400.33 38 1.9e-38 120.1 18.0 1 4 3.9e-09 3.1e-08 24.5 0.2 5 38 579 611 575 611 0.86 WD40 repeat (1089) ;mRNA; r:1948051-1960039 +Ec-00_001730.1 - 1089 WD40 PF00400.33 38 1.9e-38 120.1 18.0 2 4 1e-11 8.1e-11 32.7 0.3 3 38 662 698 660 698 0.93 WD40 repeat (1089) ;mRNA; r:1948051-1960039 +Ec-00_001730.1 - 1089 WD40 PF00400.33 38 1.9e-38 120.1 18.0 3 4 9.6e-11 7.7e-10 29.6 0.1 8 38 716 747 709 747 0.88 WD40 repeat (1089) ;mRNA; r:1948051-1960039 +Ec-00_001730.1 - 1089 WD40 PF00400.33 38 1.9e-38 120.1 18.0 4 4 1.4e-05 0.00011 13.3 0.0 4 38 754 789 751 789 0.90 WD40 repeat (1089) ;mRNA; r:1948051-1960039 +Ec-01_005390.1 - 407 WD40 PF00400.33 38 3.8e-17 52.7 3.9 1 2 7.3e-06 5.8e-05 14.1 0.0 17 37 88 108 72 109 0.77 WD40 repeat (407) ;mRNA; f:4621750-4629304 +Ec-01_005390.1 - 407 WD40 PF00400.33 38 3.8e-17 52.7 3.9 2 2 1.8e-07 1.4e-06 19.2 0.1 8 38 223 255 217 255 0.89 WD40 repeat (407) ;mRNA; f:4621750-4629304 +Ec-00_001360.1 - 442 bZIP_2 - 55 9.2e-17 51.1 12.5 1 1 1e-17 1.6e-16 50.3 12.5 2 49 183 230 182 233 0.92 aureochrome 2 (442) ;mRNA; r:1593684-1598236 +Ec-00_001360.1 - 442 bZIP_CDD - 52 4.1e-13 39.4 7.8 1 1 4.4e-14 7.1e-13 38.6 7.8 2 43 188 229 187 232 0.94 aureochrome 2 (442) ;mRNA; r:1593684-1598236 +Ec-00_001360.1 - 442 bZIP_AUREO - 52 2.2e-22 69.4 6.2 1 1 2.7e-23 4.4e-22 68.5 6.2 1 51 187 237 187 238 0.97 aureochrome 2 (442) ;mRNA; r:1593684-1598236 +# +# Program: hmmsearch +# Version: 3.3.2 (Nov 2020) +# Pipeline mode: SEARCH +# Query file: /home/saskia/code/github/galaxyproject/tools-iuc/tools/tapscan/tapscan_domains_v12.txt +# Target file: /tmp/saskia/tmp4n93kzh1/files/b/7/b/dataset_b7b39fd6-41f1-440d-bfb7-da7a3a9f070e.dat +# Option settings: hmmsearch --domtblout domtblout.txt --cut_ga /home/saskia/code/github/galaxyproject/tools-iuc/tools/tapscan/tapscan_domains_v12.txt /tmp/saskia/tmp4n93kzh1/files/b/7/b/dataset_b7b39fd6-41f1-440d-bfb7-da7a3a9f070e.dat +# Current dir: /tmp/saskia/tmp4n93kzh1/job_working_directory/000/4/working +# Date: Mon Nov 13 14:57:28 2023 +# [ok]