Regarding large dataset

[chanchalrana]

I have around 80k E. coli genomes for which I want to do pangenome analysis. Is it possible to do so? If yes, what will be the difficulties I can get to do the analyses? I have read that the largest dataset you have analyzed so far is approx. 20k.

[axbazin]

Hi,

The largest that we published was indeed around 20K, but I did manage a 100K S. enterica pangenome run once, just to see if it was possible.

It takes days to run, and you'll need a lot of RAM (it used over 1.5 To if I remember correctly), but it is possible.

A few tips on the technical aspects:

• the biggest problem in terms of RAM comes from the "partition" step. As of today, it is possible to customize the number of CPUs used for each step using a config file (explained here how to do it: https://ppanggolin.readthedocs.io/en/latest/user/practicalInformation.html#configuration-file). The higher the number of CPUs the higher the RAM consumption at the "partition" and "rarefaction" steps. It is not true for other steps, so you can use "a lot" of CPUs for the other steps (annotation clustering and so on), and something like 4 or 8 for the partition (and rarefaction if you need it) step.
• The HDF5 file resulting will be quite big, not as big as the size of your input data but still a few hundred gigas. We've vastly improved the data structure since my 100K test so maybe it's not so big today. The other output files (gexf and so on) will be very big also, so I'd recommend running ppanggolin with a "step by step", "only get the HDF5 file first" route, and then only get the other output files that you want and need for your analysis. Every possible output file and related options, and means to get them using the h5 file, should be described in this section: https://ppanggolin.readthedocs.io/en/latest/user/PangenomeAnalyses/pangenomeAnalyses.html#pangenome-outputs

In terms of the additional analysis tools that ppanggolin provides (e.g. aligning genes on the pangenome, doing a MSA, and so on), I don't quite know how well they'll cope with such a big pangenome, but I expect it should be fine as long as you manage the "main" pangenome computations.

I think @JeanMainguy did a few "big" pangenomes as well, so maybe you have some tips too?

Good luck if you try out ppanggolin for your analysis in the end, I hope it will go smoothly !

Adelme

[JeanMainguy]

Hi,

I ran ppanggolin all on 32K E. coli genomes and it finished in about 6 hours, using 727GB of memory. The final pangenome HDF5 file is 12GB.

I used 24 CPUs, except for the partition step, where I set it to 12. Partitioning is definitely the bottleneck here, so you could use even fewer CPUs for this step.

Here is my config_file to adjust the CPU usage in the different steps when using the ppanggolin all command:

annotate:
    cpu: 24

cluster:
    cpu: 24

partition:
    cpu: 12 # bottleneck step - could try lowering it further

rarefaction:
    cpu: 24

module:
    cpu: 24

write_pangenome:
    cpu: 24

write_genomes:
    cpu: 24

Otherwise, if you want to run ppanggolin step by step, as Adelme suggests, you can set the --cpu option for each of them on the command line, so you don't need a config file.

In ppanggolin all you can also choose to skip the creation of output files other than the HDF5 file with the --no_flat_files option. You can then later use the HDF5 file to produce the output you want with the ppanggolin write_genomes and ppanggolin write_pangenome commands.

Let us know how it goes.

Best,
Jean

[chanchalrana]

Thanks a lot for the help. I would definitely follow the instructions for smooth analysis.

[chanchalrana]

Hello again!!

Because the computational resources I have didn't allow me to apply ppanggolin on all 80k E. coli genomes. So, I have decided to do it in chunks and then later on compare the gene families in all the chunks as you have described in issue #68 Merge pangenome graphs.

So, I am particularly interested in studying the core genome (present in 100% of the genomes). I am using the command:
ppanggolin fasta -p pangenome.h5 --output core_genome --genes core

My question here is how the core genome is identified in this pipeline, does it have a threshold of 100% identity? What about the msa of core gnome alignment we get using the command: ppanggolin msa -p pangenome.h5 --source dna --partition core --output core_genome_aligned? How does this command works?

Also, what are the default identity percentage used for persistent and core genomes? Can this be change?

The msa core dna we get are the unique CDS regions and are not being repeated. Right? For example in the 10C_CDS_0494.aln file I have got an aligned chunk of 1500 bp so this region will not be repeat in any other CDS.aln files that we will get. Right?

[axbazin]

Hello,

All of the "core" is 100% presence in the pangenome. It is indeed a pretty bad concept when you have too many genomes as it tend to go to 0. If you with to use a threshold of presence, you can look at the "softcore" partition together with its companion option --soft_core with which you can set the threshold you wish between 0 and 1.

For the record, ppanggolin integrates another concept, which is at the heart of the tool actually, and is called the 'persistent' genome. You can read more about how it is done in the ppanggolin article, but roughly, it is a statistic-based partition which relies on the patterns of presence/absence of gene families and their organisation in the genomes, and that is much more resilient to sampling than the 'core' or 'softcore' concepts. It ends up being extremely useful when using low-quality genomes, or many genomes, or genomes whose species assignation is unsure (or basically any kind of "unperfect" genome collection).

The MSA command uses mafft basically, and the partitioning results of whichever partition you chose.

For the identity percentage, it's the same for all partitions (as it is done before partitioning) and is set during the 'cluster' step. The default is 0.80, but can of course be changed as you wish during the 'cluster' step.

For your last question, yes you are correct. With msa core dna, in a given .aln file you will get a CDS sequence that should not be seen in other .aln file.

Adelme

[chanchalrana]

Yup!! You are right for the core genome part. I just got the results after analyzing 5k sequences for the core genome and the size of each chunk is too small. I will shift for 'persistent' genome.

Thanks a lot!

[chanchalrana]

Hello again! I have a few questions to ask. Out of 80K genomes, I am currently doing the analysis on 5k genomes:

1. I want to find out the largest chunk of genes present in all the genomes. For this, I am using python code by applying the sliding window approach. For this I give the .tsv files in the table folder as input. However, this is a lengthy process. I recently discovered about your approach to classify the persistent genes based on the neighboring gene arrangement (the pangenome graph in which the thickness of the edges represent the genomes in which that gene arrangement is present, if I understood this correctly!). Is there any way I can do this kind of analysis using your pipeline?

2. As discussed earlier, I do not have the computational power to run all the 80k genomes in one go. However, I want to do the analysis on all the genomes. Is there any way to combine the gene families, RGP, Spots, modules of all the pangenome batches (If I run it in a batch of 5k)? I don't know if Merge pangenome graphs #68 will work for me. What about the projection method you have given in the docs? Or any other method?

[axbazin]

Hello,

1. To look for the largest chunks of genes conserved across genomes, I believe using the pangenome graph directly, filtering out edges based on their weight / thickness (the number of genomes that have the edge), and looking for the largest connected component is probably the way I would have done it. Using that approach, each connected component would be a "chunk of genes shared across all genomes".
   You need to be a bit familiar with graph theory concepts, and a big enough computer to load the graph object that you got, but should be feasible. The graph object is stored both in the .json file and .gexf files, so you can use any library that copes with those formats. Trying that approach using the .gexf with networkx (python package for graphs) would be my initial attempt.

2. That problem is far from trivial, especially to avoid reference bias and to guarantee some sort of reproducibility. Most likely you can perform the clustering step of your 80K genomes together, it should not take so much resources. Then, ideally, performing the other steps on chunks of 5k genomes using the same clustering as input. You would then have the same basis for all of your pangenomes, as gene family identifiers would be shared across all of them.
   Getting RGP would then be trivial. Gene families and their classification as well. Merging the spots and modules is not so easy however ... For sure you can compare them, since they will have the same gene family identifiers, but having one set of "spots" and "modules" describing all 80K genomes together is probably not as simple. The way I think of it, it would involve a lot of re-coding the methods to not load all the pangenome graphs at once, which is maybe feasible but would involve quite a bit of dev work.

I don't quite know if all of this helps, but good luck !

Adelme

[chanchalrana]

Thanks a lot for the quick response!

Further questions. Sorry to not have asked in the previous query.

1. Can you explain how modules are made through your pipeline? I did not understand the exact mechanism.

2. The Spots are the region where the RGP genes are inserted, and these are surrounded by 6 persistent genes (3 on each side). Is there any way to find out (using your pipeline) that in how many genomes these 6 genes are found consistent? For example: in this: spot nb_rgp nb_families nb_unique_family_sets mean_nb_genes stdev_nb_genes max_nb_genes min_nb_genesspot_42 4730 30 4730 10.814 2.698 18 5
   We do not know that in how many genomes the 6 persistent genes were found exactly same to have named it as spot_42.

[axbazin]

No worries !

1. The module algorithm and example application have been described thoroughly in the following preprint. To sum it up roughly and briefly, it relies on edge weights between gene families that are not persistent. There is a transitive closure of the graph to cope with gaps and annotation/assembly errors. then, the algorithm goes over all edges (the original ones and those from the transitive closure), filters out those that do not have a high weight, and searches for connected components with at least 3 gene families. Those connected components are the modules.

2. Indeed. The spot detection algorithm allows some minor variability in the bordering persistent genes (to cope with annotation/assembly errors, still getting the RGP when a single persistent gene family is sometimes missing), but roughly yes we can define a spot by being a bunch of RGP surrounded by the same persistent genes. From the default output it's not so easy to get what you want, however I see two ways of getting it in another way:

Option 1: The ppanggolin draw -p pangenome.h5 --draw_spots option, which actually does more than just drawing. It also writes a .tsv file listing all RGP, and their genomes, where the spot is found. If you take the last column, you have the list of genomes that have the spot, basically. However be wary that, in the case where some genomes do not have a RGP but do have the 6 persistent genes, it will not be listed there.

Option 2: Using 'gene table' files and some greps. You can get the list of bordering persistent families from the 'spot_borders.tsv' file (I don't remember if it is written by default or not, but if not you can refer to this option to get it.
Once you have the gene families identifiers you are looking for, you can use some grep commands to find your gene families in the files written by the ppanggolin write_genomes --table option, that you already talked about earlier. Then you can probably filter out files that do not have all gene families identifiers in them, and you get the number of genomes that have the set of gene families you are looking for this way.
You can actually also use those files to get all of the genomes that have 'spot_42', which should give the same result as 'option 1' for the number of genomes.

Hope this helps !

Adelme

[chanchalrana]

Okay, Thanks a lot.