Crash on large metatranscriptomics database

[fdelogu]

Hi,

I am trying to run mmseq-linux v1.0.9 using a large metatranscriptomics dataset as reference (>4 million ORFs). The run stops after saying “Amalgamating transcripts and calculating summary statistics...”, due to memory limits (SIGSEGV: 11) when I used the maximum RAM available to me (3TB).
Before the crash the following output files are produced along the way:
sample.k,
sample.M,
sample.gene.trace_gibbs.gz,
sample.identical.trace_gibbs.gz,
sample.prop.trace_gibbs.gz,
sample.trace_gibbs.gz.

I’ve run some tests on a E.coli dataset and observed that the memory consumption jumps during the “Amalgamation of transcripts and calculation of summary statistics” requiring several times the one used during the rest of the computation. Is there a way to circumvent this problem and complete the run?

Best regards,
Francesco

[eturro]

Hi Francesco, If you ensure that each "transcript" is assigned to a different "gene" in the and each transcript has a different sequence, then perhaps the memory requirements will be reduced. Instead of re-running bam2hits with amended FASTA files, you could try modifying the header of the .hits files using hitstools. Could you try that? bw Ernest

…

On 10 Feb 2021, at 10:51, Francesco Delogu ***@***.***> wrote: Hi, I am trying to run mmseq-linux v1.0.9 using a large metatranscriptomics dataset as reference (>4 million ORFs). The run stops after saying “Amalgamating transcripts and calculating summary statistics...”, due to memory limits (SIGSEGV: 11) when I used the maximum RAM available to me (3TB). Before the crash the following output files are produced along the way: sample.k, sample.M, sample.gene.trace_gibbs.gz, sample.identical.trace_gibbs.gz, sample.prop.trace_gibbs.gz, sample.trace_gibbs.gz. I’ve run some tests on a E.coli dataset and observed that the memory consumption jumps during the “Amalgamation of transcripts and calculation of summary statistics” requiring several times the one used during the rest of the computation. Is there a way to circumvent this problem and complete the run? Best regards, Francesco — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[fdelogu]

Hi Ernest,

thank you for your answer! In which way should I change the header with hitstools?

In the meantime I used hitstools to inspect the .hits file of one of my samples with the following results:

Hits file: mmseq.hits
Total number of transcripts: 4384648
@TranscriptMetaData FHBPBLAL_00006 45.8755 195
@TranscriptMetaData FHBPBLAL_00007 335.415 738
@TranscriptMetaData FHBPBLAL_00008 441.151 930
...
Total number of genes: 4384648
@GeneIsoforms AALFAMOJ_00005 AALFAMOJ_00005
@GeneIsoforms AALFAMOJ_00006 AALFAMOJ_00006
@GeneIsoforms AALFAMOJ_00007 AALFAMOJ_00007
...
Total number of identical transcript sets: 203329
@IdenticalTranscripts JCLHNBGN_131468 BCBBMDAJ_13859
@IdenticalTranscripts BOELJLNI_143620 MPJPILBE_63237
@IdenticalTranscripts FHBPBLAL_47327 FGMGOPFI_79015
...
Reads:
>HWI-ST201:374:BD27D5ACXX:3:1101:1279:9806
BHPLJFKF_296870
>HWI-ST201:374:BD27D5ACXX:3:1101:1280:16969
CJCGNBID_09091
>HWI-ST201:374:BD27D5ACXX:3:1101:1280:28828
EBBDFONC_09101
MCNLMFIM_217103
...
Total number of reads: 10185956

If I interpret this correctly, each transcript was assigned to a gene (which was also the initial intent).

Moreover, parsing the results of hitstools header on the same .hits file I calculated that the sequences grouped in the 203329 identical transcripts sets are 567563. The reduction of the transcript database size would be of 364234 sequences (~8% of the total), probably not enough to impact largely on the memory usage?

Best regards,
Francesco

[eturro]

