As you have run CEGMA with the --ext option to keep intermediate files, could you try running the final step separately to show what the error message from geneid-train is? i.e. run:

geneid-train -v local.geneid.selected.gff local.geneid.selected.dna geneid_params

And then report the output.

Hmm, everything is looking good initially. The output you should be seeing for the gene-id step of CEGMA is something like:

RUNNING: geneid-train  local.geneid.selected.gff local.geneid.selected.dna geneid_params 2>test.cegma.errors
DATA COLLECTED: 15 Coding sequences containing 48 introns
RUNNING: make_paramfile /Korflab/Packages/CEGMA/current/data/self.param.template \
            geneid_params/coding.initial.5.logs geneid_params/coding.transition.5.logs \
            geneid_params/start.logs geneid_params/acc.logs geneid_params/don.logs \
            geneid_params/intron.max > geneid_params/self.param

And when I run the separate geneid-train step, I see the following:

geneid-train -v local.geneid.selected.gff local.geneid.selected.dna geneid_params
DATA COLLECTED: 15 Coding sequences containing 48 introns
Intron model
Coding model
don
false_don
acc
false_acc
start
false_start

This is with geneid v1.4...I wonder if they changed something in v1.4.4 that breaks CEGMA?

I don't have much left to suggest, but let's check that the required input files for geneid-train look okay. Can you run the following please:

ls -l local.geneid.selected.gff local.geneid.selected.dna geneid_params
head local.geneid.selected.gff local.geneid.selected.dna

And then print the output here.

I'm stumped, so I can only suggest the following steps to at least try to narrow down where exactly the error is occurring.

1. Make a copy of the geneid-train script. E.g.

   cp geneid-train geneid-train2
   

2. Add __END__ after line 92. E.g.

   my $intron_initial = new geneid::SequenceModel('intron', 'FREQ', 4,
            \@intron, 10, 0);
   
   __END__
   

3. Then test the script as before (but using your modified code):

   geneid-train2 -v local.geneid.selected.gff local.geneid.selected.dna geneid_params
   

4. Do you still see the same error message? If not, move the __END__ token to after line 97. E.g.

   my $intron_transition = new geneid::SequenceModel('intron', 'MM', 5,
            \@intron, 10, 0);
   
   __END__
   

5. I'm betting that one of these two lines is producing the error message (which means digging into a very poorly documented Perl module to troubleshoot further...I didn't write this code by the way!).

I really can't imagine why everything else would work but then fail at this step. Just to double check, you didn't interrupt the code part way the first time you ran it, and then restarted?

Hmm, I'm pretty stumped by this. Not sure why everything would work until that one specific part of the geneid-train script. The only thing I can think of is to really nail down as to which part of the script is causing the problem. My suggestion:

1. Create a copy of the geneid-train script (e.g. cp /path/to/geneid-train geneid-train2)

2. Add an __END__ token around line 94. E.g.

   my $intron_initial = new geneid::SequenceModel('intron', 'FREQ', 4,
            \@intron, 10, 0);
   
   __END__
   

3. Run this new script:

   ./geneid-train2 -v local.geneid.selected.gff local.geneid.selected.dna geneid_params
   

4. If this doesn't produce an error, move the __END__ token down a few lines. E.g.

   my $intron_transition = new geneid::SequenceModel('intron', 'MM', 5,
           \@intron, 10, 0);
   __END__
   

5. Save and run script again. I'm betting that one of these two lines will report the error. Then we have to try to fathom out what is happening in the not-very-easy-to-understand geneid Perl module.

Hmm, not for the first time an answer I posted to this thread failed to post (or was removed?). 2nd attempt...

I think I have an answer. To test it, you'll need to change one line of code in two different Perl module (it's the same line of code). Make backups of geneid.pm and HMMstar.pm (in CEGMA/lib directory). Then change one line in the dna_word_table subroutine (this is line 41 in HMMstar.pm) as follows:

OLD:        $code .= "$tab foreach my \$c$i qw(A C G T) {\n";
NEW:        $code .= "$tab foreach my \$c$i (qw(A C G T)) {\n";

The addition of the 2nd set of parentheses might make all the difference and remove your problem. Please let me know if this works.

Hi,

I'm using ncbi_blast+/2.2.28, hmmer/3.0, geneid 1.4.4, genewise 2.4.1, and perl 5.10.1. Everything works fine with the sample file, but I am getting the same error as above with my datafile which is a transcriptome of 46,474 ORFs. I applied the fix above, to geneid.pm and HMMstar.pm even though I am not using perl 5.14 (though I've tried the sample file using both perl 5.14 and 5.10 and it works with both even with the fix to geneid.pm and HMMstar.pm).

Here is a sample of my input file:

