RNAseqanalysis.Rmd

---
title: "Run1_analysis"
author: "ERM"
date: "20230118"
output: 
  html_document: 
    toc: yes
    theme: flatly
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Hello! I am first loading all the beautiful libraries I will need.

```{r library loading, warning=FALSE, echo=TRUE, message=FALSE}
library(Biobase)
library(edgeR)
library(limma)
library(RColorBrewer)
library(mixOmics)
library(VennDiagram)
library(HTSFilter)
library(gridExtra)
library(reshape2)
library(devtools)
library(AnnotationHub)
library(tidyverse)
library(scales)
library(biomaRt)
library(Homo.sapiens)
library(cowplot)
library(ggrepel)
library(corrplot)
library(Hmisc)

```

The next step is to load all the data I will be using. Currently, I am
not posting the raw data, but I will release in the future.

```{r ggplot-theme-and-functions, warning=FALSE, include=FALSE}
# theme_set(theme_bw(base_size = 16))


pca_plot <- function(df, col_var = NULL, shape_var = NULL, title = "") {
  ggplot(df) + geom_point(aes_string(x = "PC1", y = "PC2", color = col_var,
                                     shape = shape_var), size = 5) +
    labs(title = title, x = "PC1", y = "PC2")
}

pca_var_plot <- function(pca) {
  # x: class == prcomp
  pca.var <- pca$sdev^2
  pca.prop <- pca.var / sum(pca.var)
  var.plot <- qplot(PC, prop, data = data.frame(PC = 1:length(pca.prop),
                                                prop = pca.prop)) +
    labs(title = 'Variance contributed by each PC',
         x = 'PC', y = 'Proportion of variance')
}

calc_pca <- function(x) {
  # Performs principal components analysis with prcomp
  # x: a sample-by-gene numeric matrix
  prcomp(x, scale. = TRUE, retx = TRUE)
}

get_regr_pval <- function(mod) {
  # Returns the p-value for the Fstatistic of a linear model
  # mod: class lm
  stopifnot(class(mod) == "lm")
  fstat <- summary(mod)$fstatistic
  pval <- 1 - pf(fstat[1], fstat[2], fstat[3])
  return(pval)
}

plot_versus_pc <- function(df, pc_num, fac) {
  # df: data.frame
  # pc_num: numeric, specific PC for plotting
  # fac: column name of df for plotting against PC
  pc_char <- paste0("PC", pc_num)
  # Calculate F-statistic p-value for linear model
  pval <- get_regr_pval(lm(df[, pc_char] ~ df[, fac]))
  if (is.numeric(df[, f])) {
    ggplot(df, aes_string(x = f, y = pc_char)) + geom_point() +
      geom_smooth(method = "lm") + labs(title = sprintf("p-val: %.2f", pval))
  } else {
    ggplot(df, aes_string(x = f, y = pc_char)) + geom_boxplot() +
      labs(title = sprintf("p-val: %.2f", pval))
  }
}

x_axis_labels=function(labels,every_nth=1,...) {
    axis(side=1,at=seq_along(labels),labels=F)
    text(x=(seq_along(labels))[seq_len(every_nth)==1],
        y=par("usr")[3]-0.075*(par("usr")[4]-par("usr")[3]),
        labels=labels[seq_len(every_nth)==1],xpd=TRUE,...)
}
```

```{r file loading, eval=TRUE, include=FALSE}
print(getwd())
design <- read.csv("data/data_outline.txt", row.names = 1)
mymatrix <- readRDS("data/mymatrix.RDS")

```

```{r  current data import,eval=FALSE, echo =FALSE}

myfiles <- list.files(path = "~/Ward Lab/Cardiotoxicity/Data/counts/", pattern ="*.txt")
design <- read.csv("data/data_outline.txt", row.names = 1)

mymatrix <- readDGE(files = myfiles, path = "~/Ward Lab/Cardiotoxicity/Data/counts/", group =as.factor(rep((c("1","2","3","4","5","6","7","8","9","10","11","12")),6)))

```

