Genetic Vectors

We address the problem of recovering regulatory networks from gene expression data. The targeted networks are directed graphs with p nodes, where each node represents a gene, and an edge directed from one gene i to another gene j indicates that gene i (directly) regulates the expression of gene j. We only consider unsigned edges; when gene i is connected to gene j, the former can be either an activator or a repressor of the latter.
The goal of (unsupervised) gene regulatory network inference is to recover the network solely from measurements of the expression of the genes in various conditions. Given the dynamic and combinatorial nature of genetic regulation, measurements of different kinds can be obtained, including steady-state expression profiles resulting from the systematic knockout or knockdown of genes or time series measurements resulting from random perturbations. In this paper, we focus on multifactorial perturbation data as generated for the DREAM4 In Silico Size 100 Multifactorial subchallenge. Multifactorial expression data are static steady-state measurements obtained by (slightly) perturbing all genes simultaneously. Multifactorial data might correspond for example to expression profiles obtained from different patients or biological replicates. Such data are easier and less expensive to obtain than knockout/knockdown or time series data and are thus more common in practice. They are however also less informative for the prediction of edge directionality [3] (link), [26] (link), [27] (link) and therefore make the regulatory network inference task more challenging.
In what follows, we define a (multifactorial) learning sample from which to infer the network as a sample of N measurements: where is a vector of expression values of all p genes in the kth experiment:
From this learning sample, the goal of network inference algorithms is to make a prediction of the underlying regulatory links between genes. Most network inference algorithms work first by providing a ranking of the potential regulatory links from the most to the less significant. A practical network prediction is then obtained by setting a threshold on this ranking. In this paper, we focus only on the first task, which is also targeted by the evaluation procedure of the DREAM4 challenge. The question of the choice of an optimal confidence threshold, although important, will be left open.
A network inference algorithm is thus defined in this paper as a procedure that exploits a LS to assign weights to putative regulatory links from any gene i to any gene j, with the aim of yielding large values for weights which correspond to actual regulatory interactions.

To determine the specificity of PopPUNK in distinguishing core sequence divergence from differences in gene content, forward-time simulations were run using Bacmeta (Sipola et al. 2018 (link)). A population of 1000 bacteria, each represented by 100 loci each 1 kb long, was simulated for 1000 generations. Insertions and deletions were fixed at a length of 100 bp. Recombinations always exchanged a complete locus and were independent of sequence divergence between donor and recipient. A sample of 25 genomes was output from the final generation of each simulation, which were analyzed using PopPUNK using default settings. Pairwise distance estimates from 50 independent simulations were then combined for plotting.
To compare PopPUNK with other clustering methods, we selected a range of previously published data sets on 10 different bacterial species (Croucher et al. 2013 (link), 2015 (link); Cohen et al. 2015 (link); Aanensen et al. 2016 (link); Grad et al. 2016 (link); Lees et al. 2016 (link), 2017 (link); Kallonen et al. 2017 (link); Koelman et al. 2017 (link); Kremer et al. 2017 (link); Alikhan et al. 2018 (link)). For each data set, as well as PopPUNK, we ran Roary (Page et al. 2015 (link)) to construct a pan-genome, using a BLAST sequence identity cutoff of 95%. We calculated core distances using the Tamura-Nei (tn93) distances (Tamura and Nei 1993 (link)) in the core genome alignment. For accessory distance, we used the Jaccard distance between the accessory gene presence/absence vectors. For comparison with another high-performance clustering algorithm, we ran RhierBAPS using between 8 and 16 cores depending on data set size (Tonkin-Hill et al. 2018 (link)). We estimated the maximum cluster size by data set, using the output from Roary and information from published analyses of these data sets.
For each species, we generated a maximum-likelihood tree from the Roary core genome alignment SNPs using IQ-TREE v1.6.3 with a GTR + I + G + ASC model (Nguyen et al. 2015 (link)). For each of these trees, we also counted polyphyly for each nonsingleton cluster. We identified all pairs of isolates from the same cluster that shared a most recent common ancestor with any isolate from a different cluster. To quantify diversity within and between clusters, we calculated the SNP distance matrices from each core alignment using pairsnp (https://github.com/gtonkinhill/pairsnp). We selected and pruned the entries in the upper triangle, which correspond to within-cluster sample comparisons, leaving the remaining entries in the upper triangle corresponding to the between-cluster distances. For E. coli and L. monocytogenes, we also performed MLST and cgMLST assignment using stringMLST (Gupta et al. 2017 (link)) and chewBBACA (Silva et al. 2018 (link)), respectively. We used the database from EnteroBase for E. coli (Alikhan et al. 2018 (link)), and from Ridom for L. monocytogenes (Ruppitsch et al. 2015 (link)). We calculated the square symmetric matrix of pairwise allelic distances for each species and each scheme which, after applying an integer allelic distance cutoff, was then used as the input to PopPUNK to produce a network in the same way as with a π-a cutoff from network refinement.