Conditional GWAS

Overview

Models

Here we present the results of a genome-wide association study (GWAS) conducted using pyseer (John A. Lees et al. 2018), a tool for identifying genetic associations with complex traits. The GWAS using the linear mixed effects model such that the minimum inhibitory concentration (MIC) of Penicillin G: $$\mathbf{y}\sim \mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{K}\mathbf{u}+ϵ$$

with

$$\begin{array}{rl}\mathbf{u}\sim \mathcal{N}(0,{\sigma }_g^2\mathbf{G}),& ϵ\sim \mathcal{N}(0,{\sigma }_g^2\mathbf{I})\end{array}$$

Where $\mathbf{G}$ is the similarity matrix that was obtained from the phylogenetic relationship of the samples, and $\mathbf{I}$ is the identity matrix.

Table 1: Information about the GWAS model

Symbol	Meaning
𝑦	A vector containing the logarithm of the MIC data for each sample.
𝑊	A design matrix containing the conditional covariates: We used a categorical variable for each unique combination of pencillin binding proteins pbp1a, pbp2b and pbp2x, coded according to the allele scheme proposed in (Li et al. 2016), and benchmarked in (Li et al. 2017). The allele codes were assigned using an inhouse script, but should be comparable (although with different names) to those of Pathogenwatch pipeline.
𝑎	Fixed effects for the covariates.
𝑋	A design matrix containing the presence or absence of the unitigs.
𝑏	Fixed effects of the unitigs.
𝐾	The similarity matrix between all pairs of samples. We used a neighbor-joining phylogeny estimated using PopPUNK (John A. Lees et al. 2019).
𝑢	Random effects for each row of the kinship matrix.
ε	Normally distributed error terms.

Thus we assume that the MIC can be predicted by the covariates with the correlation structure implied by the phylogenetic relationship of the samples. If we take the expectation of the equation we obtain $$E(\mathbf{y})\sim \mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{K}\mathbf{u}+ϵ$$ $$=E(\mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{K}\mathbf{u}+ϵ)$$ $$=\mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}$$ and variance $$V(\mathbf{y})=V(\mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{K}\mathbf{u}+ϵ)$$ $$\begin{array}{r}V(\mathbf{y})=V(\mathbf{W}\mathbf{a}+\mathbf{X}\mathbf{b}+\mathbf{K}\mathbf{u}+ϵ)\\ =V(\mathbf{K}\mathbf{u}+ϵ)\\ =V(\mathbf{K}\mathbf{u})+V(ϵ)\\ ={\sigma }_g^2\mathbf{G}+{\sigma }_g^2\mathbf{I}\end{array}$$

Filtering of unitigs and defining significance threshold

We filtered the unitigs based on a minimum allele frequency (MAF) between [0.01, 0.99] meaning that the unitigs have to be present in more than 1% and less than 99% of the isolates. We used the number of patterns in the unitigs to define a significance threshold of the p-values of the unitigs. We detected 974212 unique patterns translating to a significance threshold of $p=5.13\cdot {10}^{-08}$. All downstream results and visualizations were reported for unitigs with a smaller p-value than the threshold. Many of these hits were then annotated using the nomenclature of Bakta (Schwengers et al. 2021). The annotations can be looked up on NCBI for further information.

Results

Significant hits ordered by effect size

In the figures below the unitigs are mapped to reference genomes and. As references we used multiple genomes, and additionally we annotated our

Figure 1: Interactive plot of unitigs mapped to annotated regions of references and assemblies. This scatter plot shows the relationship between the average effect size and the maximum -log10(p-value) of genes associated with penicillin resistance. The size of the points represents the number of unitigs, and the color represents the average minor allele frequency (MAF). Points are labeled with gene name. The plot highlights genes that have high effect size and large -log10(p-value) as well as high unitig count and high MAF value.

Summary of annotated regions matching the unitigs (Table)

Show 102550100 entries

Search:

test

	gene	hits	maxp	avg_maf	avg_beta	names_for_display

	gene	hits	maxp	avg_maf	avg_beta	names_for_display