```{r setting up the files used, echo=FALSE, message=FALSE, warning=FALSE}
#samplenames <- substring(colnames(mymatrix), 0, 15)##this pulls the column names from the matrix and clips char 0-15 to rename)


time <- rep((rep(c("3h", "24h"), c(6,6))), 6) 
time <- ordered(time, levels =c("3h", "24h"))
mymatrix$samples$time <- time
indv <- as.factor(rep(c(1,2,3,4,5,6), c(12,12,12,12,12,12))) 
mymatrix$samples$indv <- indv
group <- as.factor(rep((c("1","2","3","4","5","6","7","8","9","10","11","12")),6))

drug <- rep(c("Daunorubicin","Doxorubicin","Epirubicin","Mitoxantrone","Trastuzumab", "Vehicle"),12)

mymatrix$samples$drug <- drug
RIN <- c(9.3,9.8,  9.8, 10.0,  9.6,  9.9,  7.8,  8.7,  8.6,  9.1,  9.4,  9.3,  9.4,  9.6,  9.9,  9.4,  9.9,  9.8,9.5,  9.9,  9.6,  9.9,  8.8,  9.5,  9.3, 10.0,  9.3,  9.6,  9.7,  9.4,  8.9,  7.5,  8.6,  8.5,  8.6,  9.3,9.6,  9.6,  9.4,  9.7,  9.9, 10.0,  9.4,  9.6,  9.4,  9.5,  9.8,  9.5,  9.6,  9.7,  9.6,  9.9,  9.6, 10.0,8.3, 8.4,9.4,9.7, 10.0, 10.0,  9.6,  9.6,  9.7,  9.5,  9.8,  9.5,  9.8,  8.6,  8.9,  9.1,  9.3,9.6)
mymatrix$samples$RIN <- RIN
samplenames <- c("MCW_RM_R_11_S72", "MCW_RM_R_12_S1_", "MCW_RM_R_13_S2_", "MCW_RM_R_14_S3_", "MCW_RM_R_15_S4_","MCW_RM_R_16_S5_", "MCW_RM_R_17_S6_", "MCW_RM_R_18_S7_", "MCW_RM_R_19_S8_", "MCW_RM_R_20_S9_","MCW_RM_R_21_S10", "MCW_RM_R_22_S11", "MCW_RM_R_23_S12", "MCW_RM_R_24_S13", "MCW_RM_R_25_S14","MCW_RM_R_26_S15", "MCW_RM_R_27_S16", "MCW_RM_R_28_S17", "MCW_RM_R_29_S18", "MCW_RM_R_30_S19","MCW_RM_R_31_S20", "MCW_RM_R_32_S21", "MCW_RM_R_33_S22", "MCW_RM_R_34_S23", "MCW_RM_R_35_S60","MCW_RM_R_36_S61", "MCW_RM_R_37_S62", "MCW_RM_R_38_S63", "MCW_RM_R_39_S64", "MCW_RM_R_40_S65","MCW_RM_R_41_S66", "MCW_RM_R_42_S67", "MCW_RM_R_43_S68", "MCW_RM_R_44_S69","MCW_RM_R_45_S70","MCW_RM_R_46_S71","MCW_RM_R_47_S24","MCW_RM_R_48_S25", "MCW_RM_R_49_S26", "MCW_RM_R_50_S27","MCW_RM_R_51_S28", "MCW_RM_R_52_S29", "MCW_RM_R_53_S30", "MCW_RM_R_54_S31", "MCW_RM_R_55_S32","MCW_RM_R_56_S33","MCW_RM_R_57_S34", "MCW_RM_R_58_S35", "MCW_RM_R_59_S36", "MCW_RM_R_60_S37","MCW_RM_R_61_S38", "MCW_RM_R_62_S39", "MCW_RM_R_63_S40", "MCW_RM_R_64_S41", "MCW_RM_R_65_S42","MCW_RM_R_66_S43", "MCW_RM_R_67_S44", "MCW_RM_R_68_S45", "MCW_RM_R_69_S46", "MCW_RM_R_70_S47","MCW_RM_R_71_S48","MCW_RM_R_72_S49", "MCW_RM_R_73_S50", "MCW_RM_R_74_S51", "MCW_RM_R_75_S52","MCW_RM_R_76_S53", "MCW_RM_R_77_S54", "MCW_RM_R_78_S55", "MCW_RM_R_79_S56", "MCW_RM_R_80_S57","MCW_RM_R_81_S58","MCW_RM_R_82_S59")
 
anno <- (cbind(samplenames, indv, drug, time, RIN, group))
anno <- as.data.frame(anno)
```