>m.10939
TGTCAGGgggTGGgggAGAGGTTTCCTGGCCGCTTCTGTGGGCCACCGTCGGTGTCGTCCTGCAGCTGCCCTCCTTTGCAGCCAGAGCCACTGGCACTGTCCCCGAGAAGGCCCTCCCCTCCCGTGGGCTCCCTCCTCTGCCTGCAGCTCTTTCGTTGtgtgAGGgggAGCCCGCCcccTGTGGAAGGCATGCACCTGACCCACATATTTCCTCTGAGCCCTGCCTTGCTCCTACAGTCCCCACCTGGggggggAGCAGAGGACCTTCCTCCCTctctTATTCTTtttCCTTCCTTCCTT
>m.27653
ATGTTTGCAGGAGATGAGCTGGCAGTAGTACAGAGCTGGGAAGCCCCTTCAGTGGCAGATTCAATGGAGACAGGCAAGGTTTGCACTATGGCCATGATGGTCACAGTAGAGACTTTGCAAGTGATCTCTGTTTATCCACTACAGGCAGGTGAGAAGAGTCCACTGGCCTTCAATGATGTGAAAAGGGGAGCCTCCCTCGCGAATATGAGCCTCAGGACAGCCGATTCCAGTGTCCGGCTCCCGTCGCCCACCCCATTCGACTGCCGCCCCAACTTCTCTTCATGCTCTGTCTCAGATTCA
>m.48642
ATGTTCAACAAATCACAATACTTCTCTGCTGTTATTGCATCCACAGAATGTCACCCAGTGGCCCTGTTTAGGGTTACCTGGACAGTTCTCAAGAAAATGAGCCCAGAGGCCCATCTGCAGGGCCATGTGGAGGAATTTTCCATGCATCTGGAGGAAaaaaTCACTCTGATTCATTCCCAGTTGGACGTGAAGCATACTGTCGAGCTGTCCAGGGAGCATATTTACCCAGTTGGATCTGAGGAAGTGGACAGGGTCCTCCAGAATTTAAACATTAGATCCATACTACTCCTGTCTGGTTAA
>m.2225
AACATCTGGTACGGCTTCACAACCACCCTGCTGTCCTCCGGCGTTGACGTGGTGCTGATCGCCGCCTCTTACGCTCTGATCCTCCGGGCCGTCTTCCAGCTCTCCTGCAAGGAAGCCCGGCGGAAAGCTCTCCGCACCTGCGGCTCCCACATCTCCGTCATCCTCATGTTCTACATGCCAGCCTTCTTCTCCTTCCTCACCCATCGCTTCGGCCACAACATCCCTAGCTGGATCCATATCCTGCTGGCCGACCTCTACGTGATCGTGCCGCCCATGCTCAACCCCATTGTCTACGGAGTC
>m.21342
ATGGCAGATTCCGGAGGgggAGATAGAGTTACATCAGGGGAGACTGATACACCCCGATCCACACTGGCTACAACTAGCCATGTGGAGATTGAGCGGGAGCCCCTGAgagagGCTCACCTGTCGCTGCGGGTGGTTTGGACCATAGAGGCCTTGCACAGGAAATGCACTACGAGGATCTATGATTCCTCCTGGACTTCCTTAGTGGCCCGGTTTCAAGCAGCCCAGGTAATGCCCACCAGGCCACAGTTCCACAGGTCCTTGACTTTGTCCAGGACGAGTTGGATAAGGGGCTCTCCTTAG

Here is the command issued:

cegma -v -ext -g /panfs/roc/groups/14/mcgaughs/mcgaughs/ref.fa

********************************************************************************
**                    MAPPING PROTEINS TO GENOME (TBLASTN)                    **
********************************************************************************

RUNNING: genome_map  -n genome -p 6 -o 5000 -c 2000 -t 1  -v  /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/kogs.fa /panfs/roc/groups/14/mcgaughs/mcgaughs/ref.fa 2>output.cegma.errors
Found 2287 candidate regions in /panfs/roc/groups/14/mcgaughs/mcgaughs/ref.fa

********************************************************************************
**     MAKING INITIAL GENE PREDICTIONS FOR CORE GENES (GENEWISE + GENEID)     **
********************************************************************************

RUNNING: local_map -n local -f -h /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/hmm_profiles -i KOG -v  genome.chunks.fa 2>output.cegma.errors
NOTE: created 1551 geneid predictions

********************************************************************************
**           FILTERING INITIAL PROTEINS PRODUCED BY GENEID (HMMER)            **
********************************************************************************

RUNNING: hmm_select -i KOG -o local -t 1 -v  /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/hmm_profiles local.geneid.fa /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/profiles_cutoff.tbl 2>output.cegma.errors
NOTE: Found 417 geneid predictions with scores above threshold value

********************************************************************************
**       CALCULATING GENEID PARAMETERS FROM SELECTED GENEID PREDICTIONS       **
********************************************************************************

RUNNING: geneid-train -v local.geneid.selected.gff local.geneid.selected.dna geneid_params 2>output.cegma.errors
RUNNING: make_paramfile /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/self.param.template \
        geneid_params/coding.initial.5.logs geneid_params/coding.transition.5.logs \
        geneid_params/start.logs geneid_params/acc.logs geneid_params/don.logs \
        geneid_params/intron.max >  geneid_params/self.param

