Khd4 regulates subcellular protein levels of membrane-trafficking regulators

Author

Srimeenakshi Sankaranarayanan

Published

August 13, 2024

Show the code
# Libraries
library(BSgenome.Umaydis.ENSMBL.UM1) # forged ustilago genome
library(ggplot2)
library(forcats)
library(tidyr)
library(dplyr)
library(tidyverse)
library(magrittr)
library(ggpp)
library(gginnards)
library(ggrepel)
library(ggpubr)
library(ggforce)
library(ggrastr)
library(viridis)
library(reshape2)
library(gprofiler2)
library(ggvenn)
library(lookup)
library(DESeq2)
Show the code
load("C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/Kathis_lab/hyperTRIBE/hyperTRIBE.rds")

1 Find the hyphal specificity of Khd4-bound mRNA targets

1.1 Khd4-bound mRNAs responds to hyphae morphological cue

Show the code
#Load the differential expression value
#resultstable_Fi_wtvskd<- readRDS("resultstable_Fi_wtvskd.rds")
#resultstable_kd_YevFi<- readRDS("resultstable_kd_YevFi.rds")
#resultstable_wt_YevFi<- readRDS("resultstable_wt_YevFi.rds")
#resultstable_Ye_wtvkd<- readRDS("resultstable_Ye_wtvkd.rds")


#re-enter the high-conf results
#high.conf

#check the number of genes to cross check the data.frame
#n_distinct(high.conf$gene.id)

# separate gene.Ids from high confident target list for venndiagram
hcGenes=unique(high.conf$gene.id)
hcGenes=data.frame(hcGenes)
names(hcGenes)[1] <- "gene.Id"

# select all gene.Ids from wt_yvsF
#colnames(wt_yevFi)
hy.spe=dplyr::select(wt_yevFi,c("gene.ID","regulated", "up_down")) # select rows of specific columns
nhy.spe=hy.spe %>% filter(grepl('up|down', up_down))# subset only those genes that are differentially expressed
venn.hyphae=list(hyspe=nhy.spe$gene.ID, kspe=high.conf$gene.id) # make a list. 1. hyphal specific genes 2. high conf Khd4 targets

#make venndiagram
ggvenn(venn.hyphae, c("hyspe", "kspe"),fill_color = c("gray32", "darkred")) +
      ggtitle("Overlap between DEGs (wt:hyphevs yeast)  and Khd4-bound mRNAs")#+

Show the code
      #myTheme1

1.2 Analysing the AUACCC distribution in overlapping Khd4-bound mRNAs

Although Vma21 is not differentially regulated, I have added them in the following analysis as they play an important role in vacuole biogenesis.

Show the code
# collect the overlapping genes
hy.spe_TRIBE = list(shared = nhy.spe$gene.ID[nhy.spe$gene.ID %in% intersect(nhy.spe$gene.ID, high.conf$gene.id)])

df_hy.spe_TRIBE=as_tibble(hy.spe_TRIBE) # instead of as_tibble one can also use data.frame. 

############### Diagnose our shared dataset
# doing vlookup like function to identify the gene ids in highconfident targets that are overlapping with regulated gene ids in hyphal specific data (wt_yevFi)
#lookup function is "the column value to lookup" followedby the "the reference column" followed by the " column value to be entered"
#high.conf$hyphae_specific=lookup(high.conf$gene.id, nhy.spe$gene.ID, nhy.spe$gene.ID)
# Count how many value in the column hyphae specific is not NA
#n_distinct(high.conf$hyphae_specific)## looks good!!
##########################################################

# make columns with the differential gene expression data
df_hy.spe_TRIBE$log2FC_hyphae=lookup(df_hy.spe_TRIBE$shared, wt_yevFi$gene.ID, wt_yevFi$log2FoldChange.schrink)
df_hy.spe_TRIBE$log2FC_hyphae_wtvkD=lookup(df_hy.spe_TRIBE$shared, Fi_wtvskd_DEG$gene.ID, Fi_wtvskd_DEG$log2FoldChange.schrink)

# Add AUACCC position information
UTR3_AUACCC = auaccc.mot[auaccc.mot$Position == "3UTR",]
No_UTR3_AUACCC = auaccc.mot[!auaccc.mot$Gene_ID %in% UTR3_AUACCC$Gene_ID,]

