5 PCB ~ Gene Expression

5.1 Libraries

library(tidyverse)
library(edgeR)
library(ggfortify)
library(ggrepel)
library(gplots)
library(rbioapi)
library(WebGestaltR)
library(kableExtra)
library(knitr)

5.2 Data

meta <- 
  read.csv("../hg-rna/Input Files/metadata.csv") %>% 
  filter(!sex=="male") %>% 
  mutate(sample=gsub("Hg","", sample),
         sample=gsub("a", "", sample),
         sample=gsub("b", "", sample),
         sample=gsub("11", "1", sample)) 

pcb_maxlod <- 
  read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>% 
  dplyr::rename(sample=c(1)) %>% 
  mutate(across(c(2:17), \(x) as.numeric(x))) 

sampleinfo <- 
  merge(meta, pcb_maxlod, by="sample") %>% 
  mutate(iav=ifelse(disease_stage=="acute" | disease_stage=="peak", "Positive", "Negative"),
         pcbbin=ifelse(pcbbin=="1", "Present", "Absent"))

dge_matrix <-
  read.csv("../hg-rna/Output Files/txome_genecounts_locid_nohemo.csv", row.names=1) %>%
  rename_with(~ str_remove(., "Hg"), everything()) %>% 
  rename_with(~ str_remove(., "a"), everything()) %>%
  rename_with(~ str_remove(., "b"), everything()) %>% 
  rename("1513" = "11513") %>% 
  select(all_of(sampleinfo$sample)) 
  

# Descriptive stats for DGE ~ PCBs
table(sampleinfo$pcbbin)

## 
##  Absent Present 
##      20       3

table(sampleinfo$year)

## 
## 2016 2017 2018 2019 
##    2    2   13    6

table(sampleinfo$sex)

## 
##  female unknown 
##      22       1

table(sampleinfo$location)

## 
##  monomoy muskeget 
##       21        2

table(sampleinfo$molt)

## 
##  U  V 
##  1 22

table(sampleinfo$iav, sampleinfo$pcbbin)

##           
##            Absent Present
##   Negative     15       1
##   Positive      5       2

5.3 Configure data

5.3.1 Create DGE Object

groups <- sampleinfo$pcbbin
DGE <- DGEList(counts=dge_matrix, group=groups)
DGE$samples

##        group lib.size norm.factors
## 1256  Absent  4573826            1
## 1261  Absent  1943267            1
## 1262  Absent  4653972            1
## 1265 Present  2979315            1
## 1272  Absent  4708624            1
## 1273  Absent  5011031            1
## 1282 Present  2951861            1
## 1283  Absent  4658064            1
## 1288  Absent  2823770            1
## 1351  Absent  1430600            1
## 1354  Absent  3618804            1
## 1380  Absent  3204059            1
## 1387  Absent  4158987            1
## 1397  Absent  4343361            1
## 1403  Absent  6100522            1
## 1404  Absent  3533565            1
## 1421  Absent  4508013            1
## 1513 Present  2842670            1
## 1516  Absent  4688415            1
## 1527  Absent  3763080            1
## 1528  Absent  2894695            1
## 861   Absent  5407908            1
## 862   Absent  4075756            1

5.3.2 Filter Genes with Low Counts

• Keep rows that have worthwhile counts in a min # of samples (= smallest group size = 3)
• 17,487 genes -> 12,067 genes

keep <- filterByExpr(DGE)
DGE <- DGE[keep, , keep.lib.sizes=FALSE]

write.csv(DGE$counts, "Output Files/pcb_genes_postfilter.csv")

5.3.3 Normalize, Estimate Dispersion

DGE <- calcNormFactors(DGE)

# Overall, common dispersion
DGE <- estimateCommonDisp(DGE, verbose = TRUE)

## Disp = 0.17147 , BCV = 0.4141

# Dispersion trend based on gene abundance
DGE <- estimateTrendedDisp(DGE)

# Tagwise dispersion
DGE <- estimateTagwiseDisp(DGE, verbose=TRUE)

## Using interpolation to estimate tagwise dispersion.

# Plot Dispersion
plotBCV(DGE)

# Did normalization work? 
CPM <- 
  cpm(DGE, normalized.lib.sizes = TRUE, log = TRUE)
