Section 1: Set-up

#loading all necessary packages
library(ggplot2)
library(tidyr)
library(biotools)

## Loading required package: MASS

## ---
## biotools version 4.2

library(ltm)

## Loading required package: msm

## Loading required package: polycor

Section 2: Simple Linear Regression

logFC <- read.csv("Shi_logFC_subset.csv")

BaP.plot <- ggplot(data=logFC, aes(x=combination, y=BaP))+
geom_point(shape = 1) +
geom_jitter(data= logFC, aes(x=combination, y=BaP),shape = 1, width =0.2)

LPS.plot <- ggplot(data=logFC, aes(x=combination, y=LPS))+
geom_point(shape = 1) +
geom_jitter(data= logFC, aes(x=combination, y=LPS),shape = 1, width =0.2)

BaP.plot

LPS.plot

Question 1

Answer It appears that the combination-LPS plot is positively correlated, and combination-BaP plot does not or has a poor positive correlation. Hence, the correlation of combination-LPS is stronger of the two.

res <- cor.test(~ combination + BaP, data = logFC)
res

## 
##  Pearson's product-moment correlation
## 
## data:  combination and BaP
## t = 16.453, df = 998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4117152 0.5093266
## sample estimates:
##       cor 
## 0.4619185

lm.BaP <- lm(combination ~ BaP, data = logFC)
summary(lm.BaP)

## 
## Call:
## lm(formula = combination ~ BaP, data = logFC)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.01032 -0.14972 -0.00823  0.12756  3.04406 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.009121   0.012216  -0.747    0.455    
## BaP          0.585712   0.035599  16.453   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.383 on 998 degrees of freedom
## Multiple R-squared:  0.2134, Adjusted R-squared:  0.2126 
## F-statistic: 270.7 on 1 and 998 DF,  p-value: < 2.2e-16

Question 2

Answer It appears that combination-BaP has a weak to moderate positive correlation.

Question 3

Answer If the BaP value were high for a gene, then it is somewhat likely that the combination value would be higher. With the caveat that Pearson’c correlation can’t differentiate independent vs dependent variables and the r-squared value is low, it cannot be predicted with certainty.

#‘geom_smooth()‘ using formula ’y ~ x’
BaP.plot +
geom_smooth( method = lm )

## `geom_smooth()` using formula 'y ~ x'

## ‘geom_smooth()‘ using formula ’y ~ x’
LPS.plot +
geom_smooth( method = lm )

## `geom_smooth()` using formula 'y ~ x'

Section 3: Multiple Linear Regression

Question 4

Answer

Null Hypothesis: Changes in gene expression due to combination of (LPS and BaP) are NOT DIFFERENT from the changes due to BaP treatment or LPS treatment.

OR:

All correlation coefficients in the model are zero.

Alternative Hypothesis: Changes in gene expression due to combination of (LPS and BaP) DIFFERS from either the changes due to BaP treatment or LPS treatment or both.

OR:

Not all correlation coefficients in the model are simultaneously zero.

lm.both <- lm( combination ~ BaP*LPS, data = logFC)
summary(lm.both)

## 
## Call:
## lm(formula = combination ~ BaP * LPS, data = logFC)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.47655 -0.10204 -0.00242  0.09883  1.10698 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.022277   0.006911  -3.223  0.00131 ** 
## BaP          0.208022   0.025527   8.149 1.09e-15 ***
## LPS          0.889906   0.019299  46.112  < 2e-16 ***
## BaP:LPS     -0.017127   0.027740  -0.617  0.53710    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2153 on 996 degrees of freedom
## Multiple R-squared:  0.7519, Adjusted R-squared:  0.7511 
## F-statistic:  1006 on 3 and 996 DF,  p-value: < 2.2e-16

#strength of the correlation

sqrt(summary(lm.both)$r.squared)

## [1] 0.8671021

Question 5

Answer In the simple linear regression model, based on the F statistic, significant p-value, and the R-squared value, we infer that that only 21.26% of the change in gene expression due to the combination of BaP and LPS can be explained by BaP exposure alone. Since the R-squared value is low, this is not a good model, i.e., there are other variables at play.

Thus, we need a model that accounts for the effect of LPS, too, on gene expression, and hence, we do a multiple linear regression analysis.

We find that this is a decent/good model to describe the change in gene expression (based on the R-squared value of 0.75). We infer that 75.11% of the changes in gene expression due to the combination of BaP and LPS can be explained by BaP exposure, and LPS ttreatment.

It is interesting to note that in the multiple regression analysis, the interecept value is significant.

The intercept value reflects the gene expression changes due combination of BaP and LPS when affects due to BaP and affects due to LPS, both, are zero. Negative intercept value implies downregulation.

Question 6

Answer

Multiple linear regression was carried out to investigate whether the change in gene expression due to the combination of LPS and BaP can be explained by the changes in gene expression due to the individual treatment with BaP alone or LPS alone.

There was a significant relationship between combination and BaP (p = 1.09e-15), combination and LPS (p < 2e-16).

