Getting CEGMA to run

[lukas1848]

Hi,

I try to run CEGMA on a genome assembly and after taking the first 10 obstacles in the installation, I am now stuck.


After running for ~2 h on the provided sample genome CEGMA crashes returning

geneid-train did not work properly

The cegma error output file only says

local.geneid.selected.gff file does not exist".

This local.geneid.selected.gff is a temporary file, that was definetly created during the run and should have been available to cegma.

Strangely, every application seems to work fine if I run them separately. Here's the versions I use:

- CEGMA 2.4,
- geneid 1.4.4
- genewise 2.4.1
- blast 2.2.22+

Any ideas what might cause that?

[SaMdeP]

Hi,

lukas1848, did you finally find a solution for this problem.

I have exactly the same error.
Same version of Geneid, genewise and CEGMA.

Two remarks:

1) my version of hmmer was compiled with the -disable-thread , so I had to edit the hmmselect script of CEGMA to remove the "--cpu " option.
system ("hmmsearch --cpu $THREADS $hmm_directory/$kog_id.hmm /tmp/prot$$.fa > /tmp/hmm$$.output") && die "Can't run hmmsearch\n

2) I -obviously- had problem with glib during genewise installation.
What I did is replacing "glib-config" by "pkg-config glib-2.0" in makefiles.
I don't know if it's orthodox, but it worked.

Any idea ?

SaM

[kbradnam]

Hi,

I think this is due to HMMER. For the latest version of CEGMA, you need HMMER version 3. You can test the problem like so

1) Copy the following protein to prot.fa

>KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA
MAENKKKNDDMATVEFESSEEVSIIPTFDKMGLREDLLRGIYAYGFEKPSAIQQRAIPAI
LKARDVIAQAQSGTGKTATFSISVLQSLDTQVRETQALILSPTRELAVQIQKVVLALGDY
MNVQCHACIGGTNLGEDIRKLDYGQHVVSGTPGRVFDMIRRRNLRTRAIKLLVLDEADEM
LNKGFKEQLYDIYRYLPPGAQVVLLSATLPHEILEMTSKFMTDPIRILVKRDELTLEGIK
QFFVAVDREEWKFDTLIDLYDTLTITQAVLFCNTRRKVDWLTDKMKEANFTVSSMHGDME
QKDRDEVMKEFRAGTTRVLISTDVWARGLDVPQVSLVINYDLPNNRELYIHRIGRSGRFG
RKGVAINFVKQDDVRILRDIEQYYSTQIDEMPMNIADII*

2) Now run the following HMMER command (where /path/to/CEGMA... is your correct path to your CEGMA installation)

hmmsearch /path/to/CEGMA/data/hmm_profiles/KOG0328.hmm prot.fa

3) This should produce valid HMMER output. E.g.