# Add auaccc information
df_hy.spe_TRIBE= df_hy.spe_TRIBE %>% 
              mutate(UTR3auaccc = shared %in% UTR3_AUACCC$Gene_ID)

df_hy.spe_TRIBE$UTR3auaccc = as.factor(df_hy.spe_TRIBE$UTR3auaccc)

# manually add vma21 to the plot
df_hy.spe_TRIBE_man=rbind(df_hy.spe_TRIBE, data.frame("shared"="UMAG_11418", "log2FC_hyphae"=0.209896155, "log2FC_hyphae_wtvkD"=0.833623476,  "UTR3auaccc"= TRUE))

#scatterplot diagram
scattertest_all=ggplot()+
            geom_point(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man, size=2, color="grey")+
            stat_quadrant_counts(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man,size=6,color="darkgrey") +
            geom_point(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man[df_hy.spe_TRIBE_man$UTR3auaccc==TRUE,], size=2,color="#7AB650")+
           # stat_quadrant_counts(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man[df_hy.spe_TRIBE_man$UTR3auaccc==TRUE,],size=6, position = "jitter", color="#7AB650") +
            #scale_fill_manual(values=c("TRUE" = "darkgreen", "FALSE" = "grey"))+
            coord_cartesian(xlim=c(-2,2), ylim= c(-4,4))+
            geom_hline(yintercept=0, size=0.75, color="dark grey")+
            geom_vline(xintercept=0, size=0.5, color="dark grey")+theme_bw()+stat_cor(method ="pearson", label.x = -10, label.y=5, size=8) +
            myTheme1
scattertest_all

Show the code
scattertest_3UTRauaccc=ggplot()+
            geom_point(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man, size=2, color="grey")+
            #stat_quadrant_counts(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man,size=6,color="darkgrey") +
            geom_point(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man[df_hy.spe_TRIBE_man$UTR3auaccc==TRUE,], size=2,color="#7AB650")+
           stat_quadrant_counts(aes(x=log2FC_hyphae_wtvkD, y=log2FC_hyphae),df_hy.spe_TRIBE_man[df_hy.spe_TRIBE_man$UTR3auaccc==TRUE,],size=6, color="#7AB650") +
            #scale_fill_manual(values=c("TRUE" = "darkgreen", "FALSE" = "grey"))+
            coord_cartesian(xlim=c(-2,2), ylim= c(-4,4))+
            geom_hline(yintercept=0, size=0.75, color="dark grey")+
            geom_vline(xintercept=0, size=0.5, color="dark grey")+theme_bw()+stat_cor(method ="pearson", label.x = -10, label.y=5, size=8) +
            myTheme1
scattertest_3UTRauaccc

Note: In the publication, the quadrant count of all data points in the plot was calculated excluding the manually added vma21 to maintain the count of 122 overlapping targets and to avoid potential confusion. It is explicitly stated in the figure legends that vma21 was manually included in the analysis.

Next I would like to check the enrichment of the AUACCC motif in mRNA targets responsive to morphological changes. For this i want to use the data frames in which target mRNAs are segregated based on the position of the AUACCC motif. Before using them, i need to change the datastructure in a way that the column names are same across different data frames.

Show the code
# all genes influenced by hyphal induction"
# all genes that are diff. expressed upon hyphal induction
all_GGE_hvsY = wt_yevFi[!wt_yevFi$up_down == "not",] 
auaccc.mot$DEG_hvy= lookup(auaccc.mot$Gene_ID, all_GGE_hvsY$gene.ID, all_GGE_hvsY$gene.ID)
auaccc_DEG_all = auaccc.mot[!is.na(auaccc.mot$DEG_hvy),]
auaccc_DEG_all$L2FC_DGE_hvy = lookup(auaccc_DEG_all$Gene_ID,all_GGE_hvsY$gene.ID,all_GGE_hvsY$up_down)
auaccc_DEG_all = auaccc_DEG_all[,c(1,19,23,24)]
all_hyphae = auaccc_DEG_all[auaccc_DEG_all$L2FC_DGE_hvy=="up",]
all_yeast = auaccc_DEG_all[auaccc_DEG_all$L2FC_DGE_hvy=="down",]
Show the code
df_hy.spe_TRIBE_man <- mutate(df_hy.spe_TRIBE_man, up_down = case_when(
log2FC_hyphae  > log2(1.5) ~ "up",
log2FC_hyphae  < -log2(1.5) ~ "down",
TRUE ~ "not"
))
#make it a factor 
df_hy.spe_TRIBE_man$up_down <- factor(df_hy.spe_TRIBE_man$up_down)
Show the code
# get gene ids of ORF containing genes
uo =df_hy.spe_TRIBE_man
uo$orf = uo$shared %in% auaccc.mot[auaccc.mot$Position == "ORF",]$Gene_ID

