Unsupervised MDS analysis based on combinatorial chromatin states from whole crypts in fed vs fasted mice

Christopher Terranova May 2019 based on code from Ayush Raman

Load necessary programs

library(circlize)
library(ComplexHeatmap)
library(pheatmap)
library(dplyr)
library(gplots)
library(ggplot2)
library(genefilter)
library(RColorBrewer)
library(reshape2)
library(vegan)

Identification of most variable chromatin state regions

# set working directory to combined matrix state file from chromStatesMatrixforClustering_particularStates.pl outputsummary
  setwd("~/Desktop/KS_SMINT_ChromXR/")

# set chromatin state for analysis -- can be changed for each chromatin state
  i = 3

# read in combined matrix state file from chromStatesMatrixforClustering_particularStates.pl output
dat <- read.table(paste("CombinedMatrix-",i,"-10000bps.txt",sep =""),header = FALSE,sep = "\t",quote = "",
                    row.names = 1, na.strings = FALSE, stringsAsFactors = FALSE)
  colnames(dat) <- c("NMRL-10F1-1","NMRL-11F4-1","NMRL-10S2-1","NMRL-11S5-1")
  head(dat)

##                          NMRL-10F1-1 NMRL-11F4-1 NMRL-10S2-1 NMRL-11S5-1
## chr1:10230000-10239999             1           0           0           1
## chr1:10240000-10249999             0           0           1           1
## chr1:10250000-10259999             0           0           1           1
## chr1:106930000-106939999           1           0           0           0
## chr1:106980000-106989999           0           1           1           1
## chr1:106990000-106999999           0           1           2           2

# sample type labels for plot
  condition.type <- factor(c(rep("FED",2), rep("FASTED",2)))
  
# filtering low variable regions
  dat.state1 <- varFilter(as.matrix(dat), var.cutoff = 0.9)
  print(dim(dat.state1))

## [1] 700   4

MDS plot based on active chromatin state in fed vs fasted mice

d <- dist(t(dat.state1), method = "binary")
  mds.plot <- cmdscale(d = d, eig = T, k = 2)
  mdsDist <- data.frame(genotypes = condition.type, x = mds.plot$points[,1], y = mds.plot$points[,2])
  print(ggplot(mdsDist, aes(x = x, y = y, color = genotypes)) + geom_point(size = 8) +
          ylab("MDS Coordinate 2") + xlab("MDS Coordinate 1") + theme_grey() +
          theme(legend.text = element_text(size = 18, face = "bold"),
                legend.title = element_text(size = 18, colour = "black", face = "bold"),
                axis.title = element_text(size = 18, face = "bold"),
                axis.text.x = element_text(size = 18, face = "bold", color = "black"),
                axis.text.y = element_text(size = 18, face = "bold", color = "black"),
                plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm")))

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.