This is how I retrieved the gene symbols.

```{r gene id infromation, echo=TRUE, message=FALSE, warning=FALSE}
 
 ###now we add genenames to the geneid###
 
geneid <- rownames(mymatrix) ### pulls the names we have in the counts file
genes <- select(Homo.sapiens, keys=geneid, columns=c("SYMBOL"),
  keytype="ENTREZID")
genes <- genes[!duplicated(genes$ENTREZID),]
mymatrix$genes <- genes

```

## Filtering out low count genes  

Filtering the genes that are lowly expressed using several methods.

First method, removing only those rows with zero counts across all
samples.

```{r first filter method, echo=TRUE}
# 
# old.par <- par(mar = c(0, 0, 0, 0))
# par(old.par)
# boxplot(data =RNAseqreads, total~Sample, main = "Boxplots of total reads",xaxt = "n", xlab= "")
# x_axis_labels(labels = samplenames, every_nth = 1, adj=1, srt =90, cex =0.4)
# ggplot(RNAseqreads, x = Sample, y = total)+
#   geom_boxplot()

table(rowSums(mymatrix$counts==0)==72)
```

This filtering would leave 24931 genes and remove 3464, That is too many
leftover genes!

So now to try something a little more stringent using the built in
function from the edgeR package.

```{r second filter method, echo=TRUE}
keep <- filterByExpr.DGEList(mymatrix, group = group)
filter_test <- mymatrix[keep, , keep.lib.sizes=FALSE]
dim(filter_test)

```

This method effectively uses a cutoff off that leaves 14448 genes. The
cutoff is determined by keeping genes that have a count-per-million
(CPM) above 10, (the default minimum set) in 6 samples. A set is
determined using the design matrix.\
For my design, I grouped my 72 samples into sets of 6, one set includes
each individual + a specific treatment + a specific time.

The beginning cutoff-standard in our lab is to start by using the
rowMeans \>0 cutoff on the log10 of cpm.

```{r third filter method, echo=TRUE}
cpm <- cpm(mymatrix)
lcpm <- cpm(mymatrix, log=TRUE)  ### for determining the basic cutoffs
dim(lcpm)
L <- mean(mymatrix$samples$lib.size) * 1e-6
M <- median(mymatrix$samples$lib.size) * 1e-6
c(L, M)

filcpm_matrix <- subset(lcpm, (rowMeans(lcpm)>0))

dim(filcpm_matrix)
##method 2 with rowMeans
row_means <- rowMeans(lcpm)
x <- mymatrix[row_means > 0,]
dim(x)
write.csv(x$counts, "data/raw_counts.csv")


```

Both of the above methods leave 14823 genes from 28,395. I prefer the
second method, which keeps the DGEList format of the data.

```{r plots of before and after filter, echo=FALSE, warning=FALSE}
###
hist(lcpm,  main = "Histogram of total counts (unfiltered)", 
     xlab =expression("Log"[10]*" counts-per-million"), col =4 )
hist(cpm(filter_test$counts, log = TRUE), 
     main = "Histogram of filtered counts using filterByExpr method", 
     xlab =expression("Log"[10]*" counts-per-million"), 
     col = 5 )
hist(cpm(x$counts, log = TRUE), main = "Histogram of filtered counts using rowMeans > 0 method", 
     xlab =expression("Log"[10]*" counts-per-million"), col =2 ) 

```

```{r calculations, echo= FALSE, eval = FALSE}
# x <- calcNormFactors(x, method = "TMM")
#x$samples$norm.factors
#x$samples  ###checking to make sure all is well
```