boxplot(CPM, las = 2, ylab = "log2 CPM", main = "Normalized Data")

write.csv(CPM, "Output Files/pcb_genes_postfilter_norm.csv")

5.4 Exploratory Plots

rawdist <- dist(t(CPM), method="euclidean")
plot(hclust(rawdist, method = "average"), xlab="Average Euclidean Distance")

PCA <- prcomp(t(rawdist))

autoplot(PCA, data=sampleinfo, color="pcbbin") +
  scale_color_manual(values=c("#b4b4b4", "#61dfea")) +
  theme_classic() +
  geom_text_repel(label=sampleinfo$sample, 
            size=3)

5.4.1 DGE Analysis

• DE = uncorrected p < 0.05 & absolute log2FC > 1

et <- exactTest(DGE, pair=c("Absent", "Present"))

results <- 
  as.data.frame(topTags(et, n=dim(D)))

de <- 
  results %>% 
  filter(PValue < 0.05 & abs(logFC) > 1)

write.csv(de, "Output Files/pcb_degenes.csv")

5.4.2 Heatmap

detags <- rownames(de)

# Add +/- to sample names to indicate IAV status
sampiav <- 
  sampleinfo %>% 
  select(sample, iav) %>% 
  mutate(sample=ifelse(iav=="Positive", 
                       paste0(sample, " (+)"), 
                       paste0(sample, " (-)")))

CPMiav <-
  CPM %>% 
  as.data.frame() %>%
  setNames(sampiav$sample) %>% 
  as.matrix()

# Heatmap of top 100 DE genes (50 up and 50 down)
topde <-
  data.frame(rbind(slice_min(de, logFC, n=50), 
                   slice_max(de, logFC, n=50))) %>% 
  rownames_to_column(var="gene")

jpeg(filename="Figures/heatmap_top100.jpeg", width=8, height=5, units="in", res=300)
heatmap.2(as.matrix(CPMiav[rownames(CPMiav) %in% topde$gene,]),
          scale="row", 
          trace="none", 
          dendrogram="column", 
          key=FALSE, 
          labRow=FALSE,
          lmat=rbind(c(4,3), c(2,1)),
          lwid=c(0.05,1),
          lhei=c(0.2, 1), 
          colsep=c(2,5),
          margins=c(5,2.5))
dev.off()

## png 
##   2

## All genes - gene expression patterns not clear
heatmap.2(as.matrix(CPMiav),
          scale="row", 
          trace="none", 
          dendrogram="column", 
          key=FALSE, 
          labRow=FALSE,
          lmat=rbind(c(4,3), c(2,1)),
          lwid=c(0.05,1),
          lhei=c(0.2, 1), 
          colsep=c(2,5),
          margins=c(5,2.5))

5.4.3 Heatmap - Presentation

jpeg(filename="Figures/heatmap_top100_pres.jpeg", width=8, height=5, units="in", res=300)
heatmap.2(as.matrix(CPM[rownames(CPM) %in% topde$gene,]),
          scale="row", 
          trace="none", 
          dendrogram="column", 
          key=TRUE, 
          labRow=FALSE,
           density.info="none",
          lmat = rbind(c(0,3), c(2,1), c(0,4)),
          lhei = c(0.2, 0.5, 0.2),
          lwid = c(0.05, 1),
          # lmat=rbind(c(4,3), c(2,1)),
          # lwid=c(0.05,1),
          # lhei=c(0.2, 1), 
          colsep=c(2,5),
          margins=c(5,2.5),
          xlab = "Gray seal pup",
          ylab = "Genes")
dev.off()

## png 
##   2

5.5 Gene Ontology Analysis

5.5.1 Configure data

up <- 
  de %>% 
  filter(logFC>0) %>% 
  rownames()

# test to make sure up-regulated genes are up in pups with PCBs
dge_matrix %>% 
  filter(rownames(.) %in% up) %>% 
  t() %>%
  data.frame() %>% 
  select(ACP5, CXCL8) %>% 
  mutate(pcb=sampleinfo$pcbbin) %>% 
  pivot_longer(cols=c(1:2), names_to="gene") %>% 
  ggplot(., aes(x=pcb, y=value)) +
  geom_boxplot() +
  facet_wrap(~gene, scales="free_y")

