---
title: "Gene expression, lab4"
author: "Jeff Kinne"
date: "2024-06-26"
output: 
  html_document:
    code_folding: hide
    code_download: true
---

```{r setup, include=FALSE}
#install.packages(c("rmarkdown", "knitr"))
```

# R Markdown

This is an R Markdown file. It allows us to have code, analysis, and documentation all in one file. See https://cs.indstate.edu/wiki/index.php/Bioinformatics#R_Programming for some links to look at about R markdown.

# This lab - compare results

In this lab we compare results between the two datasets GSE85331 and GSE244362. Some parts are completed for you, and you will complete the rest. 

Go back to lab2 and lab3 and save the dataframes all_info from each as these filenames: gse85331_all_info.csv and gse244362_all_info.csv. You should then load these files here as gse85331 and gse244362. Here I have given you a start by showing how I would like this done for gse85331.
```{r load data}
# change this line to point to the directory where you keep your data files
setwd("/Users/jkinne/OneDrive - Indiana State University/Courses/cs618")

# list to keep information about the gse85331 dataset
gse85331 <- list()
gse244362 <- list()

# read in the results from lab2 and save into our list
gse85331$data_org <- read.csv("gse85331_all_info.csv")

# update the column name for the gene symbol to be "gene_symbol"
colnames(gse85331$data_org)[1] <- "gene_symbol"

# let's see a summary of what we have loaded
summary(gse85331$data_org)

# lab4 - do the same for gse244361

```

Next we put a few different plots from the datasets. Let's look at a histogram of the p-values.

```{r plots}
par(mfrow=c(1,2)) # make it so plots will be side by side
hist(gse85331$data_org$p_adj, main="gse85331 adjusted p values")
# lab4 - histogram for gse244362 adjusted p values
par(mfrow=c(1,1))
```
Now let's look at genes that are in common between the two datasets.
```{r comparing results}
# a function that takes two vectors and prints out information about 
#  how many things are in common with each.
compare_results <- function(data1, text1, data2, text2) {
  print(paste(text1, "and", text2))
  print(paste(text1, "number of rows:", length(data1), "(", sum(data1 %in% data2), "also in", text2,")"))
  print(paste(text2, "number of rows:", length(data2), "(", sum(data2 %in% data1), "also in", text1, ")"))
}

# lab4 - uncomment these lines once you have gse244362 loaded properly
##compare_results(gse85331$data_org$gene_symbol, "gse85331", gse244362$data_org$gene_symbol, "gse244362")

# And now we're going to just look at the rows they have in common.
gse85331$data <- gse85331$data_org # lab4 - remove once you uncomment the next one
##gse85331$data <- gse85331$data_org[gse85331$data_org$gene_symbol %in% gse244362$data_org$gene_symbol,]
##gse244362$data <- gse244362$data_org[gse244362$data_org$gene_symbol %in% gse85331$data_org$gene_symbol,]

# Let's put them in the same order too so we can easily compare rows.
# First, gse244362 has a few genes that showed up more than once, so going to limit those to just 1.
##gse244362$data <- gse244362$data[! duplicated(gse244362$data$gene_symbol),]
##gse244362$data <- gse244362$data[order(gse244362$data$gene_symbol),]
##gse85331$data <- gse85331$data[order(gse85331$data$gene_symbol),]
##sum(duplicated(gse244362$data$gene_symbol)) # should be 0
##sum(duplicated(gse85331$data$gene_symbol))  # should be 0
##sum(gse85331$data$gene_symbol == gse244362$data$gene_symbol) # should be same as number of rows
##dim(gse85331$data)

# the columns that are the original data items
gse85331$data_cols <- (ncol(gse85331$data)-31):ncol(gse85331$data)
# lab4 - same for gse244362
```
Now we have two datasets with the same genes and in the same order. Let's look at them. We can look at the difference in the p values for each gene (which we expect to generally be close to 0 - the two different datasets are from doing the same procedure so should agree most of the time on whether a gene is signifcant or not). And we can look at boxplots of the samples, to see if the range of values is similar (if it is not, then we definitely would want to rescale the data if we were going to compare points the expression values on a gene between the two datasets).
```{r}
# lab4 - uncomment
##hist(gse85331$data$p_adj-gse244362$data$p_adj, breaks=100, main="Difference in adjusted p values")

par(mfrow=c(1,2))
boxplot(gse85331$data[,gse85331$data_cols], main = "gse85331 expression values")
# lab4 - same for gse244362
par(mfrow=c(1,1))

par(mfrow=c(1,2))
boxplot(colSums(gse85331$data_org[,gse85331$data_cols]), main="gse85331 colSums")
# lab4 - same for gse244362
par(mfrow=c(1,1))

```
We can also look at the correlation between the median values of the time points in the two different datasets. We would expect a high correlation (close to 1) for most of the genes. 
```{r}
# the columns that are the medians of the 4 time points, in the right order
gse85331$group_columns <- c(2,3,4,5)
# lab4 - same for gse244362

# lab4 - uncomment 

# compute the correlation values - which we do separately for each gene.
# note that cor compares by column, so we use t to transpose a row to be a column.
# suppressWarnings so we don't see the warning messages about 0 variance.
##cor_values <- suppressWarnings(sapply(1:nrow(gse85331$data), function(x) {
##  cor(t(gse85331$data[x,gse85331$group_columns,]),
##      t(gse244362$data[x,gse244362$group_columns]))
##  }))
##hist(cor_values, main="Correlations between the datasets")

# compare the previous histogram with what we get if we just used completely random data
cor_random_values <- suppressWarnings(sapply(1:nrow(gse85331$data), function(x) {
  # generate two vectors that are each 4 points, with a mean that is a random 
  #  value between 0 and 10 and has a variance that is a random value between 0.1 and 2
  cor(rnorm(4,runif(1,0,10), runif(1,0.1,2)), rnorm(4,runif(1,0,10), runif(1,0.1,2)))
}))
hist(cor_random_values, main="Correlations of randomly generated values")
```
We conclude from all of that, that these two datasets look pretty similar. Let's see how similar our results are from doing the very simple checks we did in labs 2 and 3. We will load the results for each of the different comparisons we did and see how many of the results are the same. I'll leave one example of how I did this. 
```{r}
gse85331$day2_up <- read.csv("day2_up.csv")
colnames(gse85331$day2_up)[1] <- "gene_symbol"

gse244362$day2_up <- read.csv("D1_up.csv")
colnames(gse244362$day2_up)[2] <- "gene_symbol"

compare_results(gse85331$day2_up$gene_symbol, "gse85331 day2_up genes", 
                gse244362$day2_up$gene_symbol, "gse244362 day2_up genes")

take_intersect <- function(name, d1, d2, name2) {
  print(paste0(name, ": ", nrow(d1)))
  d1 <- d1[d1[,"gene_symbol"] %in% d2[,"gene_symbol"],]
  print(paste0("      ", nrow(d1), " are also in ", name2))
  d1
}

# lab4 - uncomment the following after you have gse244362 loaded, and do similar
#  lines for the other results groups.

##print("gse85331...")
##gse85331$day2_up <- take_intersect("day2_up", gse85331$day2_up, gse244362$data_org, "gse244362")

##print("gse244362...")
##gse244362$day2_up <- take_intersect("day2_up", gse244362$day2_up, gse85331$data_org, "gse85331")

```