```{r boxplot counts per sample, echo=FALSE, message=FALSE, warning=FALSE}


boxplot(lcpm, main = "Boxplots of log cpm per sample",xaxt = "n", xlab= "")#, xlab=mymatrix$samples, srt =85,  adj=1, cex = 0.5)
x_axis_labels(labels = samplenames, every_nth = 1, adj=1, srt =90, cex =0.4)
boxplot(filcpm_matrix, main ="boxplots of log cpm per sample filtered",xaxt = "n", xlab="")
x_axis_labels(labels = samplenames, every_nth = 1, adj=1, srt =90, cex =0.4)
```





 now I will produce the RIN x sample plots:###
```{r sample quality checks, echo=FALSE, message=FALSE, warning=FALSE}
ggplot(x$samples, aes(x = as.factor(indv), y = RIN)) +
geom_boxplot(aes(fill=as.factor(indv)))+ theme_cowplot(font_size =12, )+
labs(x= "", fill ="Individual",y ="RNA Integrity Number")+ ggtitle("Average RIN by Individual")

ggplot(x$samples, aes(x = as.factor(time), y = RIN)) +
geom_boxplot(aes(fill=as.factor(time)))+ theme_cowplot(font_size =12, )+
labs(x= "", fill ="Time in hours",y ="RNA Integrity Number")+ ggtitle("Boxplot of RIN by time")

ggplot(x$samples, aes(x = as.factor(time), y = RIN)) +
geom_boxplot(aes(fill=as.factor(time)))+ theme_cowplot(font_size =12, )+
labs(x= "", fill ="Time in hours",y ="RNA Integrity Number")+ ggtitle("Boxplot of RIN by time and drug")+
facet_wrap(~drug)

ggplot(x$samples, aes(x = as.factor(time), y = RIN)) +
geom_boxplot(aes(fill=as.factor(time)))+ theme_cowplot(font_size =12, )+
labs(x= "", fill ="Time in hours",y ="RNA Integrity Number")+ ggtitle("Boxplot of RIN by individual and time")+
facet_wrap(~indv)

ggplot(x$samples, aes(x = as.factor(drug), y = c(RIN))) +
geom_boxplot(aes(fill=as.factor(drug)))+ 
  theme_cowplot(font_size =12, )+ labs(x= "", fill ="Drug", y ="RNA Integrity Number")+ 
  theme(axis.text.x = element_blank()) + ggtitle("Boxplot of RIN by drug and individual")+ facet_wrap(~indv)
```

# Principal Component Analysis  

PCA was done using code adopted from J. Blischak.   
```{r individual PCA by drug and time steps, echo=FALSE, warning=FALSE, fig.height=8, fig.width= 6}

label <- (interaction(substring(drug, 0, 2), indv, time))
smlabel <- (interaction(substring(drug, 0, 2), time))
x <- calcNormFactors(x, method = "TMM")
normalized_lib_size <- x$samples$lib.size * x$samples$norm.factors
dat_cpm <- cpm(x$counts, lib.size = normalized_lib_size)
colnames(dat_cpm) <- label
anno$time <- ordered(time, levels =c("3h", "24h"))
anno$group <- group
cpm_per_sample <- cbind(anno, t(dat_cpm)) 
treattype <-c("Daunorubicin","Doxorubicin","Epirubicin","Mitoxantrone","Trastuzumab","Vehicle")

for (b in treattype) { 
  dat <- dat_cpm[,anno$drug %in% c("none", b)] 
  cat("\n\n### ", b, "\n\n") 
  # PCA 
  pca <- calc_pca(t(dat))$x
  pca <- data.frame(anno[anno$drug %in% c("none", b), c("drug", "time","indv")], pca)
  pca <- droplevels(pca) 
  for (cat_var in c("drug", "time", "indv")) {
    assign(paste0(cat_var, "_pca"), arrangeGrob(pca_plot(pca, cat_var)))
  } 
  drug_time_pca <- pca_plot(pca, "drug", "time")
 
  grid.arrange(drug_time_pca, indv_pca, drug_pca, nrow =3, top =paste(b)) 
  }
write.csv(x$counts, "data/norm_counts.csv")

```

# PCA of all 72 samples  


```{r all-samples, echo=FALSE}
pca_all <- calc_pca(t(dat_cpm))
pca_all_anno <- data.frame(anno, pca_all$x)
print(pca_var_plot(pca_all))
```
## Variance contribution from treatment, extraction time, or individual on PC1 and PC2

```{r variance by treatment time and individual, echo=FALSE, message=FALSE, warning=FALSE}
facs <- c("indv", "drug", "time")
names(facs) <- c("Individual", "Treatment", "Time")
drug1 <- unique(drug)##for changing shapes and colors
for (f in facs) {
  # PC1 v PC2
  pca_plot(pca_all_anno, col_var = f, shape_var = time,
           title = names(facs)[which(facs == f)])
  print(last_plot())
  
  # Plot f versus PC1 and PC2
  f_v_pc1 <- arrangeGrob(plot_versus_pc(pca_all_anno, 1, f))
  f_v_pc2 <- arrangeGrob(plot_versus_pc(pca_all_anno, 2, f))
  grid.arrange(f_v_pc1, f_v_pc2, ncol = 2, top = names(facs)[which(facs == f)])
  #summary(plot_versus_pc(pca_all_anno, 1, f))
  #summary(plot_versus_pc(pca_all_anno, 2, f))
}

```

# Differential Expression

```{r Differential expression, echo=FALSE}
#mm <- model.matrix(~0 +drug +time)
#colnames(mm) <- gsub("group", "", colnames(mm))
group1 <- interaction(drug,time)
mm2 <- model.matrix(~0 + group1)
colnames(mm2) <- c("A3", "X3", "E3","M3","T3", "V3","A24", "X24", "E24","M24","T24", "V24")

y2 <- voom(x, mm2,plot =TRUE)

```



```{r  Differential expression2, eval=TRUE, include=FALSE}
###this is enabled when first running program, otherwise processing slows down


corfit2 <- duplicateCorrelation(y2, mm2, block = indv)

v2 <- voom(x, mm2, block = indv, correlation = corfit2$consensus)

fit2 <- lmFit(v2, mm2, block = indv, correlation = corfit2$consensus)

```

```{r Differential expression3 contrasts, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE}
# Make contrasts
#cm <- makeContrasts(
 # VvDA = drugVehicle -drugDaunorubicin,
 # VvDX = drugVehicle-drugDoxorubicin,
 # VvEP = drugVehicle-drugEpirubicin,
 # VvMT = drugVehicle-drugMitoxantrone,
 # VvTR = drugVehicle-drugTrastuzumab,
 #  levels=mm)
cm2 <- makeContrasts(
  V.DA = V3 - A3,
  V.DX = V3 - X3,
  V.EP = V3 - E3,
  V.MT = V3 - M3,
  V.TR = V3 - T3,
  V.DA24 = V24-A24,
  V.DX24= V24-X24,
  V.EP24= V24-E24,
  V.MT24= V24-M24,
  V.TR24= V24-T24,
  levels = mm2)

vfit2 <- lmFit(y2, mm2)

vfit2<- contrasts.fit(vfit2, contrasts=cm2)

efit2 <- eBayes(vfit2)

V.DA.top= topTable(efit2, coef=1, adjust="BH", number=Inf, sort.by="p")
V.DX.top= topTable(efit2, coef=2, adjust="BH", number=Inf, sort.by="p")
V.EP.top= topTable(efit2, coef=3, adjust="BH", number=Inf, sort.by="p")
V.MT.top= topTable(efit2, coef=4, adjust="BH", number=Inf, sort.by="p")
V.TR.top= topTable(efit2, coef=5, adjust="BH", number=Inf, sort.by="p")
V.DA24.top= topTable(efit2, coef=6, adjust="BH", number=Inf, sort.by="p")
V.DX24.top= topTable(efit2, coef=7, adjust="BH", number=Inf, sort.by="p")
V.EP24.top= topTable(efit2, coef=8, adjust="BH", number=Inf, sort.by="p")
V.MT24.top= topTable(efit2, coef=9, adjust="BH", number=Inf, sort.by="p")
V.TR24.top= topTable(efit2, coef=10, adjust="BH", number=Inf, sort.by="p")


sigVDA3 = V.DA.top[V.DA.top$adj.P.Val < .1 , ]
sigVDX3 = V.DX.top[V.DX.top$adj.P.Val < .1 , ]
sigVEP3 = V.EP.top[V.EP.top$adj.P.Val < .1 , ]
sigVMT3 = V.MT.top[V.MT.top$adj.P.Val < .1 , ]
sigVTR3 = V.TR.top[V.TR.top$adj.P.Val < .1 , ]
sigVDA24 = V.DA24.top[V.DA24.top$adj.P.Val < .1 , ]
sigVDX24 = V.DX24.top[V.DX24.top$adj.P.Val < .1 , ]
sigVEP24 = V.EP24.top[V.EP24.top$adj.P.Val < .1 , ]
sigVMT24 = V.MT24.top[V.MT24.top$adj.P.Val < .1 , ]
sigVTR24 = V.TR24.top[V.TR24.top$adj.P.Val < .1 , ]

siglist <- list(sigVDA3,
             sigVDX3,
             sigVEP3,
             sigVMT3,
             sigVTR3,
             sigVDA24,
             sigVDX24,
             sigVEP24,
             sigVMT24,
             sigVTR24)
names(siglist) <- c("sigVDA3",
             "sigVDX3",
             "sigVEP3",
             "sigVMT3",
             "sigVTR3",
             "sigVDA24",
             "sigVDX24",
             "sigVEP24",
             "sigVMT24",
             "sigVTR24")
siglist <- lapply(siglist, function(x) {colnames(x) <- gsub(pattern = "^siglist\\.", replacement = "", x = colnames(x))
return(x)})


#make the list, and rename the items on the list
 ##print count files to data directory (write.table needs the full file name) 
# for(i in seq_along(siglist)) {
#   fnl <- paste0("data/",names(siglist)[i],".txt")
#     write.csv(siglist[i], fnl, col.names = TRUE, row.names =FALSE, sep="") 
#              
# }

#write.csv(names(siglist),data/siglist.txt", col.names = FALSE,row.names = FALSE)
```

