Custom library and hmmer
Hi, I have been trying to use repeatmasker with hmmer as a search engine but so far I have not success. After produce de novo model predictions with RepeatModeler2, I used the stk file to produce a hmm with:
hmmbuild --dna --cpu 12 $hmm $stk
Then I used :
$srcPath/RepeatMasker -engine hmmer -lib $hmm $genome
I got, however, the following message: Search Engine: HMMER [ 3.1b2 (February 2015) ] Using Custom Repeat Library: /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/ce.hmm
analyzing file /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/caenorhabditis_elegans.PRJNA13758.WBPS14.genomic.fa
identifying Simple Repeats in batch 1 of 1732
identifying matches to ce.hmm sequences in batch 1 of 1732
sh: línea 1: 842750 Violación de segmento (core' generado) /home/isaacmartinezugalde/storage/software/hmmer-3.1b2/src/nhmmscan --cut_ga --cpu 2 /home/isaacmartinezugalde/storage/Data/Annotations/Hmm//ce.hmm /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/caenorhabditis_elegans.PRJNA13758.WBPS14.genomic.fa_batch-1.masked 2> /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/hmmerResults-1616183246-842745.err WARNING: The search engine returned an error (139, status = 139 ) Engine parameters: /home/isaacmartinezugalde/storage/software/hmmer-3.1b2/src/nhmmscan --cut_ga --cpu 2 /home/isaacmartinezugalde/storage/Data/Annotations/Hmm//ce.hmm /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/caenorhabditis_elegans.PRJNA13758.WBPS14.genomic.fa_batch-1.masked A search phase could not complete on this batch. The batch file will be re-run and if possible the program will resume. WARNING: Retrying batch ( 1 ) [ 255,, 60680]... identifying Simple Repeats in batch 1 of 1732 identifying matches to ce.hmm sequences in batch 1 of 1732 sh: línea 1: 842766 Violación de segmento (core' generado) /home/isaacmartinezugalde/storage/software/hmmer-3.1b2/src/nhmmscan --cut_ga --cpu 2 /home/isaacmartinezugalde/storage/Data/Annotations/Hmm//ce.hmm /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/caenorhabditis_elegans.PRJNA13758.WBPS14.genomic.fa_batch-1.masked 2> /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/hmmerResults-1616183247-842760.err
WARNING: The search engine returned an error (139, status = 139 )
Engine parameters: /home/isaacmartinezugalde/storage/software/hmmer-3.1b2/src/nhmmscan --cut_ga --cpu 2 /home/isaacmartinezugalde/storage/Data/Annotations/Hmm//ce.hmm /home/isaacmartinezugalde/storage/Data/Annotations/Hmm/RM_842716.FriMar191347152021/caenorhabditis_elegans.PRJNA13758.WBPS14.genomic.fa_batch-1.masked
A search phase could not complete on this batch.`
I can notice that this could be probably related to the --cut_ga parameter within nhmmer, due to my hmm file was generated from raw predictions I do not have GA lines in my hmm file. Is it possible to give an e-value as a cut-off to nhmmer instead a --cut_ga within RepatMasker?
Thanks for your support Regards,
Thanks for your feedback. I am glad to see someone building HMMs from TE family alignments. The issue in this case ( I also suspect ) is that you don't have a HMM with a threshold set. The '--cut_ga' parameter being passed to nhmmscan by RepeatMasker is indicating that the gathering threshold should be used and accordingly nhmmscan is looking for lines like:
GA 32.37; TC 32.37; NC 32.37;
Where 'GA' is the gathering threshold ( in bits ), 'TC' is the trusted cutoff, and 'NC' is the noise cutoff. These are values that are typically determined by comparing the score distributions of hits to a genome expected to contain instances, and the score distribution to a benchmark sequence that does not contain any instances. Given a false discovery rate that is tolerable, the GA threshold can be determined from these two distributions. For more detail please see the first Dfam paper: https://doi.org/10.1093/nar/gks1265
Unfortunately we don't have a way of supplying an E-value threshold to RepeatMasker for this type of search. It's a good idea though, and I will add it to our feature request list.
Dear Robert, Thanks for your support. I have been using a single hmm to tackle the ideas in the 2013 Dfam paper. I have, however, some doubts.
As you propose in the paper I used the reverse genome from which my model comes from as a benchmark sequence. Using the naïve genome and the reverse one, I got a coverage of 1270735 and 4969 bases respectively, so I have a FDR of 0.003925686. Must I directly use this value as a GA value? this is not clear at all for me. Also, I’m comparing my example with the DF0000001.4 hmm within Dfam_curatedonly.hmm file, in this hmm you have the estimation of GA, TC, and NC for three species, but you use the same value for GA, TC, and NC. If at the end I will try to use a model to annotate different genomes, should I do the same?
Again, thanks for your support