uo = uo[uo$orf,]

# get gene ids of 5´UTR containing genes
u5 =df_hy.spe_TRIBE_man
u5$utr5 = u5$shared %in% auaccc.mot[auaccc.mot$Position == "5UTR",]$Gene_ID

u5 = u5[u5$utr5,]



# get gene ids of 3´UTR containing genes
u3 =df_hy.spe_TRIBE_man
u3$utr3 = u3$shared %in% auaccc.mot[auaccc.mot$Position == "3UTR",]$Gene_ID

u3 = u3[u3$utr3,]
Show the code
# df changes
#i will try making a function
#targets that are regulated 
### hyphal specific
getregtargets = function(targ.list, posit, morpho){
  targ.list = targ.list[,c(1:3,5)]

targ.list$position = posit # make a new column
targ_hy = targ.list[targ.list$up_down == morpho,] #
targ_hy_UR = targ_hy[targ_hy$log2FC_hyphae_wtvkD >0,]
targ_hy_DR = targ_hy[targ_hy$log2FC_hyphae_wtvkD <0,]

return(list(targ_hy_UR,targ_hy_DR))
}

# 3UTR
u3_hy  = getregtargets(u3, "3UTR", "up")
u3_hy_UR = u3_hy[[1]] # only UR genes
u3_hy_DR = u3_hy[[2]] # only DR genes

# 5UTR
u5_hy = getregtargets(u5, "5UTR", "up")
u5_hy_UR = u5_hy[[1]] # only UR genes
u5_hy_DR = u5_hy[[2]] # only DR genes


#ORF
uo_hy = getregtargets(uo, "ORF", "up")
uo_hy_UR = uo_hy[[1]] # only UR genes
uo_hy_DR = uo_hy[[2]] # only DR genes

### yeast-specific
# 3UTR
u3_ye  = getregtargets(u3, "3UTR", "down")
u3_ye_UR = u3_ye[[1]] # only UR genes
u3_ye_DR = u3_ye[[2]] # only DR genes

# 5UTR
u5_ye = getregtargets(u5, "5UTR", "down")
u5_ye_UR = u5_ye[[1]] # only UR genes
u5_ye_DR = u5_ye[[2]] # only DR genes

#ORF
uo_ye = getregtargets(uo, "ORF", "down")
uo_ye_UR = uo_ye[[1]] # only UR genes
uo_ye_DR = uo_ye[[2]] # only DR genes


#rbind UR_hyphae
ur_hyphae = rbind(u3_hy_UR,u5_hy_UR, uo_hy_UR)
dr_hyphae = rbind(u3_hy_DR, u5_hy_DR, uo_hy_DR )
ur_yeast = rbind(u3_ye_UR,u5_ye_UR, uo_ye_UR)
dr_yeast = rbind(u3_ye_DR, u5_ye_DR, uo_ye_DR )
Show the code
df_ur_hy = table(ur_hyphae$position) %>% 
            as.data.frame() %>% 
            mutate(set = "Khd4_bound_increase") %>%
            mutate(stage = "hyphae") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))

df_dr_hy = table(dr_hyphae$position) %>% 
            as.data.frame() %>% 
            mutate(set = "Khd4_bound_decrease") %>%
            mutate(stage = "hyphae") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))

df_all_hy = table(all_hyphae$Position) %>% 
            as.data.frame() %>% 
            mutate(set = "all") %>%
            mutate(stage = "hyphae") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))


df_ur_ye = table(ur_yeast$position) %>% 
            as.data.frame() %>% 
            mutate(set = "Khd4_bound_increase") %>%
            mutate(stage = "yeast") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))

df_dr_ye = table(dr_yeast$position) %>% 
            as.data.frame() %>% 
            mutate(set = "Khd4_bound_decrease") %>%
            mutate(stage = "yeast") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))

df_all_ye = table(all_yeast$Position) %>% 
            as.data.frame() %>% 
            mutate(set = "all") %>%
            mutate(stage = "yeast") %>%
            set_colnames(.,c("postn", "value", "set", "stage"))

