SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
how complete is the draft assembly? cegma? chrishah De novo discovery 12 02-12-2015 11:39 AM
CEGMA error flobpf Bioinformatics 12 01-27-2015 06:54 AM
IQ run eoh001 SOLiD 1 07-05-2011 07:07 AM
Consed would not run in 10.6 YongLi Bioinformatics 0 02-10-2011 09:21 AM
25 bp/frag run vs 50 bp/frag run regarding N2S ratio and cost efficiency smallcreek SOLiD 1 12-02-2009 05:47 AM

Reply
 
Thread Tools
Old 03-27-2012, 08:24 AM   #1
lukas1848
Member
 
Location: Germany

Join Date: Jun 2011
Posts: 54
Default Getting CEGMA to run

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
Quote:
geneid-train did not work properly
The cegma error output file only says
Quote:
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?
lukas1848 is offline   Reply With Quote
Old 04-16-2012, 06:53 AM   #2
SaMdeP
Junior Member
 
Location: Montreal

Join Date: Apr 2012
Posts: 2
Default

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

Last edited by SaMdeP; 04-16-2012 at 07:03 AM.
SaMdeP is offline   Reply With Quote
Old 04-18-2012, 04:29 PM   #3
kbradnam
Member
 
Location: Davis, CA

Join Date: May 2011
Posts: 53
Default

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

Code:
>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)

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

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

Code:
# 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
kbradnam is offline   Reply With Quote
Old 04-19-2012, 07:56 AM   #4
SaMdeP
Junior Member
 
Location: Montreal

Join Date: Apr 2012
Posts: 2
Default

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 :
Code:
# 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
Code:
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
Code:
    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
SaMdeP is offline   Reply With Quote
Old 04-19-2012, 08:08 AM   #5
kbradnam
Member
 
Location: Davis, CA

Join Date: May 2011
Posts: 53
Default

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.
kbradnam is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:28 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO