clustering bins from multiple samples

[jflucier]

Need to try this:

MeShClust

[jflucier]

algo:

1. cat all sample bins together
2. run meshclust all vs all
3. parse output and select only representative bins (sequence representing the center of a cluster)
4. rebuild fasta the representating sequences
5. align sample reads against this new fasta
6. produce count matrix
7. build heatmap

[jorondo1]

Faut aussi essayer dRep qui est très répandu dans la litterature.

Une commande pour faire un primary cluster à 90% puis un secondary à 99% avec un min coverage de 30% sur des MAGs de >50% completeness et <5% contamination

dRep dereplicate out_dir/ -g /path/to/bins/*.fasta -pa 0.9 -sa 0.99 -nc 0.30 -p 24 -d -comp 50 -con 5

je m'inspire de cette étude qu'on a décortiquée en détail au labo

[jflucier]

ca l'air parfait ca!

[jorondo1]

yes d'autant plus que MeshClust se veut un improvement sur CD-hit, en contournant les limitations apportées par les algo greedy. L'affaire c'est qu'aucun des deux ne semblent être adapté à comparer des MAGs; ils ne mentionnent même pas le metagenome assembly ou les MAGs dans leur publication.

dRep en revanche a justement été design pour estimer des distances Mash et rester robuste peu importe le niveau de completeness et de contamination.

Une fois qu'on a dérepliqué, on obtient nos non-redundant MAGs (nrMAGs) qu'on peut envoyer dans quant_bins (salmon) pour obtenir une belle table d'abondance (encore là je me fie sur les méthodes de Ke et al.).

Éventuellement, le module reassemble_bins pourrait être intéressant à intégrer entre dRep et Salmon; ça réassemble les reads qui mappent au nrMAGs de manière isolée pour créer des plus gros scaffolds.

[jflucier]

metawrap reassemble_bins -o BIN_REASSEMBLY -1 CLEAN_READS/ALL_READS_1.fastq -2 CLEAN_READS/ALL_READS_2.fastq -t 96 -m 800 -c 50 -x 10 -b BIN_REFINEMENT/metawrap_50_10_bins

pense pas que lon puisse faire ca sur tout les echantillon a cause quil faut pooler l'ensemble des reads pour le reassemble bins et ca va probablement buster la memoire... A tester

[jflucier]

jai concocter le recette pour builder le container qui roule dRep

commit d3aabc803d973b15df9dbfe470d0d57ad7e95442

[jflucier]

question dRep

as-t-on besoin de cet dependance partielle: gANI (aka ANIcalculator) - Performs gANI comparison method (v1.0 confirmed works)

J'ai verifier et pour le cas de dereplicate, cela s'applique au parametre:

--S_algorithm {goANI,ANIn,gANI,ANImf,fastANI} Algorithm for secondary clustering comaprisons: fastANI = Kmer-based approach; very fast ANImf = (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regions ANIn = Align whole genomes with nucmer; compare aligned regions gANI = Identify and align ORFs; compare aligned ORFS goANI = Open source version of gANI; requires nsmimscan (default: ANImf)

Si on utilise le param par defaut, pas besoin de cette dependance. Si c'est necessaire, cela va requerir bcp de debuggage.

Cet outil n'est plus maintenu depuis 2018: NOTE : Please note, as of November 2018, IMG no longer implements this version of the ANIcalculator tool.

[jflucier]

commit checkm config in drep recipe:

3a5caf0f73adfa1bf9b114a6ca09243cd258d759

[jorondo1]

question dRep

as-t-on besoin de cet dependance partielle: gANI (aka ANIcalculator) - Performs gANI comparison method (v1.0 confirmed works)

non on va utiliser le défaut, je pense que c'est Mummer la dépendance pour ça

[jflucier]

@jorondo1

finalement reussi a installer drep dans un container.

selon toi, est-ce que l'on ajoute --ignoreGenomeQuality a la commande?

On a deja roule checkM sur ces bin, selon moi ca sert a rien de rerouler ca ou minimalement faudrait le mettre ne option au script.

Pour executer drep:

module load singularity

t="/data/GQ10_bin.10.fa /data/GQ10_bin_10.fa /data/GQ10_bin.11.fa /data/GQ10_bin.12.fa /data/GQ10_bin.13.fa /data/GQ10_bin.14.fa /data/GQ10_bin.15.fa /data/GQ10_bin.16.fa /data/GQ10_bin.17.fa /data/GQ10_bin.18.fa /data/GQ10_bin.19.fa /data/GQ10_bin_19.fa"

singularity exec --writable-tmpfs \
-B /ip29/ilafores_group/projet_PROVID19/test_data:/data \
-B /ip29/ilafores_group/projet_PROVID19/assembly_60_drep:/out \
-B /nfs3_ib/ip29-ib/ssdpool/shared/ilafores_group/checkm_db:/checkm \
-e /project/def-ilafores/common/ILL_pipelines/containers/dRep.3.4.0.sif \
dRep dereplicate /out/ \
--genomes $t \
--S_algorithm fastANI \
--P_ani 0.9 --S_ani 0.99 --cov_thresh 0.30 \
--completeness 50 --contamination 5 --length 1500 \
--processors 24 --debug

Il y a encore un bug avec l'expression permettant de catcher plusieurs genome avec wild card (i.e. --genomes /data/*.fa). Je vais travaille la dessus.

Les options -B passe a singularity lui indique ou trouve les data et le output directory. Tu doit prnedre la partyi apres le ":" de l'option -B dans la commande drep.

Serait-il possible de me concocter un commande avec les options que tu desire rouler drep? Activer si possible les options suivantes accelererait cette etape selon moi:

--ignoreGenomeQuality --checkM_method taxonomy_wf

pour la doc des option, voir ici

Jai rouler drep sur 12 bin de GQ 10 et ca tourner 20 min... j'ose pas croire combien de temps ca prendra pour un screen complet. Je vais teste ca procahinement sur toute les bins de /ip29/ilafores_group/projet_PROVID19/assembly_60/refinement//metawrap_50_10_bins/.fa

Je te reviens lorsque fait

[jorondo1]

Salut! On peut mettre --ignoreGenomeQuality pour l'instant; les auteurs disent que ça pose un problème pour les microbes eukaryotes car CheckM n'est pas bon pour les évaluer, mais pour commencer c'est correct. J'en comprends que si tu actives ça, plus besoin de spécifier la methode checkm.

C'est normal que ce soit un processus assez lourd, parce qu'on fait du pairwise entre tous les bins. Au total on en a environ 2000.

Vu qu'on a 2 compartiments (fecal et saliva), peut-être que ça aiderait de faire 2 clusterings indépendants, en fonction d'une liste préparée manuellement? Ça ne serait pas vraiment utile de comparer les bins de salive avec ceux de fecal, à mon avis, même si ça peut arriver que certaines espèces se retrouvent dans les 2.

[jflucier]

oui ca devrait etre 2 analyses independante la salive et le fecal...

faudrait que tu m'indique quel echantillon appartient a quel group.

aussi... un gros ouch:

|15:08:49|jflucier@ip29:[projet_PROVID19]> ls test_data/*.fa | wc -l
2124
|15:08:57|jflucier@ip29:[projet_PROVID19]> grep -e '>' test_data/*.fa | wc -l
525832

le comparaison est sur les 525832 entrees des fasta et non sur les 2000 bins.

Si on filtre pour la longueurs des seuqence dans le bin:

|15:24:10|jflucier@ip29:[projet_PROVID19]> grep -e '>' test_data/*.fa | perl -ne '
> my @t = split("_",$_);
> #print $_;
> if($t[5] > 5000){
>     print $_;
> }
> #print $t[5] . "\n";
> ' | wc -l
203924

|15:34:55|jflucier@ip29:[projet_PROVID19]> grep -e '>' test_data/*.fa | perl -ne '
> my @t = split(":",$_);
> my @t2 = split("_",$t[1]);
> if($t2[3] > 5000){
>     print $t[0] . "\t" . $t2[3] . "\n";
> }
> #print $t[5] . "\n";
> ' | cut -d$'\t' -f 1 | sort | uniq | wc -l

2124

on filtre 60% des contigs si on considere seulement les entree avec longueur > 5000bp et on conserve l'ensemble des 2124 bins...

on filtre 80% des sequences pour > 10K bp mais on perd quelques bins (2059 / 2124).

Baser sur ces resultats, je crois que tuner l'option --length est hyper important. La valeurs par defaut de drep est a 50000.

C quoi un cutoff minimal raisonnable.

[jorondo1]

oh daym ouais c'est du stock!

La valeur de 50 000 est un minimum de longueur de génome (i-e de bin), et non une longueur de contig. Les contigs d'un même bin ne sont pas comparés entre eux, c'est le bin qui est considéré "un génome" et chaque bin est comparé aux 2123 autres pour ANI et coverage.

Les études principales auxquelles je me suis fié on conservé la valeur par défaut --length 50000, donc je ne la changerais pas.

J'oubliais, je ne suis pas sûr de comprendre ce que je dois faire avec la partie après : dans les 3 options -B, j'inquique les directory manuellement? qu'est-ce que test_data?

[jorondo1]

^^^length sera équivalent à la colonne SIZE qu'on voit dans ces fichiers: /nfs3_ib/ip29-ib/ip29/ilafores_group/projet_PROVID19/assembly_60/refinement/GQ35/metawrap_50_10_bins.stats. À mon avis la grosse majorité des bins passent ce filtre.

[jorondo1]

L'assignment des samples est ici: https://github.com/jorondo1/PROVID-19/blob/b5ef96b4ee979c296892836078e91ed352c8a62b/data/provid19_sample_metadata.csv

le mock community MSA-3001 peut être ignoré dans le processus à partir d'ici

[jflucier]

ca temrine en 7hr l'exec de drep sur les 2124 bins des echantillons /ip29/ilafores_group/projet_PROVID19/

Les resultats sont dans /ip29/ilafores_group/projet_PROVID19/assembly_60_drep

Une fois dereplique, on obtiens 1493 / 2124 bins... me semble que c'est encore bcp. Je te laisse checker ca pendant que je dev l'alignement des reads et quantification aavec salmon sur ces bins

P-e que lon doit tuner la commande drep un peu plus.... jai tourne cette commande pour le test:

newgrp def-ilafores
cd /nfs3_ib/ip29-ib/ip29/ilafores_group/projet_PROVID19
SINGULARITY_BIND=/nfs3_ib/ip29-ib/ssdpool:/ssdpool,/nfs3_ib/ip29-ib/ip29:/ip29

# --checkm_group_size default 2000 = 35G
singularity exec --writable-tmpfs \
-B /ip29/ilafores_group/projet_PROVID19/test_data:/data \
-B /ip29/ilafores_group/projet_PROVID19/assembly_60_drep:/out \
-B /ssdpool/shared/ilafores_group/checkm_db:/checkm \
-e /project/def-ilafores/common/ILL_pipelines/containers/dRep.3.4.0.sif \
dRep dereplicate /out/ \
--genomes /data/all_fastas.test.txt \
--S_algorithm fastANI \
--P_ani 0.9 --S_ani 0.99 --cov_thresh 0.30 \
--completeness 50 --contamination 10 \
--processors 24 --debug

Ca pris env 7hr d'exec pour completer.

[jflucier]

voici les stats de la job pour ref:

|08:39:33|jflucier@ip29:[assembly_60_drep]> seff 101781
Job ID: 101781
Cluster: mp2
User/Group: jflucier/def-ilafores
State: COMPLETED (exit code 0)
Nodes: 1
Cores per node: 48
CPU Utilized: 7-00:32:07
CPU Efficiency: 48.96% of 14-08:12:00 core-walltime
Job Wall-clock time: 07:10:15
Memory Utilized: 98.65 GB
Memory Efficiency: 39.30% of 251.00 GB

[jflucier]

Selon ce que je vois dans les figures produite par drep, on devrait p-e baisser le ANI du clustering secondaire a 98.5.

J'ai pris l'algo de clustering fastANI... p-e ANIn serait mieux

[jorondo1]

très cool merci

on déréplique à 0.99 pour atteindre une résolution strain-level en fait. Typiquement 0,95 c'est species level, et 0,99 (des fois même 0,999 dépendant de ta confiance dans tes bins) pour atteindre le strain level.

Est-ce que comme ça on drep les bins tous ensemble ou tu as split salive et feces?

[jflucier]

jai tout mis ensemble finalement

[jorondo1]

c'est très raisonnable 7h je m'attendais à pire!

Pour le ratio, l'étude à laquelle je me fie est passée de 11600 à 5400 avec cette procédure. Dans notre cas, on s'attend à moins de réduction, puisqu'on sait d'emblée qu'il y a deux environnements très distincts qui ne devraient pas avoir beaucoup de similarité entre eux (Salive + Feces)

pour le clustering secondaire j'irais avec ANImf, l'option par défaut, pour les mêmes raisons (les études auxquelles je me fie). Ça va être plus long, mais si on split en 2 jobs (S + F) ça devrait être raisonnable, on a 1526 bins fecal et 588 bins salive.

Voir mon comment plus haut pour la table csv qui split les échantillons

[jflucier]

ok je repars ca

[jflucier]

saliva done

job spec:

|13:50:22|jflucier@ip29:[projet_PROVID19]> seff 102237
Job ID: 102237
Cluster: mp2
User/Group: jflucier/def-ilafores
State: COMPLETED (exit code 0)
Nodes: 1
Cores per node: 48
CPU Utilized: 1-17:11:38
CPU Efficiency: 41.79% of 4-02:34:24 core-walltime
Job Wall-clock time: 02:03:13
Memory Utilized: 42.53 GB
Memory Efficiency: 16.95% of 251.00 GB
|13:50:44|jflucier@ip29:[projet_PROVID19]> 

output dir: /nfs3_ib/ip29-ib/ip29/ilafores_group/projet_PROVID19/assembly_60_drep/saliva_samples/

447 genomes dereplicated total initial genome 588

[jflucier]

stat job drep feces:

|14:53:27|jflucier@ip29:[projet_PROVID19]> seff 102238
Job ID: 102238
Cluster: mp2
User/Group: jflucier/def-ilafores
State: COMPLETED (exit code 0)
Nodes: 1
Cores per node: 48
CPU Utilized: 5-07:45:32
CPU Efficiency: 48.73% of 10-22:10:24 core-walltime
Job Wall-clock time: 05:27:43
Memory Utilized: 67.71 GB
Memory Efficiency: 26.98% of 251.00 GB

838 genomes dereplicated 1526 initial genomes

results: /ip29/ilafores_group/projet_PROVID19/assembly_60_drep/feces_samples/

[jflucier]

@jorondo1

cadeau de noel avant le temps...

bin_abundance_table.tab.txt

Reste a monter ca dans des scripts... Va faire cela au retour des fetes.

[jorondo1]

sic !!! On clanche ça en janvier. Joyeuses fêtes et merci pour ton aide !

[jflucier]

en passant le 7hr etait pour un drep fait avec l'algo fastANI

faudrait rebenchmarker avec ANImf