This tutorial will cover the complete analysis on the pangenome/panphenome of four Sinorhizobium meliloti strains, with discussions on the results and their biological interpretation. The same results can be obtained by downloading the ductape_data package (folder smeliloti); the whole set of commands used in this tutorial are provided as well. This tutorial takes for granted that you may have at least some basic experience on the usage of the command line.

#!/bin/bash

# Setup the project
dape init
dape add Rm1021 -c red
dape add BL225C -c green
dape add AK83 -c blue
dape add AK58 -c orange

# Add the genomic data and KAAS outputs
dgenome add-dir genome/
dgenome add-ko kegg/AK58.tab kegg/AK83.tab kegg/Rm1021.tab kegg/BL225C.tab

# Pangenome estimation/KEGG mapping
dgenome start -n 4 -p SinMel_

# Annotate the orthologs w/o KEGG IDs
dgenome annotate

# Statistics and exports
dgenome stats > genome_stats.tsv
dgenome export

# Add the phenomic data
dphenome add-dir phenome

# Control signal subtraction
dphenome zero

# Parameters calculation and Elbow test to decide the value of k
dphenome start -e -g
# AV/parameters calculation (k=5)
dphenome start -f -n 5

# Curve plots
dphenome plot

# Activity rings
dphenome rings
mv ActivityRing.png ActivityOverall.png
dphenome rings -o Rm1021 -d 2

# Statistics and exports (including YAML files) 
dphenome stats -a 3 -d 2 > phenome_stats.tsv
dphenome export

# Pathway maps
dape map -d 3
dape map -a

# Combined genomic/phenomic analysis
dape start -p -t 2
dape start -a -p -t 2

For further questions and updates, please contact ductape-users@googlegroups.com or subscribe to the user group

DuctApe relies on a single file (by default called ductape.db), in which all the information will be stored and from which all the outputs can be generated. The first step in the analysis is to create this file and to add the analyzed strains to it.

dape init

dape add Rm1021 -c red -n "Sinorhizobium meliloti" -d "Reference strain"
dape add BL225C -c green -n "Sinorhizobium meliloti"
dape add AK83 -c blue -n "Sinorhizobium meliloti"
dape add AK58 -c orange -n "Sinorhizobium meliloti"

The genome is added to the project file by providing a FASTA file with all the protein sequences belonging to each strain. The easiest way to add more than one genome is to prepare a directory with one file for each genome, named ORGID.faa and use the add-dir command; just substitute ORGID with the strain names provided when setting up the project.

dgenome add-dir genome

dgenome add-ko kegg/AK58.tab kegg/AK83.tab kegg/Rm1021.tab kegg/BL225C.tab

Keep in mind that other programs may be used to infer the pangenome, and their results may be imported with the add-orth command.

dgenome start -n 4 -p SinMel_

The outer circle represents a OG, each smaller circle a protein from each strain; only the gene from strain AK58 has a KEGG annotation. The annotate command will transfer that annotation to the other proteins inside the OG.

dgenome annotate

dgenome stats

pangenome_shape.png

single.png

First output table

Even though only half of the genes of each genome has been mapped to KEGG, keep in mind that many genes are not directly related to the metabolism, and are therefore impossible to be annotated through KEGG.

pangenome_stats.png

Second output table

pangenome_reaction_shape.png

Third output table

dgenome export

#re_id  name   description    pathway(s)
R04942  L-Serine hydro-lyase (adding homocysteine)    L-Serine + Selenohomocysteine <=> L-Selenocystathionine + H2O  
R01663  dCMP aminohydrolase    dCMP + H2O <=> dUMP + Ammonia  map00240,map01100
R06979  Ectoine hydro-lyase    N-gamma-Acetyldiaminobutyrate <=> H2O + Ectoine map00260,map01100,map01120,map01210
R01290  L-serine hydro-lyase (adding homocysteine; L-cystathionine-forming)    L-Serine + L-Homocysteine <=> L-Cystathionine + H2O    map00260,map00270,map01100,map01230
R00891  L-serine hydro-lyase (adding hydrogen sulfide, L-cysteine-forming)    L-Serine + Hydrogen sulfide <=> L-Cysteine + H2O      map00260

reactions_exclusive_AK83.tsv