```{r sample code of DEG setup, eval=FALSE, message=FALSE, warning=FALSE, echo=TRUE}
mm2 <- model.matrix(~0 + group1)
##made the matrix model using the interaction between Treatment and Time
colnames(mm2) <- c("A3", "X3", "E3","M3","T3", "V3","A24", "X24", "E24","M24","T24", "V24")
y2 <- voom(x, mm2)

corfit2 <- duplicateCorrelation(y2, mm2, block = indv)

v2 <- voom(x, mm2, block = indv, correlation = corfit2$consensus)

fit2 <- lmFit(v2, mm2, block = indv, correlation = corfit2$consensus)

vfit2 <- lmFit(v2, mm2)

vfit2<- contrasts.fit(vfit2, contrasts=cm2)

efit2 <- eBayes(vfit2)

V.DA.top= topTable(efit2, coef=1, adjust="BH", number=Inf, sort.by="p")
### sorting all top expressed genes for the Vehicle and Daunorubicin 3 hour treatments
sigVDA3 = V.DA.top[V.DA.top$adj.P.Val < .1 , ]
### this helped pull only those files that were at an adjusted p value of less than 0.1
### This p-value was used as the beginning examination of the data, considering I will run multiple runs of this RNA seq library.
```

This is the example code I used to process my data.  I used two model matrix initially, one set up was /~0 +drug +time and the second was /~0+group1, then blocking by individual.
That is why you see the number 2 in the code above.  


#DEG summary  

```{r prcompish stuff, echo=FALSE, warning=FALSE, paged.print=TRUE}
prdat_cpm <- calc_pca(t(dat_cpm)) 
summary(prdat_cpm)
 sum <- summary(decideTests(efit2))

sum
 
# biplot(prdat_cpm, ind.names = FALSE)
#  biplot(prdat_cpm, y=choices = 1:2, scale = 1, pc.biplot = FALSE)

```




## Graphing DE genes  

I then created a counts table for each set of genes.  Luckily, the counts are stored in the y2 object, which is an EList class object. I can 'simplify' this process because I kept the DEGList format initially.
I first made an object called 'countstotal' from the EList y2. For ggploting later, I subsetted 'countstotal' by treatments.