1	cds-AFC94349.1	1	27.0931264652779	0.01	6
2	cds-AFC94350.2	1	27.0931264652779	0.01	6
3	cds-AFC94371.1	1	27.0931264652779	0.0109	6
4	cds-USV25787.1	3	27.0931264652779	0.0105666666666667	4.84
5	cds-AFC94373.1	2	27.0931264652779	0.01135	4.305
6	cds-SPNHU17_RS03300	11	27.0931264652779	0.0123454545454545	3.71818181818182
7	cds-WP_000766577.1|cds-SPNHU17_RS03300	1	16.5622494371796	0.0116	3.6
8	cds-WP_001032463.1	1	7.45099673797421	0.017	3.51
9	cds-USV25279.1	3	27.0931264652779	0.0130666666666667	3.45666666666667
10	pabB	31	27.0931264652779	0.0127258064516129	3.38454838709677	pabB

Showing 1 to 10 of 57 entries

Previous123456Next

Figure 2: Gene hits Annotated regions of mapped unitigs. The table summarizes the locations of unitigs in the genome, including the gene annotated region (gene), the number of unitigs in that region (HITS), the highest observed p-value (MAXP), the average minimum allele frequency (AVG_MAF), and the average effect size (AVG_BETA) of unitigs within the region. The table also includes a column (NAMES_FOR_DISPLAY) to indicate regions that have common gene names, which match the names in the interactive figures above.

Unitig sequences, significance, effect size and heritability (Table)

Show 102550100 entries

Search:

	Unitigs	af	filter-pvalue	lrt-pvalue	beta	beta-std-err	variant_h2	annotation

	Unitigs	af	filter-pvalue	lrt-pvalue	beta	beta-std-err	variant_h2	annotation
1	GAGATTCCAGTCAAGGAAGTTCGTGATGAAATGATTG	0.981	3.14e-37	1.5e-8	-1.64	0.29	0.0677	NZ_CP020549.1:840583-840619;sufD;sufD;sufD
2	CTGCTAATTTAAGCGGGCATCCTGATAGAACCAAAGCTCGAAGA	0.0124	4.38e-10	3.45e-14	2.51	0.331	0.0906	NZ_CP020549.1:670101-670144;cds-WP_000766577.1;cds-WP_000766577.1;cds-WP_000766577.1
3	CTAACGTAGTAGAAAAAAACTCCTCCGCCTAAGAC	0.846	0	1.57e-8	-2.97	0.525	0.0676	NZ_CP020549.1:390377-390411;pbp1a;pbp1a;pbp1a
4	GCATGGTCAATATCAACACAAAAAACAATTGTTTT	0.99	0.00125	1.22e-8	-0.651	0.114	0.0681	NZ_CP020549.1:523950-523984;cds-WP_000229101.1;cds-WP_000229101.1;cds-WP_000229101.1
5	ACTGAGAGGTGTTAATACTTGACGATTGACTTCTATC	0.0124	4.94e-10	1.58e-15	2.84	0.356	0.0952	NZ_CP020549.1:669462-669498;pabB;pabB;pabB
6	CATAGGCATTGTAGCCCGCCTCCTGCTCTATC	0.0108	1.94e-12	8.07e-28	6	0.547	0.13	NZ_CP020549.1:668601-668632;pabB;pabB;pabB
7	CACAATCATCATGTTTTCAGAGCGATTTTTGGGGTCCTGTTC	0.0124	1.11e-9	1.58e-15	2.84	0.356	0.0952	NZ_CP020549.1:668779-668820;pabB;pabB;pabB
8	AAGGTCAAGACTTTCTCTTCTACCTGTCCTCTTTCCAGCAAATGCTGACGGTAAATTCCTG	0.0124	1.81e-9	1.58e-15	2.84	0.356	0.0952	NZ_CP020549.1:669716-669776;pabB;pabB;pabB
9	TTGAGCGATAGCTTCTTTATTGGCAACAGGGAT	0.988	2.9e-15	2.43e-8	-1.12	0.2	0.0667	NZ_CP020549.1:1809049-1809081;groL;groL;groL
10	GTCAATATCAACACAAAAAACAATTGTTTTAT	0.99	0.000934	2.14e-8	-0.633	0.113	0.067	NZ_CP020549.1:523955-523986;cds-WP_000229101.1;cds-WP_000229101.1;cds-WP_000229101.1

Showing 1 to 10 of 394 entries

Previous12345…40Next

Figure 3: Unitig sequences. The table displays various characteristics of the significant unitigs including allele frequency (AF), p-value (filter-pvalue), likelihood (lrt-pvalue), effect size (beta), estimated standard error of the effect size (beta-std-err), estimated heritability of the unitig sequence (variant_h2), and the annotated region in which the unitig is found (annotation) as well as the nearest upstream and downstream annotated regions (gene-upstream gene-downstream gene).

Unitigs mapped to reference genomes