df_postn = rbind(df_ur_hy, df_dr_hy, df_all_hy, df_ur_ye, df_dr_ye, df_all_ye)

pl_postn = ggplot(df_postn, aes(x=factor(set, levels=c("Khd4_bound_increase","Khd4_bound_decrease", "all")), y=value, fill=factor(postn, levels=c("5UTR", "ORF", "3UTR"))))+
  geom_bar(position = "fill",stat = "identity")+
  scale_fill_manual(values=c("#C8E1D3", "#7AB695", "#20854E"))+
  facet_wrap(~ stage)+
  labs(title = "position of auaccc motif in targets influenced by hyphal induction",
       x= "",
       y= "% transcripts with AUACCC",
       fill = "")+
  myTheme1+
  theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1))

pl_postn

Show the code
ggsave(pl_postn, path="C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/TRIBE_data/Editing_events_all/Khd4/bedgraphs", filename="target_hyphae_yeast_spe_motif.pdf", width=9, height=6, unit="cm")
Show the code
#pvalue calculation
df_3UTR = data.frame(URH=c(24,29),
                     DRH=c(11,35),
                     allH=c(162,435),
                     URY = c(16,23),
                     DRY = c(2,14),
                     allY=c(123,348)) %>%
                      set_rownames(.,c("with3UTRmot", "without3UTRmot"))

chiURH = fisher.test(df_3UTR[,c(1,3)])
chiURH

    Fisher's Exact Test for Count Data

data:  df_3UTR[, c(1, 3)]
p-value = 0.006893
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.198096 4.078972
sample estimates:
odds ratio 
  2.219172 
Show the code
chiURY = fisher.test(df_3UTR[,c(4,6)])
chiURY

    Fisher's Exact Test for Count Data

data:  df_3UTR[, c(4, 6)]
p-value = 0.05975
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.9373129 4.0316139
sample estimates:
odds ratio 
  1.965331 
Show the code
df_5UTR = data.frame(URH=c(11,42),
                     DRH=c(17,29),
                     allH=c(128,469),
                     URY = c(8,31),
                     DRY = c(6,10),
                     allY=c(77,394)) %>%
                      set_rownames(.,c("with5UTRmot", "without5UTRmot"))


chiDRH5 = fisher.test(df_5UTR[,c(2,3)])
chiDRH5

    Fisher's Exact Test for Count Data

data:  df_5UTR[, c(2, 3)]
p-value = 0.02618
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.069824 4.184467
sample estimates:
odds ratio 
  2.144929 
Show the code
chiDRy5 = fisher.test(df_5UTR[,c(5,6)])
chiDRy5

    Fisher's Exact Test for Count Data

data:  df_5UTR[, c(5, 6)]
p-value = 0.03898
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.8869458 9.6220041
sample estimates:
odds ratio 
  3.060614 

1.3 normalized gene counts of mRNA targets

Conducting a DESeq2 analyse to determine the mRNA levels of arl1, vma21, hok1 and mob1 in wt yeast, wt hyphae and khd4D hyphae cells.

Show the code
# make file path
mydir<-(Path="C:/Users/Sri/Documents/Khd4/DGE_analysis_SM_PhD/DGE-bioinfofolder/20220919_yeastvhyphae/")
myfiles<-list.files(path=mydir,pattern="tabular")
myfiles
 [1] "FA_1.tabular" "FA_2.tabular" "FA_3.tabular" "FK_1.tabular" "FK_2.tabular"
 [6] "FK_3.tabular" "YA_1.tabular" "YA_2.tabular" "YA_3.tabular" "YK_1.tabular"
[11] "YK_2.tabular" "YK_3.tabular"
Show the code
## make sampleTable all
# collect information for sampleTable 
sampleFiles<- myfiles
mycondition<- substr(myfiles,1,1) # extract substring in myfiles from position 1 to 2
myreplicates <- substr(myfiles, 4,4) # extract substring in myfiles from the position 4
mysamplename<- sub(".tabular", "", myfiles)

# Make sampleTable
sampleTable<- data.frame(
  sampleName= mysamplename,
  fileName= sampleFiles,
  condition=mycondition,
  replicate=myreplicates
)
#look at the sampleTable
sampleTable
   sampleName     fileName condition replicate