```{r counts files, message=FALSE, warning=FALSE}
countstotal <- y2$E
colnames(countstotal) <- smlabel
boxplot(countstotal, xaxt = "n", xlab="")
Da3counts <- as.data.frame(as.table(countstotal[,c(1,6,13,18,25,30,37,42,49,54,61,66)]))
x_axis_labels(labels = label, every_nth = 1, adj=1, srt =90, cex =0.4)
```

```{r processing files, message=FALSE, warning=FALSE, include=FALSE}
Da3counts <- as.data.frame(as.table(countstotal[,c(1,6,13,18,25,30,37,42,49,54,61,66)]))
Dx3counts <- as.data.frame(as.table(countstotal[,c(2,6,14,18,26,30,38,42,50,54,62,66)]))
Ep3counts <- as.data.frame(as.table(countstotal[,c(3,6,15,18,27,30,39,42,51,54,63,66)]))
Mt3counts <- as.data.frame(as.table(countstotal[,c(4,6,16,18,28,30,40,42,52,54,64,66)]))
TR3counts <- as.data.frame(as.table(countstotal[,c(5,6,17,18,29,30,41,42,53,54,65,66)]))

mDa3counts <- countstotal[,c(1,6,13,18,25,30,37,42,49,54,61,66)]
mDx3counts <- countstotal[,c(2,6,14,18,26,30,38,42,50,54,62,66)]
mEp3counts <- countstotal[,c(3,6,15,18,27,30,39,42,51,54,63,66)]
mMt3counts <- countstotal[,c(4,6,16,18,28,30,40,42,52,54,64,66)]
mTR3counts <- countstotal[,c(5,6,17,18,29,30,41,42,53,54,65,66)]

Da24counts <- as.data.frame(as.table(countstotal[,c(7,12,19,24,31,36,43,48,55,60,67,72)]))
Dx24counts <- as.data.frame(as.table(countstotal[,c(8,12,20,24,32,36,44,48,56,60,68,72)]))
Ep24counts <- as.data.frame(as.table(countstotal[,c(9,12,21,24,33,36,45,48,57,60,69,72)]))
Mt24counts <- as.data.frame(as.table(countstotal[,c(10,12,22,24,34,36,46,48,58,60,70,72)]))
TR24counts <- as.data.frame(as.table(countstotal[,c(11,12,23,24,35,36,47,48,59,60,71,72)]))

mDa24counts <- countstotal[,c(7,12,19,24,31,36,43,48,55,60,67,72)]
mDx24counts <- countstotal[,c(8,12,20,24,32,36,44,48,56,60,68,72)]
mEp24counts <- countstotal[,c(9,12,21,24,33,36,45,48,57,60,69,72)]
mMt24counts <- countstotal[,c(10,12,22,24,34,36,46,48,58,60,70,72)]
mTR24counts <- countstotal[,c(11,12,23,24,35,36,47,48,59,60,71,72)]

colnames(Da24counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Dx24counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Ep24counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Mt24counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(TR24counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Da3counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Dx3counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Ep3counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(Mt3counts) <- c("ENTREZID", "Samples",  "Counts")
colnames(TR3counts) <- c("ENTREZID", "Samples",  "Counts")

write_csv(Da24counts,"data/Da24counts.txt")
write_csv(Dx24counts,"data/Dx24counts.txt")
write_csv(Ep24counts,"data/Ep24counts.txt")
write_csv(Mt24counts,"data/Mt24counts.txt")
write_csv(TR24counts,"data/TR24counts.txt")


```
library(cowplot)

# Volcano Plots  

