1.2_timepointSpecificGeneValidation.Rmd

---
title: "RNA-seq Mouse Placenta"
author: "Ha T. H. Vu"
output: html_document
---
  
```{r setup, include=FALSE}
options(max.print = "75")
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>",
  fig.path = "Files/",
  fig.width = 8,
  prompt = FALSE,
  tidy = FALSE,
  message = FALSE,
  warning = TRUE
)
knitr::opts_knit$set(width = 75)
```

Documentation for the analysis of mouse placental RNA-seq data at e7.5, e8.5 and e9.5.

## 1. Clustering analysis to confirm hierarchical clustering trends
This analysis is carried out in two steps. First, run the clustering algorithms with different number of groups (k = 3, 4 and 5). Second, get a summary table with expression levels of transcripts to validate the trends of the groups with visualization. For the second step, we can set up the following functions:
```{r}
summ <- function(trans_cluster, tpm2=tpm2) {
  tpm3 <- tpm2
  tpm3 <- cbind(rownames(tpm2), tpm3)
  rownames(tpm3) <- 1:nrow(tpm3)
  colnames(tpm3)[1] <- "transcripts"

  #e7.5
  e7.5_table <- data.frame(tpm3$transcripts, rowMeans(tpm3[,c("E7.5_2", "E7.5_3", "E7.5_4", "E7.5_5", "E7.5_6")]))
  e7.5_table$time <- c("e7.5")
  rownames(e7.5_table) <- 1:nrow(e7.5_table)
  colnames(e7.5_table) <- c("transcripts", "mean_cts_scaled", "time")

  #e8.5
  e8.5_table <- data.frame(tpm3$transcripts, rowMeans(tpm3[,c("E8.5_1", "E8.5_2", "E8.5_3", "E8.5_4", "E8.5_5", "E8.5_6")]))
  e8.5_table$time <- c("e8.5")
  rownames(e8.5_table) <- 1:nrow(e8.5_table)
  colnames(e8.5_table) <- c("transcripts", "mean_cts_scaled", "time")

  #e9.5
  e9.5_table <- data.frame(tpm3$transcripts, rowMeans(tpm3[,c("E9.5_1", "E9.5_2", "E9.5_3", "E9.5_4", "E9.5_5")]))
  e9.5_table$time <- c("e9.5")
  rownames(e9.5_table) <- 1:nrow(e9.5_table)
  colnames(e9.5_table) <- c("transcripts", "mean_cts_scaled", "time")

  summary <- rbind(rbind(e7.5_table, e8.5_table), e9.5_table)
  summary <- summary[order(summary$transcripts),]
  summary <- inner_join(summary, trans_cluster, by = c("transcripts" = "name"))
  summary <- inner_join(summary, t2g, by = c("transcripts" = "target_id"))
  
  return(summary)
}
```
```{r}
plotClus <- function(summary, title){
  ascl2 <- subset(summary, summary$ext_gene == "Ascl2") #e7.5
  gjb5 <- subset(summary, summary$ext_gene == "Gjb5") #e7.5
  dnmt1 <- subset(summary, summary$ext_gene == "Dnmt1") #e8.5
  itga4 <- subset(summary, summary$ext_gene == "Itga4") #e8.5
  gjb2 <- subset(summary, summary$ext_gene == "Gjb2") #e9.5
  igf2 <- subset(summary, summary$ext_gene == "Igf2") #e9.5
  p <- ggplot(aes(time, mean_cts_scaled), data = summary) +
    geom_line(aes(group = transcripts), alpha = 0.5, colour = "grey77") +
    geom_line(stat = "summary", fun = "median", size = 2,
              aes(group = 1, color = "Group median")) +
    labs(title = title,
         x = "Time point",
         y = "Scaled mean transcript counts", color = "Legend", linetype = "Legend") +
    theme(plot.title = element_text(size = 15, face = "bold"), legend.text=element_text(size=20)) +
    geom_line(data = ascl2, size = 2.5, aes(group = transcripts, color = "Ascl2", linetype = "Ascl2"), alpha = 1) + #e7.5
    geom_line(data = gjb5, size = 2, aes(group = transcripts, color = "Gjb5", linetype = "Gjb5" ), alpha = 1) + #e7.5
    geom_line(data = dnmt1, size = 2, aes(group = transcripts, color = "Dnmt1", linetype = "Dnmt1"), alpha = 1) + #e8.5
    geom_line(data = itga4, size = 2, aes(group = transcripts, color = "Itga4", linetype = "Itga4"), alpha = 1) + #e8.5
    geom_line(data = gjb2, size = 2, aes(group = transcripts, color = "Gjb2", linetype = "Gjb2"), alpha = 1) + #e9.5
    geom_line(data = igf2, size = 2, aes(group = transcripts, color = "Igf2", linetype = "Igf2"), alpha = 1) + #e9.5
    scale_color_manual(name = "Legend", values = c("Ascl2" = "darkolivegreen4", "Gjb5" = "yellow3",
                                                   "Dnmt1" = "dodgerblue4", "Itga4" = "deepskyblue4",
                                                   "Gjb2" = "saddlebrown", "Igf2" = "salmon3",
                                                   "Group median" = "grey22")) + 
    scale_linetype_manual(name = "Legend", values = c("Ascl2" = "solid", "Gjb5" = "solid",
                                                      "Dnmt1" = "twodash", "Itga4" = "solid",
                                                      "Gjb2" = "solid", "Igf2" = "solid",
                                                      "Group median" = "solid")) + 
      guides(linetype=F,
             colour=guide_legend(keywidth = 3, keyheight = 1)) +
      theme(text = element_text(size=15),
            legend.text=element_text(size=15),
            axis.title.x = element_text(size = 15),
            axis.title.y = element_text(size = 15),
            axis.text.x = element_text(angle=0, hjust=0.5, size = 25),
            axis.text.y = element_text(size = 15)) +
      facet_grid(cols = vars(value))
  return(p)
}
```