We predict that because of combination, the change in gene expression is 1.155 (i.e., upregulation) times the change due to BaP alone. The change in gene expression due to combination is 1.853 (i.e., upregulation) times the change due to LPS alone. Further, there is a 0.9841 times downregulation in gene expression due to combination, in the scenario where neither BaP exposure nor LPS treatment occurs (intercept value).

The adjusted R-squared value is 0.7511 so 75% of the changes in gene expression due to combination can be explained by BaP treatment and LPS treatment. The data met the assumptions of homogeneity of variance and linearity and the residuals were approximately normally distributed.

FOLD-CHANGE VALUE CALCULATIONS:

From the multiple linear regression analysis, the estimates (correlation coefficients) are the mean log2(fold-change) gene expression of the replicates for that group, relative to the control.

Therefore, to get the fold-change, we do an anti-log as follows-

Antilog2 (-0.022277) = 0.98467736495 - downregulation

Antilog2 (0.208022) = 1.1551033989 - upregulation

Antilog2 (0.889906) = 1.8530553825 - upregulation

Question 7

Answer

Both pulmonary inflammation (induced by LPS treatment) and exposure of polycyclic aromatic hydrocarbons like BaP can cause upregulation of gene expression. Whereas, there was a downregulation in gene expression in mice with pulmonary inflammation, when they were exposed to BaP.

Section 4: MANOVA

manova.genes <- read.csv("Shi_Cyp1a2_Cyp1b1.csv")

ggplot( data = manova.genes, aes( x = condition )) +
geom_jitter( aes(y = Cyp1b1, color="Cyp1b1"), width=0.15, height=0, shape=1 ) +
geom_jitter( aes(y = Cyp1a2, color="Cyp1a2"), width=0.15, height=0, shape=1 ) +
ylab("Expression")

ggplot( data = manova.genes, aes( x = condition )) +
geom_jitter( aes(y = Cyp1b1.transform, color="Cyp1b1.transform"), width=0.15, height=0, shape=1 ) +
geom_jitter( aes(y = Cyp1a2.transform, color="Cyp1a2.transform"), width=0.15, height=0, shape=1 ) +
ylab("Expression")

Question 8

Answer

Null Hypothesis: There is no difference in the [multivariate] mean expression of Cyp1a2 and Cyp1b1 genes in groups treated with BaP, combination, control, or LPS.

Alternative Hypothesis: There is a difference in the [multivariate] mean expression of Cyp1a2 and Cyp1b1 genes, under AT LEAST 2 of the conditions.

Question 9

Answer Based on the plot, the spread in each condition overlaps quite a bit, which prompts me to think that we fail to reject the null hypothesis.

# Box's M test
boxM( manova.genes[ , c("Cyp1a2.transform", "Cyp1b1.transform") ], manova.genes$condition )

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  manova.genes[, c("Cyp1a2.transform", "Cyp1b1.transform")]
## Chi-Sq (approx.) = 13.741, df = 9, p-value = 0.1318

# Cronbach alpha
cronbach.alpha( manova.genes[ , c("Cyp1a2.transform", "Cyp1b1.transform") ] )

## 
## Cronbach's alpha for the 'manova.genes[, c("Cyp1a2.transform", "Cyp1b1.transform")]' data-set
## 
## Items: 2
## Sample units: 15
## alpha: 0.69

# correlation
cor.test(manova.genes$Cyp1a2.transform, manova.genes$Cyp1b1.transform)

## 
##  Pearson's product-moment correlation
## 
## data:  manova.genes$Cyp1a2.transform and manova.genes$Cyp1b1.transform
## t = 2.2646, df = 13, p-value = 0.04128
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02697099 0.82057111
## sample estimates:
##       cor 
## 0.5318852

Question 10

Answer based on the Box’s M-test, the p-value is 0.1318, so we fail to reject the null hypothesis, which means the data passes this test!

ggplot(data = manova.genes, aes(x=Cyp1a2.transform, y=Cyp1b1.transform))+
  geom_jitter( aes(x = Cyp1a2.transform))+
  geom_jitter( aes(y = Cyp1b1.transform))

# univariate tests for each condition group
shapiro.test(manova.genes[1:3,6])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[1:3, 6]
## W = 0.86685, p-value = 0.2866

shapiro.test(manova.genes[1:3,7])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[1:3, 7]
## W = 0.79195, p-value = 0.09541

shapiro.test(manova.genes[4:7,6])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[4:7, 6]
## W = 0.91846, p-value = 0.5284

shapiro.test(manova.genes[4:7,7])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[4:7, 7]
## W = 0.88955, p-value = 0.381

shapiro.test(manova.genes[8:11,6])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[8:11, 6]
## W = 0.88215, p-value = 0.3479

shapiro.test(manova.genes[8:11,7])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[8:11, 7]
## W = 0.92171, p-value = 0.5466

shapiro.test(manova.genes[12:15,6])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[12:15, 6]
## W = 0.9155, p-value = 0.512