```{r Volcanoplots, echo=FALSE, message=FALSE, warning=FALSE}
theme_set(theme_cowplot())
plot_filenames <- c("V.DA.top","V.DX.top","V.EP.top","V.MT.top","V.TR.top","V.DA24.top",
"V.DX24.top","V.EP24.top","V.MT24.top","V.TR24.top")

volcanosig <- function(df, psig.lvl,topg) {
    df <- df %>% 
    mutate(threshold = ifelse(adj.P.Val > psig.lvl, "A", ifelse(adj.P.Val <= psig.lvl & logFC<=0,"B","C")))
      # ifelse(adj.P.Val <= psig.lvl & logFC >= 0,"B", "C")))
    
    df <- df %>% mutate(genelabels = "")
    df$genelabels[1:topg] <- df$SYMBOL[1:topg]
    
  ggplot(df, aes(x=logFC, y=-log10(P.Value))) + 
    geom_point(aes(color=threshold))+
    geom_text_repel(aes(label = genelabels), segment.curvature = -1e-20,force = 1,size=2.5,
    arrow = arrow(length = unit(0.015, "npc")), max.overlaps = Inf) +
    #geom_hline(yintercept = -log10(psig.lvl))+
    xlab(expression("Log"[2]*" FC"))+
    ylab(expression("-log"[10]*"P Value"))+
    scale_color_manual(values = c("black", "red","blue"))+
    theme_cowplot()+
    theme(legend.position = "none",
              plot.title = element_text(size = rel(1.5), hjust = 0.5),
              axis.title = element_text(size = rel(0.8))) 
}
#v1<- volcanosig(V.DA24.top, 0.01,0)
v1 <- volcanosig(V.DA.top, 0.01,10)+ ggtitle("Daunorubicin 3 hour")
v2 <- volcanosig(V.DA24.top, 0.01,10)+ ggtitle("Daunorubicin 24 hour")+ylab("")
v3 <- volcanosig(V.DX.top, 0.01,10)+ ggtitle("Doxorubicin 3 hour")
v4 <- volcanosig(V.DX24.top, 0.01,10)+ ggtitle("Doxorubicin 24 hour")+ylab("")
v5 <- volcanosig(V.EP.top, 0.01,10)+ ggtitle("Epirubicin 3 hour")
v6 <- volcanosig(V.EP24.top, 0.01,10)+ ggtitle("Epirubicin 24 hour")+ylab("")
v7 <- volcanosig(V.MT.top, 0.01,10)+ ggtitle("Mitoxatrone 3 hour")
v8 <- volcanosig(V.MT24.top, 0.01,10)+ ggtitle("Mitoxatrone 24 hour")+ylab("")
v9 <- volcanosig(V.TR.top, 0.01,10)+ ggtitle("Trastuzumab 3 hour")
v10 <- volcanosig(V.TR24.top, 0.01,10)+ ggtitle("Trastuzumab 24 hour")+ylab("")
# volcanoplot(efit2,coef = 10,style = "p-value",
#                           highlight = 8,
#                           names = efit2$genes$SYMBOL,
#                           hl.col = "red",xlab = "Log2 Fold Change",
#                           ylab = NULL,pch = 16,cex = 0.35,
#                           main = "Using Trastuzuman 24 hour data and volcanoplot function")
plot_grid(v1,v2,  rel_widths =c(.8,1))
plot_grid(v3,v4,  rel_widths =c(.8,1))
plot_grid(v5,v6,  rel_widths =c(.8,1))
plot_grid(v7,v8,  rel_widths =c(.8,1))
plot_grid(v9,v10,  rel_widths =c(.8,1))
```

```{r cormotif things, echo=FALSE, message=FALSE, warning=FALSE, paged.print=TRUE}
library(Cormotif)
library(pheatmap)


colnames(countstotal) <- label





mcor <- cor(countstotal)
#mcor
#pheatmap::pheatmap(mm2, cluster_rows = FALSE, cluster_cols = FALSE)
corrcounts <- pheatmap::pheatmap(countstotal, cluster_rows =TRUE, cluster_cols = TRUE)
#symnum(mcor)
#range(mcor)
corrplot(cor(countstotal), method = 'color',
         type = 'lower', 
         order = 'hclust',
         hclust.method = "complete",
         is.corr = TRUE,
         tl.col = 'black',
         cl.ratio = 0.2, 
         tl.srt = 45,
         tl.cex =.5, 
         col = COL2('RdYlBu', 35), 
         cl.pos = 'b',
         col.lim = c(.65,1))


pheatmap::pheatmap(mcor, cluster_rows = TRUE, cluster_cols = TRUE)



```



# [GO analysis](https://reneeisnowhere.github.io/Cardiotoxicity/GO_analysis.html)