In total we extracted 394 significant unitigs.

Note: In past version of this analysis we found 2527 significant unitigs. In the current analysis we have conditioned on all PBPs and assigned new codes to unobserved alleles. This removes most significant hits, demonstrating the importance of properly accounting for PBP composition in the data.

We mapped these unitigs to different reference genomes, summarized in the table below. Streptococcus pneumoniae has a large pan-genome, and none of the references contain all of the unitigs.

Table 2: Information about the reference strains

Reference strain	Description
ATCC70066	Multidrug-Resistant Pandemic Clone Streptococcus pneumoniae. 106 out of 394 unitigs were mapped to this genome. https://doi.org/10.1128/JB.01343-08
Hu17 - NCBI reference.	Current NCBI reference genome. 116 out of 394 unitigs were mapped to this genome. https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP020549.1
Hu17	A high-level beta-lactam resistance and penicillin sensitive phenotype. 116 out of 394 unitigs were mapped to this genome. https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_002076835.1/
EF3030	A serotype 19F isolate that colonizes the nasopharynx of mice while being mostly noninvasive. 144 out of 394 unitigs were mapped to this genome. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6509521/
ST556	A Multidrug-Resistant Isolate from an Otitis Media Patient. 124 out of 394 unitigs were mapped to this genome. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3370858/
TIGR4	From the paper: The acquisition of clinically relevant amoxicillin resistance in Streptococcus pneumoniae requires ordered horizontal gene transfer of four loci. 143 out of 394 unitigs were mapped to this genome. https://pubmed.ncbi.nlm.nih.gov/35877768/
SN75752	From the paper: The acquisition of clinically relevant amoxicillin resistance in Streptococcus pneumoniae requires ordered horizontal gene transfer of four loci. 135 out of 394 unitigs were mapped to this genome. https://pubmed.ncbi.nlm.nih.gov/35877768/

ATCC700669

Hu17 - reference

Hu17

EF3030

ST556

TIGR4

SN75752

Figure 4: Unitigs mapped to various reference genomes. The reference strains used for the mapping are described in the table above. The data displayed when hovering over the points in the figure represents the annotation attributes extracted from the genes that overlap with the starting position of the unitig (basepair position). It’s noteworthy that the “product” and “gene” attributes can provide valuable insights into the function of the genomic regions that the unitigs overlap with.

References

Lees, John A, Marco Galardini, Stephen D Bentley, Jeffrey N Weiser, and Jukka Corander. 2018. “Pyseer: A Comprehensive Tool for Microbial Pangenome-Wide Association Studies.” Edited by Oliver Stegle. Bioinformatics 34 (24): 4310–12. https://doi.org/10.1093/bioinformatics/bty539.

Lees, John A., Simon R. Harris, Gerry Tonkin-Hill, Rebecca A. Gladstone, Stephanie W. Lo, Jeffrey N. Weiser, Jukka Corander, Stephen D. Bentley, and Nicholas J. Croucher. 2019. “Fast and Flexible Bacterial Genomic Epidemiology with PopPUNK.” Genome Research 29 (2): 304–16. https://doi.org/10.1101/gr.241455.118.

Li, Yuan, Benjamin J. Metcalf, Sopio Chochua, Zhongya Li, Robert E. Gertz, Hollis Walker, Paulina A. Hawkins, et al. 2016. “Penicillin-Binding Protein Transpeptidase Signatures for Tracking and Predicting β-Lactam Resistance Levels in Streptococcus Pneumoniae.” Edited by James M. Hughes. mBio 7 (3). https://doi.org/10.1128/mbio.00756-16.

Li, Yuan, Benjamin J. Metcalf, Sopio Chochua, Zhongya Li, Robert E. Gertz, Hollis Walker, Paulina A. Hawkins, Theresa Tran, Lesley McGee, and Bernard W. Beall. 2017. “Validation of β-Lactam Minimum Inhibitory Concentration Predictions for Pneumococcal Isolates with Newly Encountered Penicillin Binding Protein (PBP) Sequences.” BMC Genomics 18 (1). https://doi.org/10.1186/s12864-017-4017-7.

Schwengers, Oliver, Lukas Jelonek, Marius Alfred Dieckmann, Sebastian Beyvers, Jochen Blom, and Alexander Goesmann. 2021. “Bakta: Rapid and Standardized Annotation of Bacterial Genomes via Alignment-Free Sequence Identification.” Microbial Genomics 7 (11). https://doi.org/10.1099/mgen.0.000685.