shapiro.test(manova.genes[12:15,7])

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes[12:15, 7]
## W = 0.98754, p-value = 0.9446

# tests for all conditions within each gene
shapiro.test(manova.genes$Cyp1a2.transform)

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes$Cyp1a2.transform
## W = 0.90773, p-value = 0.125

shapiro.test(manova.genes$Cyp1b1.transform)

## 
##  Shapiro-Wilk normality test
## 
## data:  manova.genes$Cyp1b1.transform
## W = 0.87706, p-value = 0.04289

manova.results <- manova( cbind(Cyp1a2, Cyp1b1) ~ condition, data = manova.genes)
summary(manova.results)

##           Df  Pillai approx F num Df den Df Pr(>F)
## condition  3 0.19945  0.40615      6     22 0.8669
## Residuals 11

Question 11

Answer

Assuming an alpha value of 0.05:

Since the p-value is 0.8669, we fail to reject the null hypothesis.

This indicates that there is no difference in the impact of the given treatments (control, combination, LPS, BaP) treatments on Cyp1a2 and Cyp1b1 gene expression.

Week7

Sai Krishna Bhamidipati (PID: A59010579)

2022-05-13

Section 1: Set-up

Section 2: Simple Linear Regression

Question 1

Answer It appears that the combination-LPS plot is positively correlated, and combination-BaP plot does not or has a poor positive correlation. Hence, the correlation of combination-LPS is stronger of the two.

Question 2

Answer It appears that combination-BaP has a weak to moderate positive correlation.

Question 3

Answer If the BaP value were high for a gene, then it is somewhat likely that the combination value would be higher. With the caveat that Pearson’c correlation can’t differentiate independent vs dependent variables and the r-squared value is low, it cannot be predicted with certainty.

Section 3: Multiple Linear Regression

Question 4

Answer

Null Hypothesis: Changes in gene expression due to combination of (LPS and BaP) are NOT DIFFERENT from the changes due to BaP treatment or LPS treatment.

OR:

All correlation coefficients in the model are zero.

Alternative Hypothesis: Changes in gene expression due to combination of (LPS and BaP) DIFFERS from either the changes due to BaP treatment or LPS treatment or both.

OR:

Not all correlation coefficients in the model are simultaneously zero.

Question 5

Thus, we need a model that accounts for the effect of LPS, too, on gene expression, and hence, we do a multiple linear regression analysis.

We find that this is a decent/good model to describe the change in gene expression (based on the R-squared value of 0.75). We infer that 75.11% of the changes in gene expression due to the combination of BaP and LPS can be explained by BaP exposure, and LPS ttreatment.

It is interesting to note that in the multiple regression analysis, the interecept value is significant.

The intercept value reflects the gene expression changes due combination of BaP and LPS when affects due to BaP and affects due to LPS, both, are zero. Negative intercept value implies downregulation.

Question 6

Answer

Multiple linear regression was carried out to investigate whether the change in gene expression due to the combination of LPS and BaP can be explained by the changes in gene expression due to the individual treatment with BaP alone or LPS alone.

There was a significant relationship between combination and BaP (p = 1.09e-15), combination and LPS (p < 2e-16).

The adjusted R-squared value is 0.7511 so 75% of the changes in gene expression due to combination can be explained by BaP treatment and LPS treatment. The data met the assumptions of homogeneity of variance and linearity and the residuals were approximately normally distributed.

FOLD-CHANGE VALUE CALCULATIONS:

From the multiple linear regression analysis, the estimates (correlation coefficients) are the mean log2(fold-change) gene expression of the replicates for that group, relative to the control.

Therefore, to get the fold-change, we do an anti-log as follows-

Antilog2 (-0.022277) = 0.98467736495 - downregulation

Antilog2 (0.208022) = 1.1551033989 - upregulation

Antilog2 (0.889906) = 1.8530553825 - upregulation

Question 7

Answer

Both pulmonary inflammation (induced by LPS treatment) and exposure of polycyclic aromatic hydrocarbons like BaP can cause upregulation of gene expression. Whereas, there was a downregulation in gene expression in mice with pulmonary inflammation, when they were exposed to BaP.

Section 4: MANOVA

Question 8

Answer

Null Hypothesis: There is no difference in the [multivariate] mean expression of Cyp1a2 and Cyp1b1 genes in groups treated with BaP, combination, control, or LPS.

Alternative Hypothesis: There is a difference in the [multivariate] mean expression of Cyp1a2 and Cyp1b1 genes, under AT LEAST 2 of the conditions.

Question 9

Answer Based on the plot, the spread in each condition overlaps quite a bit, which prompts me to think that we fail to reject the null hypothesis.

Question 10

Answer based on the Box’s M-test, the p-value is 0.1318, so we fail to reject the null hypothesis, which means the data passes this test!

Question 11

Answer

Assuming an alpha value of 0.05:

Since the p-value is 0.8669, we fail to reject the null hypothesis.

This indicates that there is no difference in the impact of the given treatments (control, combination, LPS, BaP) treatments on Cyp1a2 and Cyp1b1 gene expression.