down <- 
  de %>% 
  filter(logFC<0) %>% 
  rownames()

background <- 
  rownames(DGE$counts)

write.table(background, "Output Files/deg_postfilter.csv",
          quote=FALSE, row.names=FALSE, col.names=FALSE)

5.5.2 Up-regulated

up_bp <- 
  rba_panther_enrich(genes=up, organism=9606, 
                     annot_dataset="GO:0008150",
                     test_type="FISHER", correction="FDR", cutoff=0.05,
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="bp", 
         direction="up")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 318 genes from `organism 9606` against `GO:0008150` datasets.

up_mf <- 
  rba_panther_enrich(genes=up, organism=9606, 
                     annot_dataset="GO:0003674",
                     test_type="FISHER", correction="FDR", cutoff=0.05,
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="mf", 
         direction="up")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 318 genes from `organism 9606` against `GO:0003674` datasets.

up_cc <- 
  rba_panther_enrich(genes=up, organism=9606, 
                     annot_dataset="GO:0005575",
                     test_type="FISHER", correction="FDR", cutoff=0.05, 
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="cc", 
         direction="up")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 318 genes from `organism 9606` against `GO:0005575` datasets.

5.5.3 Down-regulated

down_bp <- 
  rba_panther_enrich(genes=down, organism=9606, 
                     annot_dataset="GO:0008150",
                     test_type="FISHER", correction="FDR", cutoff=0.05,
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="bp", 
         direction="down")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 256 genes from `organism 9606` against `GO:0008150` datasets.

down_mf <- 
  rba_panther_enrich(genes=down, organism=9606, 
                     annot_dataset="GO:0003674",
                     test_type="FISHER", correction="FDR", cutoff=0.05,
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="mf", 
         direction="down")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 256 genes from `organism 9606` against `GO:0003674` datasets.

down_cc <- 
  rba_panther_enrich(genes=down, organism=9606, 
                     annot_dataset="GO:0005575",
                     test_type="FISHER", correction="FDR", cutoff=0.05, 
                     ref_genes=background, ref_organism=9606)$result %>% 
  mutate(go="cc", 
         direction="down")

## Performing PANTHER over-representation analysis (Fisher's exact test) on 256 genes from `organism 9606` against `GO:0005575` datasets.

5.5.4 Compile Results

Removed cell periphery - found over-represented among up and down regulated genes.

go_all <-
  rbind(up_bp, up_mf, up_cc,
      down_bp, down_mf, down_cc) %>% 
  filter(plus_minus=="+") %>% 
  select(fold_enrichment, number_in_list, 
         term.id, term.label, 
         go, direction) %>% 
  filter(!term.label=="cell periphery")

go_all %>% 
  kable() %>% 
  kable_styling("basic")