1        FA_1 FA_1.tabular         F         1
2        FA_2 FA_2.tabular         F         2
3        FA_3 FA_3.tabular         F         3
4        FK_1 FK_1.tabular         F         1
5        FK_2 FK_2.tabular         F         2
6        FK_3 FK_3.tabular         F         3
7        YA_1 YA_1.tabular         Y         1
8        YA_2 YA_2.tabular         Y         2
9        YA_3 YA_3.tabular         Y         3
10       YK_1 YK_1.tabular         Y         1
11       YK_2 YK_2.tabular         Y         2
12       YK_3 YK_3.tabular         Y         3
Show the code
#assign pval
mypval=0.05

1.3.1 make DEseq2 object

Show the code
# make dds object for all
dds<-DESeqDataSetFromHTSeqCount(sampleTable= sampleTable,
                                 directory = mydir,
                                 design= ~condition
                                )
# inspect DEseqDataSet (dds)
dds
class: DESeqDataSet 
dim: 6765 12 
metadata(1): version
assays(1): counts
rownames(6765): UMAG_00001 UMAG_00002 ... UMAG_15102 UMAG_15103
rowData names(0):
colnames(12): FA_1 FA_2 ... YK_2 YK_3
colData names(2): condition replicate
Show the code
#building DEseq2 model
dds <-DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and te8sting
sizeFactors(dds)
     FA_1      FA_2      FA_3      FK_1      FK_2      FK_3      YA_1      YA_2 
0.6125110 0.7001448 0.6979207 0.5899424 0.5805379 0.3910669 1.9242657 2.0003358 
     YA_3      YK_1      YK_2      YK_3 
1.5032242 1.4630655 1.4727473 1.9329907 

1.3.2 Count plots

Show the code
cnts= counts(dds, normalized=TRUE) # get normalized counts

#colnames(cnts) <- c("FA","FA","FA","FK","FK","FK","YA","YA","YA","YK","YK","YK")


dat.count=as.data.frame(subset(cnts, select=-c(10:12))) # drop yeast khd4d RNA seq data

dat.count$gene.id=rownames(dat.count)

cnts.df=melt(dat.count)

colnames(cnts.df) <- c("gene.ID", "condition", "norm.count")

my.genes<-c(arl1= "UMAG_10313", vma21="UMAG_11418", hok1="UMAG_11790")

cnts.genes=subset(cnts.df, gene.ID %in% my.genes)

cnts.genes$name=names(my.genes)[match(cnts.genes$gene.ID, my.genes)]

cnts.test1=cnts.genes %>% mutate(across("condition", str_replace, "_1", ""))
cnts.test2=cnts.test1 %>% mutate(across("condition", str_replace, "_2", ""))
cnts.test3=cnts.test2 %>% mutate(across("condition", str_replace, "_3", ""))

# change the condition name
cnts.test3 = cnts.test3 %>% mutate(condition_name =  case_when(
  condition == "FA" ~ "wildtype_hyphae",
  condition == "FK" ~ "khd4D_hyphae",
  condition == "YA" ~ "Wiltype_yeast"
))

pln=ggplot(cnts.test3, aes(x=factor(condition_name, levels = c("Wiltype_yeast", "wildtype_hyphae", "khd4D_hyphae")), y=norm.count)) +
geom_point(col="darkgreen") +
scale_y_log10()+
  facet_wrap(~name, scales="free")+
  labs(x="",
       y = "mRNA expression levels")+
  myTheme1+
  theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1 )) #facet_wrap to make adjacent plots for multiple genes.
pln

Show the code
#ggsave(pln, path="C:/Users/Srimeenakshi/Desktop/DGE_analysis_SM_PhD/TRIBE data/Editing events_all/Khd4/bedgraphs", #filename="pln.pdf", width=16, height=6, unit="cm")
Show the code
sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=English_Germany.utf8  LC_CTYPE=English_Germany.utf8   
[3] LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C                    
[5] LC_TIME=English_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] DESeq2_1.42.1                     SummarizedExperiment_1.32.0      
 [3] Biobase_2.62.0                    MatrixGenerics_1.14.0            
 [5] matrixStats_1.2.0                 lookup_1.0                       
 [7] ggvenn_0.1.10                     gprofiler2_0.2.3                 
 [9] reshape2_1.4.4                    viridis_0.6.5                    