#re_id  name  description  pathway(s)
R00489  L-aspartate 1-carboxy-lyase (beta-alanine-forming)  L-Aspartate <=> beta-Alanine + CO2  map00410,map00770,map01100,map01110
R01682  3-Sulfo-L-alanine carboxy-lyase (taurine-forming)  L-Cysteate <=> Taurine + CO2  map00430,map01100
R00374  glycine:acceptor oxidoreductase (cyanide-forming)  Glycine + 2 Acceptor <=> Cyanide + CO2 + 2 Reduced acceptor + H+  map00460,map01100
R00261  L-glutamate 1-carboxy-lyase (4-aminobutanoate-forming)  L-Glutamate <=> 4-Aminobutanoate + CO2  map00250,map00650,map01100
R01521  beta-D-Glucose:NADP+ 1-oxoreductase  beta-D-Glucose + NADP+ <=> D-Glucono-1,5-lactone + NADPH + H+  map00030,map01120
R01520  beta-D-Glucose:NAD+ 1-oxoreductase  beta-D-Glucose + NAD+ <=> D-Glucono-1,5-lactone + NADH + H+  map00030,map01120
R01154  acetyl-CoA:putrescine N-acetyltransferase  Acetyl-CoA + Putrescine <=> CoA + N-Acetylputrescine  map00330,map01100
R02466  3-Sulfino-L-alanine carboxy-lyase  3-Sulfino-L-alanine <=> Hypotaurine + CO2  map00430,map01100
R05704  glycine:acceptor oxidoreductase (hydrogen-cyanide-forming)  Glycine + 2 Acceptor <=> Hydrogen cyanide + CO2 + 2 Reduced acceptor

reactions_exclusive_AK58.tsv

#re_id  name  description  pathway(s)
R06199  None  Starch(n+1) + H2O <=> alpha-D-Glucose + Starch(n)
R07659  None  UDP-L-Ara4O + L-Glutamate <=> UDP-L-Ara4N + 2-Oxoglutarate  map00520,map01110
R01790  1,4-alpha-D-Glucan glucohydrolase  Starch + H2O <=> alpha-D-Glucose + Starch  map00500,map01100
R01791  Dextrin 6-alpha-D-glucanohydrolase  Dextrin + H2O <=> alpha-D-Glucose + Dextrin  map00500,map01100
R01741  D-Gluconate:(acceptor) 2-oxidoreductase D-Gluconic acid + FAD <=> 2-Keto-D-gluconic acid + FADH2  map00030,map01100,map01120
R01621  D-Xylulose 5-phosphate D-glyceraldehyde-3-phosphate-lyase (adding phosphate; acetyl-phosphate-forming)  D-Xylulose 5-phosphate + Orthophosphate <=> Acetyl phosphate + D-Glyceraldehyde 3-phosphate + H2O  map00030,map00680,map00710,map01100
R00749  ethanolamine ammonia-lyase (acetaldehyde-forming)  Ethanolamine <=> Acetaldehyde + Ammonia map00564

reactions_exclusive_BL225C.tsv

#re_id  name    description     pathway(s)
R00783  nitric-oxide:ferricytochrome-c oxidoreductase   Nitric oxide + H2O + Ferricytochrome c <=> Nitrite + Ferrocytochrome c + H+     map00910,map01120
R06199  None    Starch(n+1) + H2O <=> alpha-D-Glucose + Starch(n)
R04858  S-Adenosyl-L-methionine:DNA (cytosine-5-)-methyltransferase     S-Adenosyl-L-methionine + DNA cytosine <=> S-Adenosyl-L-homocysteine + DNA 5-methylcytosine   map00270,map01100
R00785  nitric-oxide:oxidized azurin oxidoreductase     Nitric oxide + H2O + Oxidized azurin <=> Nitrite + Reduced azurin + H+  map00910,map01120
R07659  None    UDP-L-Ara4O + L-Glutamate <=> UDP-L-Ara4N + 2-Oxoglutarate      map00520,map01110
R01521  beta-D-Glucose:NADP+ 1-oxoreductase     beta-D-Glucose + NADP+ <=> D-Glucono-1,5-lactone + NADPH + H+   map00030,map01120
R01520  beta-D-Glucose:NAD+ 1-oxoreductase      beta-D-Glucose + NAD+ <=> D-Glucono-1,5-lactone + NADH + H+     map00030,map01120
R07127  ATP:3-dehydro-L-gulonate 6-phosphotransferase   3-Dehydro-L-gulonate + ATP <=> 3-Dehydro-L-gulonate 6-phosphate + ADP   map00040,map00053
R04420  5-Methylthio-5-deoxy-D-ribose-1-phosphate ketol-isomerase       S-Methyl-5-thio-D-ribose 1-phosphate <=> S-Methyl-5-thio-D-ribulose 1-phosphate       map00270,map01100
[...]

reactions_exclusive_dispensable.tsv

The Phenotype Microarray data from an experiment on the four strains can be added by providing a series of csv files provided directly from the Omnilog software. The easiest way to add more than one experiment is to prepare a directory with one file for each strain, named ORGID.csv and use the add-dir command; just substitute ORGID with the strain names provided when setting up the project. For these strains all the 20 bacterial PM plates have been used (PM1 to PM20).

dphenome add-dir phenome