********************************************************************************
**                           ACCURATE LOCAL MAPPING                           **
********************************************************************************

RUNNING: local_map -n local_self -g local.genewise.gff -d geneid_params/self.param -h /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/hmm_profiles -i KOG -v  genome.chunks.fa 2>output.cegma.errors
NOTE: created 0 geneid predictions

********************************************************************************
**                              FINAL FILTERING                               **
********************************************************************************

RUNNING: hmm_select -i KOG -o local_self -t 1 -v  /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/hmm_profiles local_self.geneid.fa /home/mcgaughs/mcgaughs/lcl/opt/cegma/current/data/profiles_cutoff.tbl 2>output.cegma.errors
NOTE: Foud 0 geneid predictions with scores above threshold value
sh: local_self.geneid.selected.gff: No such file or directory
Can't run sed

mcgaughs@labq03 [/panfs/roc/groups/14/mcgaughs/mcgaughs/reftwo] % geneid-train -v local.geneid.selected.gff local.geneid.selected.dna geneid_params 2>output.cegma.errors
DATA COLLECTED: 417 Coding sequences containing 679 introns 
Intron model
Coding model
don
false_don
acc
false_acc
start
false_start

mcgaughs@labq03 [/panfs/roc/groups/14/mcgaughs/mcgaughs/reftwo] % ls -l local.geneid.selected.gff local.geneid.selected.dna geneid_params
-rw-rw-r--. 1 mcgaughs mcgaughs 563370 Apr  4 12:12 local.geneid.selected.dna
-rw-rw-r--. 1 mcgaughs mcgaughs 147601 Apr  4 12:12 local.geneid.selected.gff

geneid_params:
total 1104
-rw-rw-r--. 1 mcgaughs mcgaughs    981 Apr  4 13:06 acc.logs
-rw-rw-r--. 1 mcgaughs mcgaughs  57501 Apr  4 13:06 coding.initial.5.logs
-rw-rw-r--. 1 mcgaughs mcgaughs 249017 Apr  4 13:06 coding.transition.5.logs
-rw-rw-r--. 1 mcgaughs mcgaughs    567 Apr  4 13:06 don.logs
-rw-rw-r--. 1 mcgaughs mcgaughs     76 Apr  4 13:06 intron.max
-rw-rw-r--. 1 mcgaughs mcgaughs 310082 Apr  4 12:12 self.param
-rw-rw-r--. 1 mcgaughs mcgaughs    290 Apr  4 13:06 start.logs

mcgaughs@labq03 [/panfs/roc/groups/14/mcgaughs/mcgaughs/reftwo] % head local.geneid.selected.gff local.geneid.selected.dna
==> local.geneid.selected.gff <==
# Sequence KOG0003.2 - Length = 471 bps
KOG0003.2    geneid_v1.4    First    1    376    159.74    +    0    KOG0003.2_1
KOG0003.2    geneid_v1.4    Exon    1    376    159.74    +    .    KOG0003.2_1
KOG0003.2    geneid_v1.4    Terminal    414    439    -10.69    +    2    KOG0003.2_1
KOG0003.2    geneid_v1.4    Exon    414    439    -10.69    +    .    KOG0003.2_1
# Sequence KOG0018.2 - Length = 3702 bps
KOG0018.2    geneid_v1.4    First    1    1113    751.23    +    0    KOG0018.2_1
KOG0018.2    geneid_v1.4    Exon    1    1113    751.23    +    .    KOG0018.2_1
KOG0018.2    geneid_v1.4    Internal    1138    2355    822.39    +    0    KOG0018.2_1
KOG0018.2    geneid_v1.4    Exon    1138    2355    822.39    +    .    KOG0018.2_1

==> local.geneid.selected.dna <==
>KOG0003.2
ATGCAGATCTTCGTGAAGACCCTGACGGGGAAAACCATCACCCTCGAGGTTGAACCATCT
GATACTATAGAAAATGTCAAAGCCAAGATCCAGGATAAAGAAGGAATTCCTCCTGACCAA
CAGCGTTTAATTTTTGCTGGCAAACAATTGGAGGATGGACGTACTCTGTCTGATTATAAC
ATCCAAAAGGAGTCAACACTTCATCTTGTGCTGAGACTTCGTGGTGGGGCAAAGAAAAGA
AAGAAGAAGTCTTACACCACTCCTAAGAAGAACAAACATAAGAGAAAGAAGGTGAAACTT
GCTGTCCTTAAATACTACAAGGTGGATGAGAATGGCAAAATCAGCCGTTTACGCCGTGAG
TGTCCCTCTGAAGAATGTGGGGCTGGAGTGTTCATGGCTAGCCACTTTGATAGACATTAT
TGTGGCAAGTGCTGCTTGACATACTGCTTCAACAAACCAGAAGACAAGTGA
>KOG0018.2

Any advice would be appreciated! Thanks!!

Suzanne