# hmmsearch :: search profile(s) against a sequence database
# HMMER 3.0 (March 2010); http://hmmer.org/
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query HMM file:                  Work/Bioinformatics_packages/CEGMA/current/data/hmm_profiles/KOG0328.hmm
# target sequence database:        a.fa
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       KOG0328  [M=394]
Scores for complete sequences (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Sequence                                           Description
    ------- ------ -----    ------- ------ -----   ---- --  --------                                           -----------
          0 1013.0   5.9          0 1012.8   4.1    1.0  1  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 


Domain annotation for each sequence (and alignments):
>> KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA  
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ! 1012.8   4.1         0         0       2     394 .]       7     399 ..       6     399 .. 1.00

  Alignments for each domain:
  == domain 1    score: 1012.8 bits;  conditional E-value: 0
                                             KOG0328   2 keeddekvefetseevevistFeemnlkedlLRGiYaYGfEkPsaiqqRaitqiikgrD 60 
                                                         k++d+++vefe+seev++i+tF++m+l+edlLRGiYaYGfEkPsaiqqRai +i+k+rD
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA   7 KNDDMATVEFESSEEVSIIPTFDKMGLREDLLRGIYAYGFEKPSAIQQRAIPAILKARD 65 
                                                         689******************************************************** PP

                                             KOG0328  61 viAqaqsGtGKtatfsisvlqslDlsvretqaLiLsPtRELAvqiqkvvlalGdymnvq 119
                                                         viAqaqsGtGKtatfsisvlqslD++vretqaLiLsPtRELAvqiqkvvlalGdymnvq
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA  66 VIAQAQSGTGKTATFSISVLQSLDTQVRETQALILSPTRELAVQIQKVVLALGDYMNVQ 124
                                                         *********************************************************** PP

                                             KOG0328 120 ahaciGGtslgeDirKldyGqhvvsGtPGRvlDmikrrsLrtRaiklLvLDEaDElLnk 178
                                                         +haciGGt+lgeDirKldyGqhvvsGtPGRv+Dmi+rr+LrtRaiklLvLDEaDE+Lnk
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 125 CHACIGGTNLGEDIRKLDYGQHVVSGTPGRVFDMIRRRNLRTRAIKLLVLDEADEMLNK 183
                                                         *********************************************************** PP

                                             KOG0328 179 GFKeqiYDiyryLPpatqvvlvsAtlpkeiLEmtsKFmtdPvriLvKRDEltLEGiKqf 237
                                                         GFKeq+YDiyryLPp++qvvl+sAtlp+eiLEmtsKFmtdP+riLvKRDEltLEGiKqf
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 184 GFKEQLYDIYRYLPPGAQVVLLSATLPHEILEMTSKFMTDPIRILVKRDELTLEGIKQF 242
                                                         *********************************************************** PP

                                             KOG0328 238 fvavekEEWKFDtLcDlYDtLtitqaviFCntkrKvDwLteklreanFtvssmHGdmpq 296
                                                         fvav++EEWKFDtL+DlYDtLtitqav+FCnt+rKvDwLt+k++eanFtvssmHGdm+q
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 243 FVAVDREEWKFDTLIDLYDTLTITQAVLFCNTRRKVDWLTDKMKEANFTVSSMHGDMEQ 301
                                                         *********************************************************** PP

                                             KOG0328 297 keRdeimkeFRsGesRvListDvWARGiDvqqvsLvinYDLPnnrElYiHRiGRsGRfG 355
                                                         k+Rde+mkeFR+G++RvListDvWARG+Dv+qvsLvinYDLPnnrElYiHRiGRsGRfG
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 302 KDRDEVMKEFRAGTTRVLISTDVWARGLDVPQVSLVINYDLPNNRELYIHRIGRSGRFG 360
                                                         *********************************************************** PP

                                             KOG0328 356 RKGvAinFvkkdDiriLRdiEqyYstqidemPmnvadli 394
                                                         RKGvAinFvk+dD+riLRdiEqyYstqidemPmn+ad+i
  KOG0328.7_1|geneid_v1.4_predicted_protein_1|400_AA 361 RKGVAINFVKQDDVRILRDIEQYYSTQIDEMPMNIADII 399
                                                         *************************************97 PP



Internal pipeline statistics summary:
-------------------------------------
Query model(s):                            1  (394 nodes)
Target sequences:                          1  (400 residues)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.04u 0.00s 00:00:00.04 Elapsed: 00:00:00.04
# Mc/sec: 3.94
//

If you don't get this output, something is probably wrong with your HMMER set up.

Regards,

Keith

[SaMdeP]

Hi,

Thank you Keith for your reply.

So, after further investigation, you pointed me to the right direction ! It was a problem with hmmer. Well, not hmmer itself, but the alias for hmmer in our system that points to version 2.3.2 instead of 3.0

Specifically, what happen is that hmm_select parses the output of hmmsearch and look for this line :

# match score from hmmsearch output
 if (/\s+\S+\s+([\d\.]+).*$full_kog_id/){

where $full_kog_id is 'KOG3497.9_1' in the following exemple.

In hmmer 2.3.2, the output of hmmsearch is

KOG3497.9_1|geneid_v1.4_predicted_protein_1|68_AA               173.7      5e-53   1

were KOG3497.9_1 is before the score

and in hmmer 3.0, the output is

    2.8e-56  175.3   1.8      3e-56  175.1   1.3    1.0  1  KOG3497.9_1|geneid_v1.4_predicted_protein_1|68_AA

were KOG3497.9_1 is after the score

That's why hmm results were not found.

Thanks again for your help and thank you for the work you did on moving CEGMA to ncbi blast

SaM

[kbradnam]

Glad we fixed the problem. When I was upgrading CEGMA to work with NCBI BLAST rather than WU-BLAST/AB-BLAST I also took the decision to make CEGMA only work with the new version of HMMER (v3.0). This required several changes to the code due to the new output format of HMMER. It did cross my mind that I could try to keep backwards compatibility with HMMER v2.0 as well, but as my available time to work on CEGMA is so limited, I decided to only support the latest version.