## Background


Today we will perform an RNASeq analysis of the effects of a common steroid on airway
cells.

In particular, dexamethasone (hereafter just called “dex”) on different airway smooyh muscle
cell lines (ASM cells).

## Data Import

We need two different inputs:

• countData: with genes in rows and experiments in columns
• colData: meta data that describes the columns in countData

```{r}
counts <- read.csv("airway_scaledcounts.csv", row.names=1)
metadata <- read.csv("airway_metadata.csv")
head(counts)
```

```{r}
metadata
```

Q1. How many genes are in this dataset?

```{r}
nrow(counts)
```

Q2. How many ‘control’ cell lines do we have?

```{r}
table(metadata$dex)
```

## Differential gene expression

We have 4 replicate drug treated and control (no drug) clumns?experments in our ‘counts’
object.
We want one “mean” value for each gene (rows) in “treated” (drug) and one mean value for
each gene in “control” cols.
Step 1. find all “control” colums Step 2. Extract these colums to a new object called ‘con-
trol.counts’ Step 3. Then caculate the mean value fo each gene

```{r}
contol.inds <- metadata$dex=="control"
```

Step 2.
```{r}
control.counts <- counts[ ,contol.inds]
```

Step 3.

```{r}
control.mean <-rowMeans(control.counts)
```


```{r}
metadata
```

Now do the same thing for the “treated” colums/experiments

```{r}
treated.mean <- rowMeans(counts[ ,metadata$dex=="treated"])
meancounts<- data.frame(control.mean, treated.mean)
plot(meancounts)
```

Let’s log transform this count data:

```{r}
plot(meancounts, log="xy")
```


N.B. We most often use log2 for this type of data as it makes the interpretation much more straightforward.
Treated/Control is often called “fold-change”.
if there was no change we would have a log2-fc of zero:
```{r}
log2(10/10)
```

if we had double the amount of trancript around we would have a log2-fc of 1:

```{r}
log2(20/10)
```

if we had half as much trancript around we would have a log2-fc of -1

```{r}
log2(5/10)
```

Q. Caculate a log2 fold change value for all our genes and add it as a new column to our ‘meanscounts’ object.

```{r}
meancounts$log2f<- log2 (meancounts$treated.mean /
meancounts$control.mean)
head(meancounts)
```

There are some “funky” values that come about when ever we have 0 mean count values. Typically remove these genes from any further analysis- as we can’t say anything about them
if we have no data for them.

```{r}
log2(40/10)
```

## DESeq anlaysis

let’s do this anlaysis with an estimate of statiscal significance using the DESeq2 package.


```{r}
library(DESeq2)
```
DESeq (like bioconductor packages) want it’s input data in a very specific way.

```{r}
dds <-DESeqDataSetFromMatrix(countData= counts,
colData = metadata,
design = ~dex)
```

## Run the DESeq anlaysis pipline
The main function ‘DESeq()’

```{r}
dds <- DESeq(dds)
```

```{r}
res<- results(dds)
head(res)
```
```{r}
36000*0.05
```


## Volcano plot
This is a main summary results figure from these kinds of studies. it is a plot of log2-foldchnage
vs of (Adjusted) P-value.

```{r}
plot(res$log2FoldChange,
res$padj)
```

Again this y-axis is a highly needs log transforming and we can flip the y-axis with a minus sign so it looks like every other volcano plot.

```{r}
plot(res$log2FoldChange,
-log(res$padj))
abline(v=-2, col="red")
abline(v=+2, col="red")
abline(h=-log(0.05), col="red")
```

### Adding some colors annotation
Start with default base color “gray”

```{r}
# Start with a default base color "gray"
mycols <- rep("gray", nrow(res))
mycols[ res$log2FoldChange > 2 ] <- "blue"
mycols[ res$log2FoldChange < -2 ] <- "darkgreen"
mycols[res$padj >=0.05]<-"pink"
plot(res$log2FoldChange,
-log(res$padj),
col = mycols)
#Volcano
plot(res$log2FoldChange,
-log(res$padj),
col=mycols)
abline(v=c(-2,+2), lty=2)
abline(h=-log(0.05), lty=2)
```

Q.Make a presentation quailty ggplot version quality of this plot. Include clear axis labels, a clean theme , your custom colors, cut-off lines and a plot title.

```{r}
library(ggplot2)
```


```{r}
df <- as.data.frame(res)
df <- df[!is.na(df$log2FoldChange) & !is.na(df$padj), ]
# same color logic as your base plot
df$color <- "gray"
df$color[df$log2FoldChange > 2] <- "blue"
df$color[df$log2FoldChange < -2] <- "darkgreen"
df$color[df$padj >= 0.05] <- "pink"

ggplot(df, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = color), size = 1.5) +
  geom_vline(xintercept = c(-2, 2), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  scale_color_identity() +
  labs(
    title = "Volcano Plot",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-value"
  ) +
  theme_classic()

```

## Save our results
Write a CVS file

```{r}
write.csv(res, file="results.csv")
```