fold_enrichment	number_in_list	term.id	term.label	go	direction
5.734258	14	GO:0030218	erythrocyte differentiation	bp	up
2.134628	49	GO:0006952	defense response	bp	up
5.247387	14	GO:0034101	erythrocyte homeostasis	bp	up
2.407892	38	GO:0098542	defense response to other organism	bp	up
2.112821	46	GO:0051707	response to other organism	bp	up
2.110381	46	GO:0043207	response to external biotic stimulus	bp	up
1.994489	50	GO:0044419	biological process involved in interspecies interaction between organisms	bp	up
3.648693	18	GO:0051607	defense response to virus	bp	up
19.865108	5	GO:0046501	protoporphyrinogen IX metabolic process	bp	up
2.039721	46	GO:0009607	response to biotic stimulus	bp	up
9.932554	7	GO:0089718	amino acid import across plasma membrane	bp	up
26.486811	4	GO:0015669	gas transport	bp	up
4.213811	14	GO:0002262	myeloid cell homeostasis	bp	up
2.332498	30	GO:0140546	defense response to symbiont	bp	up
1.304037	152	GO:0050896	response to stimulus	bp	up
39.730216	3	GO:0015670	carbon dioxide transport	bp	up
1.512108	91	GO:0006950	response to stress	bp	up
8.179750	7	GO:0006778	porphyrin-containing compound metabolic process	bp	up
7.946043	7	GO:0045071	negative regulation of viral genome replication	bp	up
9.932554	6	GO:0006779	porphyrin-containing compound biosynthetic process	bp	up
9.932554	6	GO:0033014	tetrapyrrole biosynthetic process	bp	up
19.865108	4	GO:0006785	heme B biosynthetic process	bp	up
19.865108	4	GO:0046492	heme B metabolic process	bp	up
2.293703	28	GO:0045087	innate immune response	bp	up
17.657874	4	GO:0006782	protoporphyrinogen IX biosynthetic process	bp	up
7.131064	7	GO:0048821	erythrocyte development	bp	up
2.883645	18	GO:0048872	homeostasis of number of cells	bp	up
5.885958	8	GO:0140374	antiviral innate immune response	bp	up
2.860575	18	GO:0030099	myeloid cell differentiation	bp	up
29.797662	3	GO:0072488	ammonium transmembrane transport	bp	up
2.755015	19	GO:0009615	response to virus	bp	up
3.844860	12	GO:2001236	regulation of extrinsic apoptotic signaling pathway	bp	up
6.621703	7	GO:0033013	tetrapyrrole metabolic process	bp	up
5.297362	8	GO:0045069	regulation of viral genome replication	bp	up
39.730216	6	GO:0170014	ankyrin-1 complex	cc	up
33.108513	5	GO:0014731	spectrin-associated cytoskeleton	cc	up
1.502012	113	GO:0005886	plasma membrane	cc	up
4.786773	10	GO:0030863	cortical cytoskeleton	cc	up
2.255427	26	GO:0098794	postsynapse	cc	up
3.730502	22	GO:0002250	adaptive immune response	bp	down
15.097048	10	GO:0042101	T cell receptor complex	cc	down
5.283967	14	GO:0098802	plasma membrane signaling receptor complex	cc	down
3.071306	21	GO:0098797	plasma membrane protein complex	cc	down
3.400062	17	GO:0043235	receptor complex	cc	down
35.100636	3	GO:0071821	FANCM-MHF complex	cc	down
15.600283	4	GO:0044194	cytolytic granule	cc	down
6.240113	6	GO:0001772	immunological synapse	cc	down

write.csv(go_all, "Output Files/goterm.csv",
          row.names=FALSE)

5.5.5 Plot

goplotdata <- 
  go_all %>% 
  mutate(fold_enrichment=ifelse(direction=="down", -fold_enrichment, 
                                fold_enrichment)) %>%
  arrange(fold_enrichment) %>% 
  mutate(term.label=factor(term.label, levels=term.label),
         go=ifelse(go=="bp", "Biological Process", 
                   ifelse(go=="cc", "Cellular Component", "Molecular Function")))

go_plot <-
  ggplot(goplotdata, aes(x=term.label, y=fold_enrichment, fill=direction)) +
  geom_bar(stat = "identity", position="dodge") +
  scale_fill_manual(values =c(up = "#7eaaac", 
                              down = "#b4b4b4")) +
  facet_grid(vars(go), scales="free_y", space="free_y") +
  labs(x=NULL, y="Fold Enrichment") +
  scale_y_continuous(limits=c(-40,45), expand=c(0,0), n.breaks=10) +
  coord_flip() +
  theme_bw() +
  theme(panel.grid=element_blank(), text=element_text(size=13),
        legend.position = "none")
go_plot

ggsave("Figures/go_foldenrichment.jpeg", go_plot, width=10, height=9, units="in")

5.5.6 Production figure

go_plot_prod <-
  go_plot + 
  theme(text = element_text(size = 9))

ggsave("Figures/Figure_3.jpeg", go_plot_prod, 
        width = 7480, height = 6732, units = "px", dpi = 1000)

5.6 KEGG Pathway Analysis

KEGG <- WebGestaltR(enrichMethod = "ORA", organism = "hsapiens", 
                    enrichDatabase = "pathway_KEGG", interestGene = rownames(de),
                    interestGeneType = "genesymbol", referenceGene = background, 
                    referenceGeneType = "genesymbol", 
                    sigMethod = "fdr", fdrThr = 0.05,
                    isOutput = FALSE)   

## Loading the functional categories...
## Loading the ID list...
## Loading the reference list...
## Performing the enrichment analysis...

## Warning in oraEnrichment(interestGeneList, referenceGeneList, geneSet, minNum =
## minNum, : No significant gene set is identified based on FDR 0.05!