Transcriptomics Analysis - Acinetobacter baumannii

Ampicillin/sulbactam + colistin resistance

1. Introduction

Acinetobacter baumannii is a WHO priority critical pathogen for being capable of acquiring multiples mechanisms of antibiotic resistance and also being one of the biggest causes of hospital acquired infections. In this sense, alternative treatments to monotherapy, such as combination therapy are being used to tackle multidrug resistant A. baumannii. One of the last resort therapies used in hospitals is the ampicillin/sulbactam + colistin triple therapy. However, the underlying mechanisms of resistance against triple therapies remains cloudy. This work aimed to understand how A. baumannii resistance to ampicillin/sulbactam + colistin triple therapy naturally evolves from an already MDR A. baumannii in an attempt to reproduce clinical usage environment, characterizing this acquired resistance through growth parameters and molecular mechanisms under the eyes of genomics and transcriptomics.

2. Project overview

3. RNA-seq biological questions

Has the acquired resistance changed A. baumannii expression profile after 80 days of evolution? Is the expression profile of two different replicates after 16 days of evolution already different enough to determine which one thrived?

4. RNA-seq results

4.1. Pre-processing

I used an A. baumannii clinical isolate reference, but I will also use my parental strain. Kept only genes that has >=10 counts in total in at least 2 replicates and using only treatment as the analysis’ design

strain_list <- as.character(unique(coldata$strain))
for (strain in strain_list){
  coldata_temp <- coldata[coldata$strain == strain, ]
  dds_strain <- DESeqDataSetFromHTSeqCount(sampleTable = coldata_temp, 
                                  directory = here("1_Inputs_bowtie2_HTseq/counts_HTseq"), 
                                  design = ~ treatment)

  keep <- rowSums(counts(dds_strain) >= 10) >= 2
  dds_filtered <- dds_strain[keep,]
}

4.2. Sanity checks + DESeq2

4.2.1 Identifying confounding and batch factors related to extraction in different days

At first, I was analyzing all strains together within the same dataset. However, the sanity checks and metadata analysis revealed a batch effect and confounding data due to different days of extraction, making it impossible to compare different strains. Knowing that, I decided to analyze each strain separately. Since all biological replicates from each strain were extracted on the same day (strain 1 in day 1, strain 2 in day 2), it is possible to compare the treatment effect.

Figure 1. PCA containing all samples that passed the QC. It is possible to see that the control replicates have differences related to batch.

Figure 2. Sample distance heatmap containing all sample that passed the QC. It is possible to see that control replicates are not totally close to each other.

Figure 3. Variance partitioning plot addressing the observed variance in specific variables. It is possible to observe that batch variance correspond to ~15% of the observed variance.

4.2.2 Looking at each strain individually

Figure 4-7. Sample PCAs of each separate strain. It is possible to observe that the majority are well divided, except the control strain.

4.2.3 Correcting control-exclusive batch effect

To proceed with the analysis, I had to do a batch correction in the control strain using RUVSeq.

if (strain == "Control") {
    
    rep_matrix <- makeGroups(dds_filtered$treatment)
    set <- newSeqExpressionSet(as.matrix(cts))
    set_ruvs <- RUVs(set, rownames(cts), k = 2, rep_matrix)
    
    dds_filtered$W_1 <- pData(set_ruvs)$W_1
    dds_filtered$W_2 <- pData(set_ruvs)$W_2
    
    design(dds_filtered) <- ~ W_1 + W_2 + treatment
    
    vsd_ctrl <- vst(dds_filtered, blind = FALSE)
    assay(vsd_ctrl) <- removeBatchEffect(assay(vsd_ctrl), covariates = cbind(dds_filtered$W_1, dds_filtered$W_2))
    
    pcaData_ctrl <- plotPCA(vsd_ctrl, intgroup = "treatment", ntop = 500, returnData = TRUE)
    percentVar_ctrl <- round(100 * attr(pcaData_ctrl, "percentVar"))
}

After batch correction, the control samples are separaed between two groups (with treatment and without it).