Dear Francesco, So there is a one to one mapping between transcripts and genes but you do have quite a lot of duplicated transcript sequences in your FASTA. You don't lose anything by removing those beforehand from your FASTA. Could you please try that? So, for example, if you see @IdenticalTranscripts A B C Pick a transcript to retain and remove the other two. E.g. remove B and C and retain A (perhaps renaming it A-B-C). I appreciate that doesn't seem like a tremendous reduction but by doing this you should in principle sidestep the amalgamation step entirely (where all these big MCMC traces get added up). If that's the specific step that is causing issues then perhaps this will resolve the problem. Best wishes Ernest

…

On 11 Feb 2021, at 08:47, Francesco Delogu ***@***.***> wrote:  Hi Ernest, thank you for your answer! In which way should I change the header with hitstools? In the meantime I used hitsools to inspect the .hits file of one of my samples with the following results: Hits file: mmseq.hits Total number of transcripts: 4384648 @TranscriptMetaData FHBPBLAL_00006 45.8755 195 @TranscriptMetaData FHBPBLAL_00007 335.415 738 @TranscriptMetaData FHBPBLAL_00008 441.151 930 ... Total number of genes: 4384648 @GeneIsoforms AALFAMOJ_00005 AALFAMOJ_00005 @GeneIsoforms AALFAMOJ_00006 AALFAMOJ_00006 @GeneIsoforms AALFAMOJ_00007 AALFAMOJ_00007 ... Total number of identical transcript sets: 203329 @IdenticalTranscripts JCLHNBGN_131468 BCBBMDAJ_13859 @IdenticalTranscripts BOELJLNI_143620 MPJPILBE_63237 @IdenticalTranscripts FHBPBLAL_47327 FGMGOPFI_79015 ... Reads: >HWI-ST201:374:BD27D5ACXX:3:1101:1279:9806 BHPLJFKF_296870 >HWI-ST201:374:BD27D5ACXX:3:1101:1280:16969 CJCGNBID_09091 >HWI-ST201:374:BD27D5ACXX:3:1101:1280:28828 EBBDFONC_09101 MCNLMFIM_217103 ... Total number of reads: 10185956 If I interpret this correctly, each transcript was assigned to a gene (which was also the initial intent). Moreover, parsing the results of hitstools header on the same .hits file I calculated that the sequences grouped in the 203329 identical transcripts sets are 567563. The reduction of the transcript database size would be of 364234 sequences (~8%) of the total, probably not enough to impact largely on the memory usage? Best regards, Francesco — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[fdelogu]

Dear Ernest,

I modified the .hits file as suggested, but the run halted in the same step with the same error.

Best regards,
Francesco

[eturro]

Do you mind making one of your hits files and the corresponding FASTA file available to me so I can try it on my server?

…

On 18 Feb 2021, at 03:42, Francesco Delogu ***@***.***> wrote: Dear Ernest, I modified the .hits file as suggested, but the run halted in the same step with the same error. Best regards, Francesco — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[fdelogu]

Dear Ernest,

sure! Please find the .hits file and the reduced version as well as the FASTA file of the transcripts at the following link: https://www.dropbox.com/sh/uhncijsbt9n02fw/AAAsL9m8Txg6k1sedPE53mm1a?dl=0

Best regards,
Francesco

[eturro]

Dear Francesco, Your analysis includes more than 4M transcripts. This is a very large number of transcripts and is likely causing you to run out of memory. As a workaround, you can try making the following change to mmseq.cpp and recompiling: - int trace_length=1024; + int trace_length=256; Hopefully this will solve the issue, but note you will still require around 30GB of memory. Best wishes, Ernest

…

On 22 Feb 2021, at 06:02, Francesco Delogu ***@***.***> wrote: Dear Ernest, sure! Please find the .hits file and the reduced version as well as the FASTA file of the transcripts at the following link: https://www.dropbox.com/sh/uhncijsbt9n02fw/AAAsL9m8Txg6k1sedPE53mm1a?dl=0 Best regards, Francesco — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.