Next, load some necessary files:
```{r}
library("ggplot2", suppressMessages())
library("dplyr", suppressMessages())
library("tidyverse", suppressMessages())

t2g <- read.table("Files/t2g.txt", header = T, sep = "\t")
t2g <- t2g[order(t2g$target_id),]
coding <- read.table("Files/Mus_musculus_grcm38_coding_transcripts.txt", header = F)
tpm2 <- read.table("/work/LAS/geetu-lab/hhvu/project1_2/rna-seq/PlacentaRNA-seq/Files/tpmForClustering.txt", header = T)
```

### k = 3:
Hierarchical clustering, k = 3:
```{r}
hc <- hclust(dist(tpm2, method = "euclidean"), "complete")
trans_cluster <- cutree(hc, k = 3) %>% enframe()
table(trans_cluster$value)
summK3 <- summ(trans_cluster, tpm2)
summK3$value <- ifelse(summK3$value == "1", "e8.5",
                        ifelse(summK3$value == "2", "e9.5",
                               ifelse(summK3$value == "3", "e7.5", summK3$value)))
summK3$value <- ifelse(summK3$value == "e8.5", "2",
                        ifelse(summK3$value == "e9.5", "3",
                               ifelse(summK3$value == "e7.5", "1", summK3$value)))
plotClus(summK3, "Hierarchical Clustering of Transcripts, k = 3")

summK32 <- dplyr::distinct(summK3[,c("value", "transcripts")])
summK32 <- summK32[order(summK32$transcripts),]
```

For K-means clustering: <br>
```{r}
set.seed(123)
km <- kmeans(tpm2, centers = 3)
trans_cluster_kmeans <- km$cluster %>% enframe()
summKmeans <- summ(trans_cluster_kmeans, tpm2)
summKmeans$value <- ifelse(summKmeans$value == "1", "e9.5",
                        ifelse(summKmeans$value == "2", "e7.5",
                               ifelse(summKmeans$value == "3", "e8.5", summKmeans$value)))
summKmeans$value <- ifelse(summKmeans$value == "e8.5", "2",
                        ifelse(summKmeans$value == "e9.5", "3",
                               ifelse(summKmeans$value == "e7.5", "1", summKmeans$value)))
plotClus(summKmeans, "K-means Clustering of Transcripts, k = 3")

summKmeans2 <- dplyr::distinct(summKmeans[,c("value", "transcripts")])
summKmeans2 <- summKmeans2[order(summKmeans2$transcripts),]

t <- summKmeans2
for (i in 1:3) {
  for (j in 1:3) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent agreement between hierCluster", i, "and Kmeans cluster", j, "is:", q/length(unique(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"]))*100))
  }
}
```