Figure 8. Control sample PCA after batch correction.

4.3. DEseq2 analysis (DEGs table)

complete_results <- list()
sig_results <- list()
translate <- readRDS(here("1_Inputs_bowtie2_HTseq/rds", "annotation.rds"))
for (strain in strain_list){
  
  dds_final <- readRDS(here("1_Inputs_bowtie2_HTseq/rds", paste0("dds_final_",strain,".rds")))

# All_ and significant_gene tables
  res_temp <- results(dds_final, contrast = c("treatment", "ON", "OFF"), alpha = PADJ_CUTOFF)
  res_df <- as.data.frame(res_temp) %>% rownames_to_column("gene_id")
  res_with_names <- left_join(res_df, translate, by="gene_id")
  
  sigDE <- res_with_names %>%
    filter(padj < PADJ_CUTOFF & abs(log2FoldChange) > LFC_CUTOFF) %>%
    arrange(padj)
  
  write_csv(res_with_names, here("3_Results_bowtie2_HTseq/Gene_Tables", paste0("All_genes_", strain, ".csv")))
  write_csv(sigDE, here("3_Results_bowtie2_HTseq/Gene_Tables", paste0("Sig_genes_", strain, ".csv")))
  
  complete_results[[strain]] <- res_with_names
  sig_results[[strain]] <- sigDE
  
  saveRDS(res_with_names, file = here("1_Inputs_bowtie2_HTseq/rds", paste0("complete_results_",strain,".rds")))
  saveRDS(sigDE, file = here("1_Inputs_bowtie2_HTseq/rds", paste0("sig_genes_",strain,".rds")))
}

# Display Master Table
master_table <- dplyr::bind_rows(sig_results, .id = "Comparison")
clean_table <- master_table %>%
  dplyr::select(Comparison, gene_id, gene, product, log2FoldChange, padj) %>%
  mutate(log2FoldChange = round(log2FoldChange, 2), padj = signif(padj, 3))

write_csv(clean_table, here("3_Results_bowtie2_HTseq/Gene_Tables", "Master_Table_All_Significant_Genes.csv"))

Table 1. DEGs master table containing all significant genes along its padj and log2FC. The DEG cutoff parameters were padj <= 0.05 and log2FC of ± 0.585 (± 1.5) While the Control strain has around 100 DEGs, R1p16 and L2p16 have between 200-300 DEGs. However, even R1 having accumulated 80 days of evolution and mutations, it has only 40 DEGs, indicating that the acquired expression mutations are now constitutional. This way, it is even more important to compare R1 with the Control strain in a basal state.

4.5. Volcano plots

Figures 9-12. Volcano plots showing DEGs related to the antibiotics effect in all four strains.

4.6. Plot Counts

Figures 13-16. Plot counts of DEGs related to antibiotic response. The expression profile in R1 for this colistin stress protein greatly differs from the other strains. Also, R1 did not show any expression changes for this gene, while exhibiting a higher count.

4.7. KEGG enrichment

Figures 17-20. KEGG enrichment plot of DEGs related to antibiotic response in the different strains. While R1 and R1p16 (strains that got resistant) had energy and protein metabolisms pathways enriched, L2p16 and the Control strain had more and different pathways enriched, such as fatty acid metabolism and antibiotic resistance.

4.8. GO enrichment

Figures 21-23. GO enrichment plot of DEGs related to antibiotic response in the different strains. For this database, R1 did not present any enrichment, indicating that the strain is adapted to a point that the transcriptional change caused by the antibiotic is inherent and it is the basal state already (my hypothesis). On the other hand, R1p16 GO enrichment revealed a lot of protein metabolism pathway regulation, which can indicate a metabolic reprogramming or cell structural modification related to antibiotic stress.

5. What the R1-Control comparison can show

Since R1 strain (that survived 80 days of increasing antibiotic concentrations) is does not significantly change its transcriptome while in antibiotic stress, comparing R1 and the Control strain basal state can reveal what mechanisms are changed and always activated in the resistance strain.