Error in rule run spades metatranscriptomes

[Sofie8]

Hi Silas,

I am trying to run atlas genecatalog on metatranscriptomics data from soil.

I ran into an error in the rule run_spades. It is looking for contig and scaffolds.fasta, but the output of spades RNA-mode is named transcripts.fasta..

rule run_spades:
    input: S16/assembly/reads/QC.errorcorr.merged_R1.fastq.gz, S16/assembly/reads/QC.errorcorr.merged_R2.fastq.gz, S16/assembly/reads/QC.errorcorr.merged_me.fastq.gz
    output: S16/assembly/contigs.fasta, S16/assembly/scaffolds.fasta
    log: S16/logs/assembly/spades.log
    jobid: 54
    benchmark: logs/benchmarks/assembly/spades/S16.txt
    wildcards: sample=S16
    threads: 8
    resources: tmpdir=/tmp/e6e7418e14, mem=170, time=48, mem_mb=170000, time_min=2880

Activating conda environment: ../atlas/conda_envs/2cca4e8d0d1322a33b3dca276ded2b0f
Waiting at most 5 seconds for missing files.
MissingOutputException in line 420 of /vsc-hard-mounts/leuven-data/314/vsc31426/mambaforge/envs/mga29/lib/python3.8/site-packages/atlas/workflow/rules/assemble.smk:
Job Missing files after 5 seconds. This might be due to filesystem latency. If that is the case, consider to increase the wait time with --latency-wait:
S16/assembly/contigs.fasta
S16/assembly/scaffolds.fasta completed successfully, but some output files are missing. 54
Job failed, going on with independent jobs.
Select jobs to execute...

These files it generates:

 * Assembled transcripts are in /lustre1/project/stg_00091/metatrans_out/S15/assembly/transcripts.fasta
 * Paths in the assembly graph corresponding to the transcripts are in /lustre1/project/stg_00091/metatrans_out/S15/assembly/transcripts.paths
 * Hard filtered transcripts are in /lustre1/project/stg_00091/metatrans_out/S15/assembly/hard_filtered_transcripts.fasta
 * Soft filtered transcripts are in /lustre1/project/stg_00091/metatrans_out/S15/assembly/soft_filtered_transcripts.fasta
 * Assembly graph is in /lustre1/project/stg_00091/metatrans_out/S15/assembly/assembly_graph.fastg
 * Assembly graph in GFA format is in /lustre1/project/stg_00091/metatrans_out/S15/assembly/assembly_graph_with_scaffolds.gfa

======= SPAdes pipeline finished.

SPAdes log can be found here: /lustre1/project/stg_00091/metatrans_out/S15/assembly/spades.log

Thank you for using SPAdes!

Where should I add the --latency-wait flag, in the assemble.smk file and then in rule run_spades?
But I think there is just no contig.fasta, instead transcript.fasta? So should the check for missing files be changed? Or in the config.yaml, genecatalog asks for 'contigs' or 'genomes'
genecatalog:
source: contigs

I attach the log files:
atlas_output2.log
spades.log

** Atlas version** 2.9.1
You briefly mentioned it also in #143

Thanks!

[SilasK]

Hey Sofie, happy new year.
Unfortunately, Atlas does not support transcriptome assembly. I never had transcriptome data to work with until recently.

One option is if you have transcriptome and metagenome data from the same samples or at least from the same sample type. Is to create MAGs and or genecatlog from the metagenome and then simply map the transcriptome reads to the MAGs / Genecatalog.

There is an option to set genome_dir in a config file to point to a directory with dereplicated MAGs. and then you can run atlas with atlas run quantify_genomes

[Sofie8]

Hi Silas,

Happy new year to you too!

Yes, I agree with that approach, which was used I think by @jmtsuji, so both gathering DNA metagenomes and transcriptomes from the same sample.

In this case, the phd students already brought me the metatranscriptome data and no DNA shotgun metagenomes. But the experimental setup and research question is the following:

• pot experiment with/without inoculation with a consortium of bacteria (inoculated once, or 2 times) and monitoring the effect on oil dissipation. So linking oil degradation results with metatranscriptome data, so fold changes compared to the non-inoculated control etc.

So, indeed, my idea was the suggestion you give below, no mags, but gene prediction. Will reply there.

[SilasK]

On the other hand it should not be so difficult to implement the assembly of transcriptomes. As you mention the spades rule is already adapted.

One option is to create symbolic links from {sample}/assembly/transcripts.fasta to {sample}/assembly/{sample}_final_contigs.fasta
then you could continue to run atlas.

What do you want to do with the assembled transcripts? e.g. Genecatalog?

Questions I have:

• If it would make sense to use some of the filtered transcripts
• Can one simply run prodigal on assembled transcripts?

[Sofie8]

Ok, I renamed the transcripts.fasta file to contigs.fasta and scaffolds.fasta (because spades keeps on searching both).. there is no scaffold_transcripts.fasta..

Now it is running, see the log so far:
atlas_output.log

Questions I have:

• If it would make sense to use some of the filtered transcripts
• Can one simply run prodigal on assembled transcripts?
  What is atlas using now, prodigal? See here:
  "Along with GeneMarkS-T we assessed the performance of specialized tools for gene prediction in transcripts: ESTscan, v2.1 (4), TransDecoder (http://transdecoder.sourceforge.net/), as well as prokaryotic Prodigal, v2.60 (8) used in the special mode of predicting ‘intronless genes’. Source: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4499116/

We are interested in prokaryotes for now, some other pipelines that we will give a try:

• HUMAnN 3.0

For eukaryotes, maybe these:

• FINDER
• META-trans

Happy to discuss this further if other atlas users are also running metatranscriptomes and have ideas?

[SilasK]

Is atlas continue with the renamed files?

I think the transcripts are the cds of the genes. They are the genes. So we could actually only try to translate them and do the gene catalog.

[Sofie8]

Hi Silas,
Yes it worked by renaming the output files!

Now, looking into the genecatalog output.. is there somewhere the merged abundance table also? (merged in the sense, counts.tsv with eggnog.tsv)?
And in the orf2gene file, I have a question about the interpretation, so cytP450 (Gene08141) is found in all samples (good), and then in the orf2gene file when I search for the Gene08141, it says only S46_2410_1, so in sample 46...

And another interpretation question: the eggnog.tsv; gives also a column with a taxonomical information, like for CYP450, 'Actinobacteria', this is an estimated taxonomical placement (it could maybe be any bacterium?).

Excel with the output: metatrans_out.xlsx

• first tab: eggnogg.tsv
• 2nd: counts
• 3rd: orf2gene

Can I point atlas_analyze to the genecatalog output folder?

Cheers!
Sofie

[SilasK]

I see I need some more documentation about that,

As background: atlas gene catalog takes all the orfs predicted from all samples and clusters them at 90% id on the amino acid level. The representative gene is then annotated with eggNOG mapper.
reads are mapped to the cds of the representative to get the counts.

This means the eggNOG annotation e.g. taxonomy is one of the representative sequences, which might not fully transferable to all the genes in the cluster. If you want you can cluster at higher identity level. But, there are 44 annotated cytP450 genes in your data from different phyla.

How to analyze the data?

Bravo you managed to open the files. They were in parquet isn't it?
You don't need the orf2gene table. You can combine the other two to get what you want.
There is the gene name in the first column of the eggnog annotations, which you can combine with the gene counts table.

What I would do:

• Check mapping rate of the gene catalog (sum of read counts divided by total reads in
• Load median coverages and calculate gene copies per million.
• Select all genes for annotation of interest and sum the gene abundance of this annotation.

You can do the last step as matrix multiplication.
@Sofie8 Do you prefer R or python?

You got quite some annotation from the plant (host) isn't it?
Do you think dram annotation on the genecatalog that gives you KEgg/CAZy would be more informative?

[SilasK]

See my PR #615