For self-organizing maps: <br>
```{r, eval=F}
library("kohonen")
set.seed(123)
som <- som(as.matrix(tpm2), grid = somgrid(3, 1, "rectangular"))
trans_cluster_som <- som$unit.classif %>% enframe()
trans_cluster_som$name <- row.names(tpm2)
```
Since this run could take time, I saved the results to reproduce the plot here. <br>
```{r}
trans_cluster_som <- read.table("Files/trans_cluster_som.txt", header = T)
summSOM <- summ(trans_cluster_som, tpm2)
summSOM$value <- ifelse(summSOM$value == "1", "e8.5",
                        ifelse(summSOM$value == "2", "e9.5",
                               ifelse(summSOM$value == "3", "e7.5", summSOM$value)))
summSOM$value <- ifelse(summSOM$value == "e8.5", "2",
                        ifelse(summSOM$value == "e9.5", "3",
                               ifelse(summSOM$value == "e7.5", "1", summSOM$value)))
plotClus(summSOM, "Self-organizing Maps of Transcripts, k = 3")

summSOM2 <- dplyr::distinct(summSOM[,c("value", "transcripts")])
summSOM2 <- summSOM2[order(summSOM2$transcripts),]

t <- summSOM2
for (i in 1:3) {
  for (j in 1:3) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in SOM cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

For spectral clustering: <br>
This clustering method required high memory and computing time, so we should run this analysis on a HPC cluster. If you are reproducing this analysis from an Iowa State cluster such as pronto, the R library `Rclustool` would require to run on an *interactive* node with X11 forwarding enabled (see https://researchit.las.iastate.edu/x-forwarding-mac-and-windows). Additionally, request high memory (for example, 128G) to run the analysis. It can take a few hours to finish running. <br>
```{r, eval=F}
library("RclusTool")
set.seed(123)
sim <- computeGaussianSimilarity(tpm2, 1)
res <- spectralClustering(sim, K=3)
trans_cluster_sc <- data_frame(row.names(tpm2), res$label)
colnames(trans_cluster_sc) <- c("name", "value")
```
Plot of transcripts: <br>
```{r}
trans_cluster_sc <- read.table("Files/trans_cluster_sc.txt", header = T)
summSC <- summ(trans_cluster_sc, tpm2)
summSC$value <- ifelse(summSC$value == "1", "e7.5",
                        ifelse(summSC$value == "2", "e8.5",
                               ifelse(summSC$value == "3", "e9.5", summSC$value)))
summSC$value <- ifelse(summSC$value == "e8.5", "2",
                        ifelse(summSC$value == "e9.5", "3",
                               ifelse(summSC$value == "e7.5", "1", summSC$value)))
plotClus(summSC, "Spectral Clustering of Transcripts")

summSC2 <- dplyr::distinct(summSC[,c("value", "transcripts")])
summSC2 <- summSC2[order(summSC2$transcripts),]

