Quantification of genes present in the gene catalog among samples

[fconstancias]

Dear atlas developers/users,

I have generated a gene catalog using atlas. Do you have any guidline how to map back and quantify the proportion of genes (I guess noramlizing by gene length) in the different samples used to generate the gene catalog?

Thanks a lot.

[SilasK]

My personal view is that one should quantify genomes and not genes. The whole point of assembly-based metagenomics is to work with the genomes instead of genes.

What do you think @fconstancias : Why would you quantify genes instead of genomes?

I made some scripts on how to use the genome output of atlas.

But obviously you can do everything with atlas:
The following command:

atlas run genecatalog combine_gene_coverages <other parameters> 

gives you the two tables:

• "Genecatalog/counts/Nmapped_reads.tsv.gz"
• "Genecatalog/counts/median_coverage.tsv.gz"

Where the second one is the gene counts normalized to the gene length.

However, I would take the "Nmapped reads" and use a CODA approach. For more information:

Microbiome Datasets Are Compositional: And This Is Not Optional https://doi.org/10.3389/fmicb.2017.02224
https://github.com/ggloor/CoDa_microbiome_tutorial

The tutorial is for 16S sequencing but the same logic applies also for genome- or gene-based metagenomics.

Last thing: In the next days I will work on an extension for the gene catalog workflow on https://github.com/metagenome-atlas/genecatalog_atlas

E.g. Make genecatalogs at 100id, 90id and 50id.
Any other ideas?

[fconstancias]

@SilasK thanks for your feedback.

My personal view is that one should quantify genomes and not genes. The whole point of assembly-based metagenomics is to work with the genomes instead of genes.

Well, it all depend on your question afterall but it is good to have both because most of the case reconstructing all the MAGs from all the populations of all our sample is impossible. And if you want to compare oral samples from healthy patients or from patient with dental caries for instance, you migh want to check if genes present and/or expressed are specific. I agree it is even better to do it at the MAG level.
In the current dataset I am working on, we are interested in genes involved in N cycle in the context of wastewater. We are getting MAGs as well which are very interesting but we would like to compare proportion of genes involved in N cycle is different among our disturbance we induce on our system.

But obviously you can do everything with atlas:
Excellent !

The following command:

atlas run genecatalog combine_gene_coverages

I have already sucessfully generated a genecatalog with atlas genecatalog .... but without the combine_gene_coverages option. Can I start from there or do I have to redo the entire process?

Last thing: In the next days I will work on an extension for the gene catalog workflow on https://github.com/metagenome-atlas/genecatalog_atlas
E.g. Make genecatalogs at 100id, 90id and 50id.
Yes, that would be great. It might be great as well to be able to choose read mapping software (e.g., Bowtie2, bwa).

thanks

[SilasK]

I have already sucessfully generated a genecatalog with atlas genecatalog .... but without the combine_gene_coverages option. Can I start from there or do I have to redo the entire process?

Yes, of course. Just run atlas in the same working directory.

Yes, that would be great. It might be great as well to be able to choose read mapping software (e.g., Bowtie2, bwa).

Do you have experience with bwa. I've heard you have to tweak the parameters to get a mapping done in a reasonable time.

About the mapping, now it's done with bbmap, and I map the paired-end reads. After reflection, this is maybe not the best way to do it. As the two ends of paired-end reads might mao to different genes.

[fconstancias]

Do you have experience with bwa. I've heard you have to tweak the parameters to get a mapping done in a reasonable time.

I usually use Bowtie2 but I haven't benchmark the different options.

Can I start from there or do I have to redo the entire process?
Yes, of course. Just run atlas in the same working directory.

I did but it seems atlas did took the already processed assembly but somehow is running again the annotation :

Removing temporary output file Genecatalog/subsets/genes/subset3.emapper.seed_orthologs.
[Fri Jan 24 16:13:57 2020]
Finished job 548.
14 of 184 steps (8%) done

[SilasK]

Sorry, this shouldn't happen.
You can always run atlas with --dryrun first to see what's he will do.

You have the combined annotation ? eggnog.tsv

Then stop the running process.
And you can run

atlas run None combine_gene_coverages <other parameters> 

This shouldn't start the annotation.

[soilmicrobiome]

I am having some issues running "genecatalog combine_gene_coverages". Atlas can't seem to find a sam file. I can see while its running bam and sam files are created, but then the sam is removed and only the bam remains. However atlas seems to want to find the missing sam:

Currently running this via singularity (atlas, version 2.6a2) on an HPC

[2021-07-29 10:46 INFO] Executing: snakemake --snakefile /usr/local/lib/python3.7/site-packages/atlas/Snakefile --directory /scratch/user/metagenome-atlas --rerun-incomplete --configfile '/scratch/user/metagenome-atlas/config.yaml' --nolock --use-conda --conda-prefix /scratch/user/metagenome-atlas/databases/conda_envs --scheduler greedy genecatalog combine_gene_coverages --cores 24
localrules directive specifies rules that are not present in the Snakefile:
verify_eggNOG_files

Building DAG of jobs...
Updating job combine_egg_nogg_annotations.
Using shell: /usr/bin/bash
Provided cores: 24
Rules claiming more threads will be scaled down.
Singularity containers: ignored
Job stats:
job count min threads max threads

---

align_reads_to_Genecatalog 3 24 24
combine_gene_coverages 1 1 1
convert_sam_to_bam 3 4 4
pileup_Genecatalog 3 24 24
total 10 1 24

[Thu Jul 29 10:46:44 2021]
rule align_reads_to_Genecatalog:
input: HKS2/sequence_quality_control/HKS2_QC_R1.fastq.gz, HKS2/sequence_quality_control/HKS2_QC_R2.fastq.gz, HKS2/sequence_quality_control/HKS2_QC_se.fastq.gz, Genecatalog/gene_catalog.fna
output: Genecatalog/alignments/HKS2.sam
log: logs/Genecatalog/alignment/HKS2_map.log
jobid: 71
wildcards: sample=HKS2
threads: 24
resources: tmpdir=/tmp, mem=250, java_mem=212, time=72