[11] viridisLite_0.4.2                 ggrastr_1.0.2                    
[13] ggforce_0.4.2                     ggpubr_0.6.0                     
[15] ggrepel_0.9.5                     gginnards_0.2.0                  
[17] ggpp_0.5.7                        magrittr_2.0.3                   
[19] lubridate_1.9.3                   stringr_1.5.1                    
[21] purrr_1.0.2                       readr_2.1.5                      
[23] tibble_3.2.1                      tidyverse_2.0.0                  
[25] dplyr_1.1.4                       tidyr_1.3.1                      
[27] forcats_1.0.0                     ggplot2_3.5.1                    
[29] BSgenome.Umaydis.ENSMBL.UM1_2.0.0 BSgenome_1.70.2                  
[31] rtracklayer_1.62.0                BiocIO_1.12.0                    
[33] Biostrings_2.70.1                 XVector_0.42.0                   
[35] GenomicRanges_1.54.1              GenomeInfoDb_1.38.8              
[37] IRanges_2.36.0                    S4Vectors_0.40.2                 
[39] BiocGenerics_0.48.1              

loaded via a namespace (and not attached):
  [1] rstudioapi_0.16.0        jsonlite_1.8.8           ggbeeswarm_0.7.2        
  [4] GenomicFeatures_1.54.4   farver_2.1.2             rmarkdown_2.27          
  [7] ragg_1.3.2               zlibbioc_1.48.0          vctrs_0.6.5             
 [10] memoise_2.0.1            Rsamtools_2.18.0         RCurl_1.98-1.13         
 [13] rstatix_0.7.2            progress_1.2.3           htmltools_0.5.8.1       
 [16] S4Arrays_1.2.0           polynom_1.4-1            curl_5.2.1              
 [19] broom_1.0.6              SparseArray_1.2.3        htmlwidgets_1.6.4       
 [22] plyr_1.8.9               cachem_1.0.8             plotly_4.10.4           
 [25] GenomicAlignments_1.38.0 mime_0.12                lifecycle_1.0.4         
 [28] pkgconfig_2.0.3          Matrix_1.6-4             R6_2.5.1                
 [31] fastmap_1.1.1            GenomeInfoDbData_1.2.11  shiny_1.8.1.1           
 [34] digest_0.6.33            colorspace_2.1-0         AnnotationDbi_1.64.1    
 [37] textshaping_0.4.0        crosstalk_1.2.1          RSQLite_2.3.7           
 [40] labeling_0.4.3           filelock_1.0.3           fansi_1.0.6             
 [43] timechange_0.3.0         httr_1.4.7               polyclip_1.10-6         
 [46] abind_1.4-5              compiler_4.3.1           bit64_4.0.5             
 [49] withr_3.0.0              backports_1.5.0          BiocParallel_1.36.0     
 [52] carData_3.0-5            DBI_1.2.3                biomaRt_2.58.2          
 [55] ggsignif_0.6.4           MASS_7.3-60              rappdirs_0.3.3          
 [58] DelayedArray_0.28.0      rjson_0.2.21             tools_4.3.1             
 [61] vipor_0.4.7              beeswarm_0.4.0           httpuv_1.6.15           
 [64] glue_1.6.2               restfulr_0.0.15          promises_1.3.0          
 [67] generics_0.1.3           gtable_0.3.5             tzdb_0.4.0              
 [70] data.table_1.15.4        hms_1.1.3                xml2_1.3.6              
 [73] car_3.1-2                utf8_1.2.4               pillar_1.9.0            
 [76] later_1.3.2              tweenr_2.0.3             BiocFileCache_2.10.2    
 [79] lattice_0.22-5           bit_4.0.5                tidyselect_1.2.1        
 [82] locfit_1.5-9.10          knitr_1.47               gridExtra_2.3           
 [85] xfun_0.45                stringi_1.8.3            lazyeval_0.2.2          
 [88] yaml_2.3.7               evaluate_0.24.0          codetools_0.2-20        
 [91] cli_3.6.1                systemfonts_1.1.0        xtable_1.8-4            
 [94] munsell_0.5.1            Rcpp_1.0.12              dbplyr_2.5.0            
 [97] png_0.1-8                XML_3.99-0.16.1          parallel_4.3.1          
[100] blob_1.2.4               prettyunits_1.2.0        bitops_1.0-7            
[103] scales_1.3.0             crayon_1.5.3             rlang_1.1.1             
[106] KEGGREST_1.42.0          cowplot_1.1.3