11 PCB ~ RNA-Cyto
11.2 Data
No differentially expressed genes related to cytokines.
rna_meta <-
read.csv("../hg-rna/Input Files/metadata.csv") %>%
mutate(sample=gsub("Hg","", sample),
sample=gsub("a", "", sample),
sample=gsub("b", "", sample),
sample=gsub("11", "1", sample),
sample=as.factor(sample))
pcb <-
read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>%
rename("sample" = "X") %>%
mutate(pcbbin = ifelse(pcbbin=="1", "PCB+", "PCB-"),
sample = as.factor(sample))
cyto <-
read.csv("../hg-cyto/Output Files/cleaned_cytokine.csv") %>%
rename("sample" = "Sample") %>%
mutate(sample=as.factor(sample))
cytobin <-
read.csv("../hg-cyto/Output Files/cleaned_cytokine_bin.csv") %>%
rename("sample" = "Sample") %>%
mutate(sample=as.factor(sample))
samples <-
intersect(intersect(cyto$sample, pcb$sample), rna_meta$sample)
rna <-
read.csv("../hg-rna/Output Files/txome_genecounts_locid_nohemo.csv", row.names = 1) %>%
setNames(gsub("Hg", "", names(.))) %>%
setNames(gsub("a", "", names(.))) %>%
setNames(gsub("b", "", names(.))) %>%
setNames(gsub("11", "1", names(.))) %>%
select(all_of(samples))
pcb_cyto <-
pcb %>%
filter(sample %in% samples) %>%
merge(., cyto, by="sample")
pcb_cytobin <-
pcb %>%
filter(sample %in% samples) %>%
merge(., cytobin, by="sample")11.2.1 Filter & Normalize RNAseq
Normalize based on library size, includes all samples - males and females.
# Order RNAseq by pcb_cyto
rna <-
rna %>%
select(pcb_cyto$sample)
# Create DGE
dge <- DGEList(counts = rna, group = pcb_cyto$pcbbin)
dge$samples## group lib.size norm.factors
## 1256 PCB- 4573826 1
## 1261 PCB- 1943267 1
## 1262 PCB- 4653972 1
## 1265 PCB+ 2979315 1
## 1272 PCB- 4708624 1
## 1273 PCB- 5011031 1
## 1274 PCB- 3026478 1
## 1282 PCB+ 2951861 1
## 1283 PCB- 4658064 1
## 1288 PCB- 2823770 1
## 1351 PCB- 1430600 1
## 1354 PCB- 3618804 1
## 1380 PCB- 3204059 1
## 1402 PCB- 6560784 1
## 1403 PCB- 6100522 1
## 1404 PCB- 3533565 1
## 1407 PCB- 4238107 1
## 1513 PCB+ 2842670 1
## 1516 PCB- 4688415 1
## 1527 PCB- 3763080 1
## 1528 PCB- 2894695 1
## 861 PCB- 5407908 1
## 862 PCB- 4075756 1# Calculate normalization factors
dge <- calcNormFactors(dge)
# Normalize
cpm <-
cpm(dge, normalized.lib.sizes = TRUE, log = TRUE)
# Visuzalize to check
boxplot(cpm, las = 2)
11.2.2 Cytokine genes normalized
cyto_genes_recep <-
cpm %>%
data.frame() %>%
setNames(gsub("X", "", names(.))) %>%
rownames_to_column(var = "gene") %>%
filter(grepl("IFNG", gene) |
grepl("^IL2[A-Z]", gene) |
grepl("^IL2$", gene) |
grepl("^IL6", gene) |
grepl("^IL7", gene) |
grepl("CXCL8", gene) |
grepl("IL15", gene) |
grepl("IL10", gene) |
grepl("IL18", gene) |
grepl("CXCL1$", gene))
cyto_genes <-
cyto_genes_recep %>%
filter(gene=="IFNG" |
gene=="IL2" |
gene=="IL6" |
gene=="IL7" |
gene=="CXCL8" |
gene=="IL15" |
gene=="IL10" |
gene=="IL18" |
gene=="CXCL1") %>%
pivot_longer(!gene, names_to="sample", values_to="count") %>%
merge(., pcb_cyto[,c("sample", "pcbbin")], by="sample")
cyto_genes_raw <-
rna %>%
select(cyto_genes$sample) %>%
rownames_to_column(var = "gene") %>%
filter(gene %in% cyto_genes$gene) %>%
mutate(gene = gsub("IFNG", "IFNg", gene),
gene = gsub("IL2", "IL.2", gene),
gene = gsub("IL6", "IL.6", gene),
gene = gsub("IL7", "IL.7", gene),
gene = gsub("CXCL8", "IL.8", gene),
gene = gsub("IL15", "IL.15", gene),
gene = gsub("IL10", "IL.10", gene),
gene = gsub("IL18", "IL.18", gene),
gene = gsub("CXCL1", "KC.like", gene))
cyto_genes <-
cyto_genes %>%
mutate(gene = gsub("IFNG", "IFNg", gene),
gene = gsub("IL2", "IL.2", gene),
gene = gsub("IL6", "IL.6", gene),
gene = gsub("IL7", "IL.7", gene),
gene = gsub("CXCL8", "IL.8", gene),
gene = gsub("IL15", "IL.15", gene),
gene = gsub("IL10", "IL.10", gene),
gene = gsub("IL18", "IL.18", gene),
gene = gsub("CXCL1", "KC.like", gene))11.2.3 Cytokine genes raw
cyto_genes_rawbin <-
cyto_genes_raw %>%
mutate(across(!gene, \(x) ifelse(x == "0", "0", "1"))) %>%
column_to_rownames(var = "gene") %>%
t() %>%
data.frame() %>%
rownames_to_column(var = "sample") %>%
merge(., pcb, by = "sample")
cyto_genes_raw <-
cyto_genes_raw %>%
column_to_rownames(var = "gene") %>%
t() %>%
data.frame() %>%
rownames_to_column(var = "sample") %>%
merge(., pcb, by = "sample")11.3 Cytokine stats
cyto_stat <-
apply(pcb_cyto[,25:33], 2, function(x) sum(x>0)) %>%
data.frame() %>%
mutate(perc=round((./23)*100, digits=2))
cyto_stat## . perc
## IFNg 12 52.17
## IL.2 8 34.78
## IL.6 2 8.70
## IL.7 15 65.22
## IL.8 0 0.00
## IL.15 1 4.35
## KC.like 6 26.09
## IL.10 4 17.39
## IL.18 23 100.0011.4 Gene Stats
gene_stat <-
apply(cyto_genes_rawbin[,2:10], 2, function(x) sum(x>0)) %>%
data.frame() %>%
mutate(perc=round((./23)*100, digits=2))
gene_stat## . perc
## KC.like 23 100.00
## IL.8 22 95.65
## IFNg 3 13.04
## IL.10 19 82.61
## IL.15 22 95.65
## IL.18 23 100.00
## IL.2 12 52.17
## IL.6 7 30.43
## IL.7 22 95.6511.4.1 Remove cytokines detected in few samples
remove <-
cyto_stat %>%
rownames_to_column(var="cytokine") %>%
filter(.<3)
pcb_cyto <-
pcb_cyto %>%
select(-remove$cytokine)
pcb_cytobin <-
pcb_cytobin %>%
select(-remove$cytokine)
cyto_genes <-
cyto_genes %>%
filter(!gene %in% remove$cytokine)
cyto_genes_raw <-
cyto_genes_raw %>%
select(-remove$cytokine)
cyto_genes_rawbin <-
cyto_genes_rawbin %>%
select(-remove$cytokine)
cytokine_order <-
c("IFNg", "IL.2", "IL.7", "KC.like", "IL.10", "IL.18")11.5 Individual Cytokines ~ PCBs
11.5.1 IFNg
11.5.1.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(pcb_cytobin$IFNg, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$IFNg, pcb_cytobin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(pcb_cytobin$IFNg, pcb_cytobin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.08754753 129.48831785
## sample estimates:
## odds ratio
## 1.94222611.5.1.2 Concentration
Not sig
data_ifng <-
pcb_cyto %>%
select(IFNg, pcbbin) %>%
filter(!IFNg == 0)
anova_ifng <- aov(log(IFNg) ~ pcbbin, data = data_ifng)
plot(anova_ifng)
##
## Shapiro-Wilk normality test
##
## data: anova_ifng$residuals
## W = 0.89369, p-value = 0.1315## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 0.006 0.0058 0.004 0.952
## Residuals 10 15.227 1.522711.5.2 IL-2
11.5.2.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(pcb_cytobin$IL.2, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$IL.2, pcb_cytobin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(pcb_cytobin$IL.2, pcb_cytobin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.01380025 21.03362878
## sample estimates:
## odds ratio
## 0.931533511.5.2.2 Concentration
Not sig
data_il2 <-
pcb_cyto %>%
select(IL.2, pcbbin) %>%
filter(!IL.2 == 0)
anova_il2 <- aov(IL.2 ~ pcbbin, data = data_il2)
plot(anova_il2)## Warning: not plotting observations with leverage one:
## 6
##
## Shapiro-Wilk normality test
##
## data: anova_il2$residuals
## W = 0.95404, p-value = 0.7518## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 257 257.4 0.459 0.523
## Residuals 6 3366 561.111.5.3 IL-7
11.5.3.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(pcb_cytobin$IL.7, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$IL.7, pcb_cytobin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(pcb_cytobin$IL.7, pcb_cytobin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.04754291 72.46247693
## sample estimates:
## odds ratio
## 1.07349911.5.3.2 Concentration
Not sig
data_il7 <-
pcb_cyto %>%
select(IL.7, pcbbin) %>%
filter(!IL.7 == 0)
anova_il7 <- aov(log(IL.7) ~ pcbbin, data = data_il7)
plot(anova_il7)
##
## Shapiro-Wilk normality test
##
## data: anova_il7$residuals
## W = 0.94932, p-value = 0.5137## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 0.81 0.8124 0.305 0.59
## Residuals 13 34.66 2.666311.5.4 KC-like
11.5.4.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(pcb_cytobin$KC.like, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$KC.like, pcb_cytobin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(pcb_cytobin$KC.like, pcb_cytobin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.02122636 34.49703185
## sample estimates:
## odds ratio
## 1.47189111.5.4.2 Concentration
Not sig
data_kc <-
pcb_cyto %>%
select(KC.like, pcbbin) %>%
filter(!KC.like == 0)
anova_kc <- aov(KC.like ~ pcbbin, data = data_kc)
plot(anova_kc)## Warning: not plotting observations with leverage one:
## 3
##
## Shapiro-Wilk normality test
##
## data: anova_kc$residuals
## W = 0.82152, p-value = 0.09096## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 0.0822 0.08216 0.378 0.572
## Residuals 4 0.8695 0.2173811.5.5 IL-10
11.5.5.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(pcb_cytobin$IL.10, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$IL.10, pcb_cytobin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(pcb_cytobin$IL.10, pcb_cytobin$pcbbin)
## p-value = 0.4529
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.03610715 69.42659296
## sample estimates:
## odds ratio
## 2.67506111.5.5.2 Concentration
Not sig
data_il10 <-
pcb_cyto %>%
select(IL.10, pcbbin) %>%
filter(!IL.10 == 0)
anova_il10 <- aov(IL.10 ~ pcbbin, data = data_il10)
plot(anova_il10)## Warning: not plotting observations with leverage one:
## 2
##
## Shapiro-Wilk normality test
##
## data: anova_il10$residuals
## W = 0.82743, p-value = 0.1612## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 15.8 15.8 0.25 0.667
## Residuals 2 126.4 63.211.5.6 IL-18
11.5.6.1 Presence/Absence
Doesn’t make sense to do, detected in all samples?
tryCatch({
chisq.test(table(pcb_cytobin$IL.18, pcb_cytobin$pcbbin))
}, warning = function(w) {
fisher.test(table(pcb_cytobin$IL.18, pcb_cytobin$pcbbin))
})##
## Chi-squared test for given probabilities
##
## data: table(pcb_cytobin$IL.18, pcb_cytobin$pcbbin)
## X-squared = 12.565, df = 1, p-value = 0.00039311.5.6.2 Concentration
Not sig
data_il18 <-
pcb_cyto %>%
select(IL.18, pcbbin) %>%
filter(!IL.18 == 0)
anova_il18 <- aov(IL.18 ~ pcbbin, data = data_il18)
plot(anova_il18)
##
## Shapiro-Wilk normality test
##
## data: anova_il18$residuals
## W = 0.95199, p-value = 0.3216## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 196 196.2 0.508 0.484
## Residuals 21 8110 386.211.6 Individual Cytokine Genes ~ PCBs
11.6.1 IFNg
11.6.1.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(cyto_genes_rawbin$IFNg, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$IFNg, cyto_genes_rawbin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(cyto_genes_rawbin$IFNg, cyto_genes_rawbin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.00000 20.26128
## sample estimates:
## odds ratio
## 011.6.2 IL-2
11.6.2.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(cyto_genes_rawbin$IL.2, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$IL.2, cyto_genes_rawbin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(cyto_genes_rawbin$IL.2, cyto_genes_rawbin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.08754753 129.48831785
## sample estimates:
## odds ratio
## 1.94222611.6.2.2 Concentration
Not sig
rnadata_il2 <-
cyto_genes_raw %>%
select(IL.2, pcbbin) %>%
filter(!IL.2 == 0)
anova_il2 <- aov(IL.2 ~ pcbbin, data = data_il2)
plot(anova_il2)## Warning: not plotting observations with leverage one:
## 6
##
## Shapiro-Wilk normality test
##
## data: anova_il2$residuals
## W = 0.95404, p-value = 0.7518## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 257 257.4 0.459 0.523
## Residuals 6 3366 561.111.6.3 IL-7
11.6.3.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(cyto_genes_rawbin$IL.7, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$IL.7, cyto_genes_rawbin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(cyto_genes_rawbin$IL.7, cyto_genes_rawbin$pcbbin)
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.003861719 Inf
## sample estimates:
## odds ratio
## Inf11.6.3.2 Concentration
Not sig
rnadata_il7 <-
cyto_genes_raw %>%
select(IL.7, pcbbin) %>%
filter(!IL.7 == 0)
anova_il7 <- aov(log(IL.7) ~ pcbbin, data = data_il7)
plot(anova_il7)
##
## Shapiro-Wilk normality test
##
## data: anova_il7$residuals
## W = 0.94932, p-value = 0.5137## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 0.81 0.8124 0.305 0.59
## Residuals 13 34.66 2.666311.6.4 KC-like
11.6.4.1 Presence/Absence
**Sig
tryCatch({
chisq.test(table(cyto_genes_rawbin$KC.like, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$KC.like, cyto_genes_rawbin$pcbbin))
})##
## Chi-squared test for given probabilities
##
## data: table(cyto_genes_rawbin$KC.like, cyto_genes_rawbin$pcbbin)
## X-squared = 12.565, df = 1, p-value = 0.00039311.6.4.2 Concentration
Not sig
rnadata_kc <-
cyto_genes_raw %>%
select(KC.like, pcbbin) %>%
filter(!KC.like == 0)
anova_kc <- aov(KC.like ~ pcbbin, data = data_kc)
plot(anova_kc)## Warning: not plotting observations with leverage one:
## 3
##
## Shapiro-Wilk normality test
##
## data: anova_kc$residuals
## W = 0.82152, p-value = 0.09096## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 0.0822 0.08216 0.378 0.572
## Residuals 4 0.8695 0.2173811.6.5 IL-10
11.6.5.1 Presence/Absence
Not sig
tryCatch({
chisq.test(table(cyto_genes_rawbin$IL.10, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$IL.10, cyto_genes_rawbin$pcbbin))
})##
## Fisher's Exact Test for Count Data
##
## data: table(cyto_genes_rawbin$IL.10, cyto_genes_rawbin$pcbbin)
## p-value = 0.4529
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.0144037 27.6953488
## sample estimates:
## odds ratio
## 0.373823311.6.5.2 Concentration
Not sig
rnadata_il10 <-
cyto_genes_raw %>%
select(IL.10, pcbbin) %>%
filter(!IL.10 == 0)
anova_il10 <- aov(IL.10 ~ pcbbin, data = data_il10)
plot(anova_il10)## Warning: not plotting observations with leverage one:
## 2
##
## Shapiro-Wilk normality test
##
## data: anova_il10$residuals
## W = 0.82743, p-value = 0.1612## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 15.8 15.8 0.25 0.667
## Residuals 2 126.4 63.211.6.6 IL-18
11.6.6.1 Presence/Absence
**Sig Doesn’t make sense to do, detected in all samples?
tryCatch({
chisq.test(table(cyto_genes_rawbin$IL.18, cyto_genes_rawbin$pcbbin))
}, warning = function(w) {
fisher.test(table(cyto_genes_rawbin$IL.18, cyto_genes_rawbin$pcbbin))
})##
## Chi-squared test for given probabilities
##
## data: table(cyto_genes_rawbin$IL.18, cyto_genes_rawbin$pcbbin)
## X-squared = 12.565, df = 1, p-value = 0.00039311.6.6.2 Concentration
Not sig
rnadata_il18 <-
cyto_genes_raw %>%
select(IL.18, pcbbin) %>%
filter(!IL.18 == 0)
anova_il18 <- aov(IL.18 ~ pcbbin, data = data_il18)
plot(anova_il18)
##
## Shapiro-Wilk normality test
##
## data: anova_il18$residuals
## W = 0.95199, p-value = 0.3216## Df Sum Sq Mean Sq F value Pr(>F)
## pcbbin 1 196 196.2 0.508 0.484
## Residuals 21 8110 386.211.7 Plot Cytokines & Genes
11.7.1 Presence/Absence
11.7.1.1 Cytokines
pcb_cytobin %>%
select(sample, pcbbin, IFNg, IL.2, IL.7, KC.like, IL.10, IL.18) %>%
pivot_longer(cols = -c(sample, pcbbin), names_to = "cytokine", values_to = "presence") %>%
mutate(presence = as.numeric(presence)) %>%
group_by(pcbbin, cytokine) %>%
summarize(Present =sum(presence)/n(),
Absent = (1-Present)) %>%
pivot_longer(cols = -c(cytokine, pcbbin),
names_to = "presence", values_to = "prop") %>%
ggplot(., aes(pcbbin, y = prop, fill = presence)) +
geom_col(position = "fill") +
facet_wrap(~ cytokine, ncol = 3) +
scale_fill_manual(values = c("gray60", "slategray3"), name = "Cytokine") +
labs(x = "PCB Status", y = "Proportion of Samples") +
theme_bw() +
theme(panel.grid = element_blank())
11.7.1.2 Genes
cyto_genes_rawbin %>%
select(sample, pcbbin, IFNg, IL.2, IL.7, KC.like, IL.10, IL.18) %>%
pivot_longer(cols = -c(sample, pcbbin), names_to = "cytokine", values_to = "presence") %>%
mutate(presence = as.numeric(presence)) %>%
group_by(pcbbin, cytokine) %>%
summarize(Present = sum(presence)/n(),
Absent = (1-Present)) %>%
pivot_longer(cols = -c(cytokine, pcbbin),
names_to = "presence", values_to = "prop") %>%
ggplot(., aes(pcbbin, y = prop, fill = presence)) +
geom_col(position = "fill") +
facet_wrap(~ cytokine, ncol = 3) +
scale_fill_manual(values = c("gray60", "slategray3"), name = "Gene") +
labs(x = "PCB Status", y = "Proportion of Samples") +
theme_bw() +
theme(panel.grid = element_blank())
11.7.2 Concentration
11.7.2.1 Cytokines
pcb_cyto %>%
select(sample, pcbbin, IFNg, IL.2, IL.7, KC.like, IL.10, IL.18) %>%
mutate_at(vars(c(3:8)), scale) %>%
pivot_longer(-c(sample, pcbbin), names_to = "cytokine", values_to = "count") %>%
mutate(cytokine = factor(cytokine, levels = cytokine_order)) %>%
ggplot(., aes(pcbbin, y = count)) +
geom_boxplot() +
facet_wrap(~ cytokine, ncol = 3, scales = "free") +
labs(x = "PCB Status", y = "Normalized Cytokine Concentration") +
theme_bw() +
theme(panel.grid = element_blank())