dphenome also automatically imports Phenotype Microarray files in YAML/JSON format

dphenome automatically detect replicas

dphenome zero

The experimental data for this well stops at 70 hours, but the plateau point is predicted to occur at 180.

In order to facilitate the comparison between distinct curves and map the compound metabolic activity inside the KEGG metabolic reconstruction, a single summarized value, termed AV (ActiVity index) is computed. Five parameters (max, area, height, lag and slope) for each curve are used to perform a k-means clustering; similar curve shapes will belong to the same cluster. Since the clusters are ranked by average area, a lower AV value will indicate a lower metabolic activity, while an higher AV value indicatates the presence of metabolic activity.

The default number of clusters used in the k-means clusterization is 10, but the optimal value can be determined performing an elbow test.

dphenome start -e -g

The -g option allows us to skipp the KEGG mapping.

elbow.png

For three parameters (max, area and height) the choice of 5 clusters gives the best improve in curves separation: we will therefore clusterize the PM data in five clusters. The minimum and maxmum AV value will therefore be 0 and 4, respectively.

dphenome start -n 5 -f

The -f option outputs a scatter plot of the five parameters used in the clusterization.

To get an overall picture of the whole Phenotype Microarray analysis, the rings command can be used; a more detailed comparison, with the single curves, can be obtained with the plot command.

dphenome rings -o Rm1021 -d 2
dphenome plot

ActivityRing.png

Each strip corresponds to a single PM curve, colored according to their AV value.

PM01_B12.png

Each strain is colored according to the dape add command.

dphenome stats -a 3 -d 2
dphenome export

ActivityCategBoxplot.png

First output table: Active wells (% of wells with activity >= 3)

Third output table: Unique metabolic functions (% of wells with AV delta >= 2 as compared to the other strains)

The curve parameters for each strain/replica, generated by the export command, may be used for further analysis; the YAML/JSON files generated may be used by other PM analysis programs, such as opm, which has additional analysis steps.

#plate_id well_id chemical category moa co_id replica activity min max height plateau slope lag area source
PM01 A01 Negative Control carbon C-Source, negative control  1 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 DuctApe
PM01 A02 L-Arabinose carbon C-Source, carbohydrate C00259 1 4 0.1 212.876351661 88.4481873788 217.820042396 8.17836437167 26.5905182274 6222.53364031 DuctApe
PM01 A03 N-Acetyl-D-Glucosamine carbon C-Source, carbohydrate C00140 1 3 0.1 175.750723251 72.2820872802 180.106135342 6.19972356011 25.4451576156 5084.32121119 DuctApe
PM01 A04 D-Saccharic acid carbon C-Source, carboxylic acid C00818 1 0 0.1 0.1 0.1 0.0 0.0 0.0 7.05 DuctApe
PM01 A05 Succinic acid carbon C-Source, carboxylic acid C00042 1 4 0.1 215.615477063 89.6254004782 225.986240731 7.30545393964 24.904552184 6305.05627692 DuctApe
PM01 A06 D-Galactose carbon C-Source, carbohydrate C00124 1 3 0.1 192.507967471 70.0603851535 196.340666445 7.67622266928 30.6220829804 4919.67860561 DuctApe
PM01 A07 L-Aspartic acid carbon C-Source, amino acid C00049 1 3 0.1 138.538227924 52.3515997237 152.92481527 4.05461063968 25.5909464203 3678.06214484 DuctApe
[...]

phenome_Rm1021.tsv

Once that both the genome and the phenome of our four strains have been analyzed, the dape module will be used to highlight the sources of genetic variability that may be related to the observed metabolic variability.

dape start -p -t 2

As suggested in the howto page, we have set a AV threshold of 2, half of the maximum value (4).

combined_variability.png (modified file)

The matrix is also saved in text format in the combined_matrix.tsv file

dape map -d 3
dape map -a

the map command will create a series of colored KEGG pathway maps to facilitate the search of putative genotype/phenotype links

pangenome_carbon/map00650.png (detail)

pangenome_carbon/map00250.png (detail)

As suggested in the I/O page, blue and orange boxes represent conserved and variable reactions, respectively; circles filled with shades of purple are compounds with significant phenotypic variability.

AK58_carbon/map00650.png and Rm1021_carbon/map00650.png(details)

AK58_carbon/map00250.png and Rm1021_carbon/map00250.png(details)

As suggested in the I/O page, colored boxes represent reactions, the yellow ones being those reactions present only in that particular strain; compounds are colored according to their AV value, from red (0 AV) to dark green (4 AV).

The start command outputs a series of other outputs, which may be useful for additional analysis steps; the metNet directory will contain a series of graph files that may be used with Gephi or Cytoscape for visualization purposes or topological studies on the metabolic pathways. Additional information may be found in the I/O page