ESC[33mlocalrules directive specifies rules that are not present in the Snakefile:
verify_eggNOG_files
ESC[0m
[Thu Jul 29 11:28:05 2021]
Finished job 71.
1 of 10 steps (10%) done

[Thu Jul 29 11:28:05 2021]
rule convert_sam_to_bam:
input: Genecatalog/alignments/HKS2.sam
output: Genecatalog/alignments/HKS2.bam, Genecatalog/alignments/HKS2.forpileup.sam
jobid: 72
wildcards: file=Genecatalog/alignments/HKS2
threads: 4
resources: tmpdir=/tmp, mem=8, time=72

ESC[33mlocalrules directive specifies rules that are not present in the Snakefile:
verify_eggNOG_files
ESC[0m
[bam_sort_core] merging from 60 files...
Removing temporary output file Genecatalog/alignments/HKS2.forpileup.sam.
[Thu Jul 29 12:04:55 2021]
Finished job 72.
2 of 10 steps (20%) done
Traceback (most recent call last):
File "/usr/local/lib/python3.7/site-packages/snakemake/init.py", line 780, in snakemake
keepmetadata=keep_metadata,
File "/usr/local/lib/python3.7/site-packages/snakemake/workflow.py", line 1054, in execute
success = scheduler.schedule()
File "/usr/local/lib/python3.7/site-packages/snakemake/scheduler.py", line 477, in schedule
run = self.job_selector(needrun)
File "/usr/local/lib/python3.7/site-packages/snakemake/scheduler.py", line 809, in job_selector_greedy
c = list(map(self.job_reward, jobs)) # job rewards
File "/usr/local/lib/python3.7/site-packages/snakemake/scheduler.py", line 891, in job_reward
temp_size = self.dag.temp_size(job)
File "/usr/local/lib/python3.7/site-packages/snakemake/dag.py", line 597, in temp_size
return sum(f.size for f in self.temp_input(job))
File "/usr/local/lib/python3.7/site-packages/snakemake/dag.py", line 597, in
return sum(f.size for f in self.temp_input(job))
File "/usr/local/lib/python3.7/site-packages/snakemake/io.py", line 251, in wrapper
return func(self, *args, **kwargs)
File "/usr/local/lib/python3.7/site-packages/snakemake/io.py", line 569, in size
return self.size_local
File "/usr/local/lib/python3.7/site-packages/snakemake/io.py", line 574, in size_local
self.check_broken_symlink()
File "/usr/local/lib/python3.7/site-packages/snakemake/io.py", line 579, in check_broken_symlink
if not self.exists_local and os.lstat(self.file):
FileNotFoundError: [Errno 2] No such file or directory: 'Genecatalog/alignments/HKS2.sam'
[2021-07-29 12:04 CRITICAL] Command 'snakemake --snakefile /usr/local/lib/python3.7/site-packages/atlas/Snakefile --directory /scratch/user/metagenome-atlas --rerun-incomplete --configfile '/scratch/user/metagenome-atlas/config.yaml' --nolock --use-conda --conda-prefix /scratch/user/metagenome-atlas/databases/conda_envs --scheduler greedy genecatalog combine_gene_coverages --cores 24 ' returned non-zero exit status 1.
(END)

thoughts?

[SilasK]

@soilmicrobiome Thank you for the reporting. I solved the issue in v. 2.6a3 (master branch). I'm testing it and then releasing it via conda soon.

[soilmicrobiome]

Great thanks very much! Will this be available in the singularity/docker as well?

[SilasK]

Yes, yes. As you do in #359 You install atlas outside the singularity.

[soilmicrobiome]

Until a new conda comes with the gene abundance part arrives, wondering about running this externally - to generate e.g FPKM (fragments per kilobase per million mappable reads) or TPM (transcripts per million).

I assume the correct files to use would be:
qc-ed reads to map
/path_to_samples/sample1/sequence_quality_control/sample1_QC_R1.fastq.gz
/path_to_samples/sample1/sequence_quality_control/sample1_QC_R2.fastq.gz

/path_to_samples/sample2/sequence_quality_control/sample2_QC_R1.fastq.gz
/path_to_samples/sample2/sequence_quality_control/sample2_QC_R2.fastq.gz

genes to map to
/path_to_samples/Genecatalog/gene_catalog.fna

But if i have more than one sample - the genecatalog is combined for all samples? Would this lead to issues? or just map each sample one at a time against it?

[SilasK]

New version should be on conda ! 👍🏻

[soilmicrobiome]

Fantastic will try it out on our HPC! thanks

[jmtsuji]

@soilmicrobiome I haven't tried this ATLAS rule in quite some time, although I might be interested in using it for a new set of soil metagenomes our lab received.

[soilmicrobiome]

just re-testing on our hpc - might be an issue with that rather than atlas (i don't have much control with the HPC!).

[soilmicrobiome]

got the version running (was issue with mamba and our HPC again). I have tested it with the sample reads and everything thing runs including the gene abundance part. Thanks!

[SilasK]

This calculation is implemented in atlas v2.8

[slambrechts]

Hi everyone,

Great new feature! Before atlas I searched the raw or clean reads using diamond whenever I needed to quantify some specific proteins for which I had the sequences in fasta files among samples. Usually, these genes are not yet in the KEGG database etc. How would you look for these genes using the atlas Gene Catalog? Would you blast them against the atlas_directory/Genecatalog/gene_catalog.faa file to find them and then check the results of the new align_reads_to_Genecatalog rule?

[SilasK]

Yes, I would do that.

search against the gene catalog. (be aware of the clustering id)
And then combine with the quantification output. You can even check if your gene is part of a MAG (#413 ).

Have also a look at #480