t <- summSC2
for (i in 1:3) {
  for (j in 1:3) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent agreement between hierCluster", i, "and SC cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```
We can see the cluster trends agree with these of hierarchical clustering. <br>

### k = 4:
Hierarchical clustering, k = 4:
```{r}
trans_cluster <- cutree(hc, k = 4) %>% enframe()
write.table(trans_cluster, "Files/trans_cluster_hc_k4.txt", quote = F, row.names = F, sep = "\t")
table(trans_cluster$value)
summK3 <- summ(trans_cluster, tpm2)
summK3$value <- ifelse(summK3$value == "1", "e8.5",
                        ifelse(summK3$value == "2", "e9.5",
                               ifelse(summK3$value == "3", "e7.5", summK3$value)))
summK3$value <- ifelse(summK3$value == "e8.5", "2",
                        ifelse(summK3$value == "e9.5", "3",
                               ifelse(summK3$value == "e7.5", "1", summK3$value)))
plotClus(summK3, "Hierarchical Clustering of Transcripts, k = 4")
summK32 <- dplyr::distinct(summK3[,c("value", "transcripts")])
summK32 <- summK32[order(summK32$transcripts),]
```

For K-means clustering: <br>
```{r}
set.seed(123)
km <- kmeans(tpm2, centers = 4)
trans_cluster_kmeans <- km$cluster %>% enframe()
write.table(trans_cluster_kmeans, "Files/trans_cluster_kmeans_k4.txt", quote = F, row.names = F, sep = "\t")
summKmeans <- summ(trans_cluster_kmeans, tpm2)

plotClus(summKmeans, "K-means Clustering of Transcripts, k = 4")

summKmeans2 <- dplyr::distinct(summKmeans[,c("value", "transcripts")])
summKmeans2 <- summKmeans2[order(summKmeans2$transcripts),]

t <- summKmeans2
for (i in 1:4) {
  for (j in 1:4) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in Kmeans cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

For self-organizing maps: <br>
```{r, eval=F}
library("kohonen")
set.seed(123)
som <- som(as.matrix(tpm2), grid = somgrid(4, 1, "rectangular"))
trans_cluster_som <- som$unit.classif %>% enframe()
trans_cluster_som$name <- row.names(tpm2)
```
Since this run could take time, I saved the results to reproduce the plot here. <br>
```{r}
trans_cluster_som <- read.table("Files/trans_cluster_som_k4.txt", header = T)
summSOM <- summ(trans_cluster_som, tpm2)

plotClus(summSOM, "Self-organizing Maps of Transcripts, k = 4")

summSOM2 <- dplyr::distinct(summSOM[,c("value", "transcripts")])
summSOM2 <- summSOM2[order(summSOM2$transcripts),]

t <- summSOM2
for (i in 1:4) {
  for (j in 1:4) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in SOM cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

For spectral clustering: <br>
```{r, eval=F}
library("RclusTool")
set.seed(123)
sim <- computeGaussianSimilarity(tpm2, 1)
res <- spectralClustering(sim, K=4)
trans_cluster_sc <- data_frame(row.names(tpm2), res$label)
colnames(trans_cluster_sc) <- c("name", "value")
```
Plot of transcripts: <br>
```{r}
trans_cluster_sc <- read.table("Files/trans_cluster_sc_k4.txt", header = T)
summSC <- summ(trans_cluster_sc, tpm2)

plotClus(summSC, "Spectral Clustering of Transcripts")

summSC2 <- dplyr::distinct(summSC[,c("value", "transcripts")])
summSC2 <- summSC2[order(summSC2$transcripts),]

t <- summSC2
for (i in 1:4) {
  for (j in 1:4) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in SC cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

### k = 5:
Hierarchical clustering, k = 5:
```{r}
trans_cluster <- cutree(hc, k = 5) %>% enframe()
write.table(trans_cluster, "Files/trans_cluster_hc_k5.txt", quote = F, row.names = F, sep = "\t")
table(trans_cluster$value)
summK3 <- summ(trans_cluster, tpm2)
summK3$value <- ifelse(summK3$value == "1", "e8.5",
                        ifelse(summK3$value == "2", "e9.5",
                               ifelse(summK3$value == "3", "e7.5", summK3$value)))
summK3$value <- ifelse(summK3$value == "e8.5", "2",
                        ifelse(summK3$value == "e9.5", "3",
                               ifelse(summK3$value == "e7.5", "1", summK3$value)))
plotClus(summK3, "Hierarchical Clustering of Transcripts, k = 5")
summK32 <- dplyr::distinct(summK3[,c("value", "transcripts")])
summK32 <- summK32[order(summK32$transcripts),]
```

For K-means clustering: <br>
```{r}
set.seed(123)
km <- kmeans(tpm2, centers = 5)
trans_cluster_kmeans <- km$cluster %>% enframe()
write.table(trans_cluster_kmeans, "Files/trans_cluster_kmeans_k5.txt", quote = F, row.names = F, sep = "\t")
summKmeans <- summ(trans_cluster_kmeans, tpm2)

plotClus(summKmeans, "K-means Clustering of Transcripts, k = 5")

summKmeans2 <- dplyr::distinct(summKmeans[,c("value", "transcripts")])
summKmeans2 <- summKmeans2[order(summKmeans2$transcripts),]

t <- summKmeans2
for (i in 1:5) {
  for (j in 1:5) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in Kmeans cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

For self-organizing maps: <br>
```{r, eval=F}
library("kohonen")
set.seed(123)
som <- som(as.matrix(tpm2), grid = somgrid(5, 1, "rectangular"))
trans_cluster_som <- som$unit.classif %>% enframe()
trans_cluster_som$name <- row.names(tpm2)
```
Since this run could take time, I saved the results to reproduce the plot here. <br>
```{r}
trans_cluster_som <- read.table("Files/trans_cluster_som_k5.txt", header = T)
summSOM <- summ(trans_cluster_som, tpm2)

plotClus(summSOM, "Self-organizing Maps of Transcripts, k = 5")

summSOM2 <- dplyr::distinct(summSOM[,c("value", "transcripts")])
summSOM2 <- summSOM2[order(summSOM2$transcripts),]

t <- summSOM2
for (i in 1:5) {
  for (j in 1:5) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in SOM cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```

For spectral clustering: <br>
```{r, eval=F}
library("RclusTool")
set.seed(123)
sim <- computeGaussianSimilarity(tpm2, 1)
res <- spectralClustering(sim, K=5)
trans_cluster_sc <- data_frame(row.names(tpm2), res$label)
colnames(trans_cluster_sc) <- c("name", "value")
```
Plot of transcripts: <br>
```{r}
trans_cluster_sc <- read.table("Files/trans_cluster_sc_k5.txt", header = T)
summSC <- summ(trans_cluster_sc, tpm2)

plotClus(summSC, "Spectral Clustering of Transcripts, k = 5")

summSC2 <- dplyr::distinct(summSC[,c("value", "transcripts")])
summSC2 <- summSC2[order(summSC2$transcripts),]

t <- summSC2
for (i in 1:5) {
  for (j in 1:5) {
    q <- length(intersect(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])) #number of observed successes
    print(paste("Percent of hierCluster", i, "in SC cluster", j, "is:", q/length(unique(c(summK32[summK32$value == i, "transcripts"], t[t$value == j, "transcripts"])))*100))
  }
}
```


## 2. Marker counts in hierarchical groups
Markers were collected from 5 review papers (see References). <br>
We first load necessary files: <br>
```{r}
t2g <- read.table("Files/t2g.txt", header = T, sep = "\t")

e7.5specific <- read.table("Files/e7.5specific_ensGenes.txt", header = F)
e7.5hier <- read.table("Files/e7.5Group.txt", header = T)
e8.5specific <- read.table("Files/e8.5specific_ensGenes.txt", header = F)
e8.5hier <- read.table("Files/e8.5Group.txt", header = T)
e9.5specific <- read.table("Files/e9.5specific_ensGenes.txt", header = F)
e9.5hier <- read.table("Files/e9.5Group.txt", header = T)
```
Trophoblast giant cell differentiation markers: <br>
```{R}
tgc <- read.table("Files/0_e7.5_TGCdiffMarkers.txt", header = F)
tgc <- dplyr::inner_join(tgc, t2g[,2:3], by = c("V1" = "ext_gene"))
tgc <- dplyr::distinct(tgc)

tgc[tgc$ens_gene %in% e7.5specific$V1,]
tgc[tgc$ens_gene %in% e7.5hier$ens_gene,]

tgc[tgc$ens_gene %in% e8.5specific$V1,]
tgc[tgc$ens_gene %in% e8.5hier$ens_gene,]
dim(tgc[tgc$ens_gene %in% e8.5hier$ens_gene,])

tgc[tgc$ens_gene %in% e9.5specific$V1,]
tgc[tgc$ens_gene %in% e9.5hier$ens_gene,]
dim(tgc[tgc$ens_gene %in% e9.5hier$ens_gene,])
```
Ectoplacental cone and spongiotrophoblast maintainance markers: <br>
```{r}
epcSpt <- read.table("Files/0_e7.5_EPC-SGT-markers.txt", header = F)
epcSpt <- dplyr::inner_join(epcSpt, t2g[,2:3], by = c("V1" = "ext_gene"))
epcSpt <- dplyr::distinct(epcSpt)

epcSpt[epcSpt$ens_gene %in% e7.5specific$V1,]
epcSpt[epcSpt$ens_gene %in% e7.5hier$ens_gene,]
dim(epcSpt[epcSpt$ens_gene %in% e7.5hier$ens_gene,])

epcSpt[epcSpt$ens_gene %in% e8.5specific$V1,]
epcSpt[epcSpt$ens_gene %in% e8.5hier$ens_gene,]
dim(epcSpt[epcSpt$ens_gene %in% e8.5hier$ens_gene,])

epcSpt[epcSpt$ens_gene %in% e9.5specific$V1,]
epcSpt[epcSpt$ens_gene %in% e9.5hier$ens_gene,]
dim(epcSpt[epcSpt$ens_gene %in% e9.5hier$ens_gene,])
```
Chorioallantoic attachment markers: <br>
```{r}
chorioAll <- read.table("Files/0_e8.5_chorioallantoicAttachment.txt", header = F)
chorioAll <- dplyr::left_join(chorioAll, t2g[,2:3], by = c("V1" = "ext_gene"))
chorioAll <- dplyr::distinct(chorioAll)

chorioAll[chorioAll$ens_gene %in% e7.5specific$V1,]
chorioAll[chorioAll$ens_gene %in% e7.5hier$ens_gene,]
dim(chorioAll[chorioAll$ens_gene %in% e7.5hier$ens_gene,])

chorioAll[chorioAll$ens_gene %in% e8.5specific$V1,]
chorioAll[chorioAll$ens_gene %in% e8.5hier$ens_gene,]
dim(chorioAll[chorioAll$ens_gene %in% e8.5hier$ens_gene,])

chorioAll[chorioAll$ens_gene %in% e9.5specific$V1,]
chorioAll[chorioAll$ens_gene %in% e9.5hier$ens_gene,]
dim(chorioAll[chorioAll$ens_gene %in% e9.5hier$ens_gene,])
```
Labyrinth branching and vascularization - Syncytiotrophoblast markers: <br>
```{r}
laby <- read.table("Files/0_e9.5_branching-labyrinthMarkers.txt", header = F)
laby <- dplyr::left_join(laby, t2g[,2:3], by = c("V1" = "ext_gene"))
laby <- dplyr::distinct(laby)

laby[laby$ens_gene %in% e7.5specific$V1,]
head(laby[laby$ens_gene %in% e7.5hier$ens_gene,])
dim(laby[laby$ens_gene %in% e7.5hier$ens_gene,])

laby[laby$ens_gene %in% e8.5specific$V1,]
head(laby[laby$ens_gene %in% e8.5hier$ens_gene,])
dim(laby[laby$ens_gene %in% e8.5hier$ens_gene,])

laby[laby$ens_gene %in% e9.5specific$V1,]
head(laby[laby$ens_gene %in% e9.5hier$ens_gene,])
dim(laby[laby$ens_gene %in% e9.5hier$ens_gene,])
```

## References
1. E. D. Watson and J. C. Cross, “Development of structures and transport functions in the mouse placenta,” Physiology, vol. 20, no. 3, pp. 180–193, 2005. <br>
2. J. C. Cross, “How to Make a Placenta: Mechanisms of Trophoblast Cell Differentiation in Mice – A Review,” vol. 26, 2005. <br>
3. M. Hemberger and J. C. Cross, “Genes governing placental development,” vol. 12, no. 4, pp. 162–168, 2001. <br>
4. D. Hu and J. C. Cross, “Development and function of trophoblast giant cells in the rodent placenta,” Int. J. Dev. Biol., vol. 54, no. 2–3, pp. 341–354, 2010. <br>
5. J. Rossant, J. C. Cross, and S. Lunenfeld, “Placental Development: Lessons from Mouse Mutants,” vol. 2, no. July, pp. 538–548, 2001. <br>

```{r}
sessionInfo()
```