---
title: "Connexinplexity_2021"
author: "Rachel Lukowicz-Bedford"
date: "10/20/2021"
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
rm(list=ls()) 
library(ggridges)
#install.packages('Seurat')
library(Seurat)
library(Matrix)
library(sctransform)
library(stringr)
library(BiocManager)
library(ggplot2)
library(RColorBrewer)
```

All pre-processing was done before this. 

```{r, import data}
all  <- readRDS('/Users/rachellukowicz/Desktop/Connexplexity/v3/re-run_updatedGTF/Connexinplexity_Final.rds')
DimPlot(all, label = T, label.size = 1.5) + NoLegend()
#original  <- readRDS('/Users/rachellukowicz/Desktop/Connexplexity/v3/Farnsworth_2019.rds')
```

Figure 1 - Spatial information - plots of the dataset with labels and highlight Cx. 
```{r, fig1_spatial}
DimPlot(all, label = F) + NoLegend()
DimPlot(all, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93")) 
FeaturePlot(all, label = F, features = c("gjc4b"), order = T, cols = c("gray93", "#FE4A49"))
FeaturePlot(all, label = F, features = c("gja8b"), order = T, cols = c("gray93", "#FE4A49"))
FeaturePlot(all, label = F, features = c("gja1b"), order = T, cols = c("gray93", "#FE4A49"))
```

Figure  1 - dotplots for all Cx, but seperate the broadly expressed ones. 
```{r, fig1_Dotplots, fig.height = 2, fig.width = 7}
DotPlot(all, features = c("gja1a",  "gja2", "gja3", "gja4", "gja5a", "gja5b", "gja8a", "gja8b", "gja9a", "gja9b", "gja10a", "gja10b", "gja11", "gja12.1", "gja12.2", "gja13.2", "gja13.1", "gjb1a", "gjb1b", "gjb3", "gjb7", "gjb8", "gjb9a", "gjb9b", "gjb10", "gjc1", "gjc2", "gjc4a.1", "gjc4a.2", "gjd1a", "gjd1b", "gjd2a", "gjd2b", "gjd4", "gjd5", "gjd6", "gje1b", "gje1a", "gjz1"), scale.by = "size", group.by = "orig.ident", cols = c("gray93", "#FE4A49")) +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 15)) +
  scale_x_discrete(limits=c("gja1a", "gja2", "gja3", "gja4", "gja5a", "gja5b", "gja8a", "gja8b", "gja9a", "gja9b", "gja10a", "gja10b", "gja11", "gja12.1", "gja12.2", "gja13.2", "gja13.1", "gjb1a", "gjb1b", "gjb3", "gjb7", "gjb8", "gjb9a", "gjb9b", "gjb10", "gjc1", "gjc2", "gjc4a.1", "gjc4a.2", "gjd1a", "gjd1b", "gjd2a", "gjd2b", "gjd4", "gjd5", "gjd6", "gje1b", "gje1a", "gjz1"))  + scale_y_discrete(limits=c("5dpf", "2dpf", "1dpf"))

DotPlot(all, features = c("gjc4b", "gja1b"), scale.by = "size", group.by = "orig.ident", cols = c("gray93", "#FE4A49")) + scale_y_discrete(limits=c("5dpf", "2dpf", "1dpf"))

```

Figure 2- dot plot of all Cx over all clusters. Cx arranged by expression, and clusters orgnaized by germ layer and tissue. 
```{r, fig2_dotplot, fig.height = 13, fig.width = 10}
DotPlot(all, features = c("gjc4b", "gja1b", "gje1a", "gja8a", "gja8b", "gja3", "gjd5", "gjd2b", "gjd1b", "gja9a",  "gja10b", "gjd2a", "gjd1a",  "gjb1a", "gjb1b", "gjc2", "gja4", "gja5b",  "gjb3", "gjb8", "gjb9a",  "gjb10", "gjc4a.1", "gjc4a.2", "gjb7", "gjb9b", "gja11", "gja12.1", "gja13.2", "gja1a","gja2", "gje1b", "gjd4",  "gja5a", "gja10a", "gjc1", "gjd6", "gja13.1", "gja12.2", "gjz1", "gja9b"), scale.by = "size", cols =c("gray93", "#FE4A49")) +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 15)) +
  scale_x_discrete(limits=c("gjc4b", "gja1b", "gjz1", "gje1a", "gja8a", "gja8b", "gja3", "gjd5", "gjd1b", "gjd2b",   "gja10b", "gja9b", "gjd2a", "gjd1a",  "gjb1a", "gjb1b", "gjc2", "gja1a", "gja4", "gja5b",  "gjb3", "gjb10", "gjb8", "gjb9a",   "gjc4a.1", "gjc4a.2", "gjb7", "gjb9b", "gja11", "gja12.1", "gja13.1", "gja2", "gje1b", "gjd4",  "gja5a", "gja10a", "gjc1", "gjd6", "gja13.2", "gja12.2", "gja9a" )) + 
  scale_y_discrete(limits=c("156", "138", "136", "135", "134", "133", "132", "131", "131", "128", "128", "59", "59", "21", "22", "22", "117", "167", "66", "65", "64", "64", "226", "225", "120", "120", "230", "229", "228", "227", "227", "18", "18", "20", "19", "17", "17", "72", "72",  "232", "165", "101", "101", "231", "71", "6", "6", "234", "5", "4", "233", "11", "10", "9", "8", "7", "7", "83", "119", "84", "123", "114", "113", "118", "91", "92", "88", "87","87", "63", "63", "168", "89", "70", "70", "236", "235", "235", "116", "115", "58", "58", "68", "69", "93", "93", "82", "81", "81", "224", "223", "222", "222", "16", "219", "219", "2", "1", "1", "140", "67", "60", "3","3", "75", "75", "15", "15", "94", "94", "49", "46", "45", "44", "43", "42", "41", "40", "221", "33", "48", "39", "38", "47", "37", "36", "35", "34", "27", "25", "23", "51", "32", "31", "30", "29", "28", "26", "24", "24", "53", "52", "57", "56", "54", "90", "169", "166", "163", "130", "122", "99", "98", "97", "96", "95", "62", "61", "86", "85", "164", "162", "76", "76", "50", "100", "121", "121", "173", "145", "145", "74", "177","177","172", "171", "161", "190", "189", "182", "127", "197", "183", "195", "181", "176", "154", "152", "151", "150", "148", "207", "218", "188", "147", "186", "160", "184", "194", "192", "187", "175", "174", "174", "178", "170", "139", "105", "196", "180", "179", "153", "149", "104", "103", "212", "204", "191", "205", "199", "125", "124", "143", "146","146", "144", "198", "111", "73", "211", "209", "208", "155", "142", "141", "110", "108", "106", "107", "238", "202", "109", "102", "220", "220", "129", "137", "157", "158", "159", "185", "78", "79", "80", "126", "200", "201", "213", "214", "215", "216", "217", "77", "203", "206", "210", "12", "13", "14", "14", "237", "193", "112", "112")) + theme(axis.text.y = element_text(angle = 45, hjust=1, size = 5))

DotPlot(all, features = c("gjc4b", "gja1b", "gje1a", "gja8a", "gja8b", "gja3", "gjd5", "gjd2b", "gjd1b", "gja9a",  "gja10b", "gjd2a", "gjd1a",  "gjb1a", "gjb1b", "gjc2", "gja4", "gja5b",  "gjb3", "gjb8", "gjb9a",  "gjb10", "gjc4a.1", "gjc4a.2", "gjb7", "gjb9b", "gja11", "gja12.1", "gja13.2", "gja1a","gja2", "gje1b", "gjd4",  "gja5a", "gja10a", "gjc1", "gjd6", "gja13.1", "gja12.2", "gjz1", "gja9b"), scale.by = "size", cols =c("gray93", "#FE4A49")) +  theme(axis.text.x = element_text(angle = 45, hjust=1, size = 15)) +
  scale_x_discrete(limits=c("gjc4b", "gja1b", "gjz1", "gje1a", "gja8a", "gja8b", "gja3", "gjd5", "gjd1b", "gjd2b",   "gja10b", "gja9b", "gjd2a", "gjd1a",  "gjb1a", "gjb1b", "gjc2", "gja1a", "gja4", "gja5b",  "gjb3", "gjb10", "gjb8", "gjb9a",   "gjc4a.1", "gjc4a.2", "gjb7", "gjb9b", "gja11", "gja12.1", "gja13.1", "gja2", "gje1b", "gjd4",  "gja5a", "gja10a", "gjc1", "gjd6", "gja13.2", "gja12.2", "gja9a" )) + 
  scale_y_discrete(limits=c("156", "138", "136", "135", "134", "133", "132", "131", "131", "128", "128", "59", "59", "21", "22", "22", "117", "167", "66", "65", "64", "64", "226", "225", "120", "120", "230", "229", "228", "227", "227", "18", "18", "20", "19", "17", "17", "72", "72",  "232", "165", "101", "101", "231", "71", "6", "6", "234", "5", "4", "233", "11", "10", "9", "8", "7", "7", "83", "119", "84", "123", "114", "113", "118", "91", "92", "88", "87","87", "63", "63", "168", "89", "70", "70", "236", "235", "235", "68", "69", "93", "93", "116", "115", "58", "58", "82", "81", "81", "224", "223", "222", "222", "16", "219", "219", "2", "1", "1", "140", "67", "60", "3","3", "75", "75", "15", "15", "94", "94", "50", "100", "121", "121", "49", "46", "45", "44", "43", "42", "41", "40", "221", "33", "48", "39", "38", "47", "37", "36", "35", "34", "27", "25", "23", "51", "32", "31", "30", "29", "28", "26", "24", "24", "53", "52", "57", "56", "54", "90", "169", "166", "163", "130", "122", "99", "98", "97", "96", "95", "62", "61", "86", "85", "164", "162", "76", "76", "173", "145", "145", "74", "177","177","172", "171", "161", "190", "189", "182", "127", "197", "183", "195", "181", "176", "154", "152", "151", "150", "148", "207", "218", "188", "147", "186", "160", "184", "194", "192", "187", "175", "174", "174", "178", "170", "139", "105", "196", "180", "179", "153", "149", "104", "103", "212", "204", "191", "205", "199", "125", "124", "143", "146","146", "144", "198", "111", "73", "211", "209", "208", "155", "142", "141", "110", "108", "106", "107", "238", "202", "109", "102", "220", "220", "129", "137", "157", "158", "159", "185", "78", "79", "80", "126", "200", "201", "213", "214", "215", "216", "217", "77", "203", "206", "210", "12", "13", "14", "14", "237", "193", "112", "112")) + theme(axis.text.y = element_text(angle = 45, hjust=1, size = 5))

```

Figure 3 - integument umaps for marker genes and cx 
```{r, fig3_Umaps}
skin_sub <- subset(all, idents = c("24", "26", "28", "29", "30", "31", "32", "51", "39", "38", "40", "41", "42", "43", "44", "45", "46", "47", "48", "54", "52", "53", "56", "57", "23", "25", "27", "222", "223", "224"))

DimPlot(object = skin_sub, reduction = 'umap', cols = c('222' = '#BF6C78', '223' = '#BF6C78', '224' = '#BF6C78', '23' = '#BF6C78', '25' = '#BF6C78', '27' = '#BF6C78', '24' = '#BF6C78', '26' = '#BF6C78', '28' = '#BF6C78', '29' = '#BF6C78', '30' = '#BF6C78', '31' = '#BF6C78', '32' = '#BF6C78', '51' = '#BF6C78', '52' = '#61A899', '53' = '#61A899', '54' = '#61A899', '56' = '#61A899', '57' = '#61A899', '39' = '#999944', '38' = '#999944', '47' = '#999944', '48' = '#999944', '40' = '#9E4A95', '41' = '#9E4A95', '42' = '#9E4A95', '43' = '#9E4A95', '44' = '#9E4A95', '45' = '#9E4A95', '46' = '#9E4A95'), label = F, label.size = 3, repel = F)

#Marker genes 
FeaturePlot(skin_sub, features = c("ppl"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("foxi3a"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("tp63"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("sox10"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))

#Cx Clusters 
FeaturePlot(skin_sub, features = c("gjb7"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjb8"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjb9a"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjc4a.2"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjc4a.1"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjb3"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjb10"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjc4b"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gja5b"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gjb9b"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(skin_sub, features = c("gja4"), label = F, label.size = 3, repel = T, order = T, cols =c("gray93", "#FE4A49"))
```

Figure 3 - gjb8 (periderm and pigment), and gjb7 (ionocytes)
```{r, fig3_plots}
pig_sub <- subset(all, idents = c("52", "53", "54", "56", "57"))
peri_sub <- subset(all, idents = c("40", "41", "42","43","44","45", "46"))
iono_sub <- subset(all, idents = c("39", "38", "47", "48"))

                                                                                                           
#########pigment cells
p1 <-  VlnPlot(pig_sub, features = c("gjb8"), cols = c("#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB"), pt.size = 0)
p1 + geom_jitter(mapping = aes(color = '#7E7877')) +
    scale_color_manual(values = c('#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877')) + NoLegend()

p2 <-  VlnPlot(pig_sub, features = c("aox5"), cols = c("#61A899", "#61A899", "#61A899", "#61A899", "#61A899", "#61A899", "#61A899"), pt.size = 0)
p2 + geom_jitter(mapping = aes(color = '#488477')) +
    scale_color_manual(values = c('#488477', '#488477', '#488477', '#488477', '#488477', '#488477', '#488477')) + NoLegend()

p3 <-  VlnPlot(pig_sub, features = c("dct"), cols = c("#3A695F", "#3A695F", "#3A695F", "#3A695F", "#3A695F", "#3A695F", "#3A695F"), pt.size = 0)
p3 + geom_jitter(mapping = aes(color = '#61A899')) +
    scale_color_manual(values = c('#61A899', '#61A899', '#61A899', '#61A899', '#61A899', '#61A899', '#61A899')) + NoLegend()

#####periderm 
p <-  VlnPlot(peri_sub, features = c("gjb8"), cols = c("#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB", "#B0ACAB"), pt.size = 0)
p + geom_jitter(mapping = aes(color = '#7E7877'), data = p$data) +
    scale_color_manual(values = c('#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877', '#7E7877')) + NoLegend()

p <-  VlnPlot(peri_sub, features = c("ppl"), cols = c("#9E4A95", "#9E4A95", "#9E4A95", "#9E4A95", "#9E4A95", "#9E4A95", "#9E4A95"), pt.size = 0)
p + geom_jitter(mapping = aes(color = '#6F3469'), data = p$data) +
    scale_color_manual(values = c('#6F3469', '#6F3469', '#6F3469', '#6F3469', '#6F3469', '#6F3469', '#6F3469')) + NoLegend()

#########ionocytes
p <-  VlnPlot(iono_sub, feature = "gjb7", cols = c("#7E7877", "#7E7877", "#7E7877", "#7E7877"), pt.size = 0)
p + geom_jitter(mapping = aes(color = '#7E7877'), data = p$data) +
    scale_color_manual(values = c('#7E7877', '#7E7877', '#7E7877', '#7E7877')) + 
  NoLegend()

p <-  VlnPlot(iono_sub, feature = "atp1b1b", cols = c("#F4CB15", "#F4CB15", "#F4CB15", "#F4CB15"), pt.size = 0)
p + geom_jitter(mapping = aes(color = '#C3A109'), data = p$data) +
    scale_color_manual(values = c('#C3A109', '#C3A109', '#C3A109', '#C3A109')) + 
  NoLegend()

p <-  VlnPlot(iono_sub, feature = "atp6v1aa", cols = c("#00B800", "#00B800", "#00B800", "#00B800"), pt.size = 0)
p + geom_jitter(mapping = aes(color = '#007A00'), data = p$data) +
    scale_color_manual(values = c('#007A00', '#007A00', '#007A00', '#007A00')) + 
  NoLegend()
```
Supplemental figure 1 was made in geneious, figure 2 also made in genious. 

Import Farnsworth dataset first. 

```{r, suppfig_3}
original  <- readRDS('/Users/rachellukowicz/Desktop/Connexplexity/v3/Farnsworth_2019.rds')
library(Seurat)
original <- UpdateSeuratObject(original)

#FeaturePlot(original, features = c("GJD4"), label = T, label.size = 2, order = T)
#FeaturePlot(all, features = c("gjd4"), label = T, label.size = 2, order = T)
original_slowmuscle <- subset(original, idents = ("57"))
new_slowmuscle <- subset(all, idents = ("4"))
VlnPlot(original_slowmuscle, features = c("smyhc2", "GJD4"))
VlnPlot(new_slowmuscle, features = c("smyhc2", "gjd4"))

original_schwann <- subset(original, idents = ("173"))
new_schwann <- subset(all, idents = ("76"))
VlnPlot(original_schwann, features = c("mbpa", "cx27.5", "cx31.7", "cx47.1"), cols = c("skyblue3"))
VlnPlot(new_schwann, features = c("mbpa", "gjb1a", "gjb1b", "gjc2"), cols = c("skyblue3"))
VlnPlot(original_schwann, features = c("mbpa", "cx27.5"), cols = c("skyblue3"))
VlnPlot(new_schwann, features = c("mbpa", "gjb1a"), cols = c("skyblue3"))

original_heart <- subset(original, idents = ("205"))
new_heart <- subset(all, idents = ("232"))
VlnPlot(original_heart, features = c("myl7", "cx36.7"), cols = c("olivedrab3"))
VlnPlot(new_heart, features = c("myl7", "gjd6"), cols = c("olivedrab3"))

```

Supp fig 4- all Cx in umap space
```{r, supp_fig4_umap}
FeaturePlot(all, features = c("gja1a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja1b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja2"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja3"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja4"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja5a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja5b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja8a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja8b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja9a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja9b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja10a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja10b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja11"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja12.1"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja12.2"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja13.1"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gja13.2"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb1a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb1b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb3"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb7"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb8"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb9a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb9b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjb10"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjc1"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjc2"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjc4a.1"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjc4a.2"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjc4b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd1a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd1b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd2a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd2b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd4"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd5"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjd6"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gje1a"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gje1b"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)
FeaturePlot(all, features = c("gjz1"), order = T, cols = c("gray93", "#FE4A49"), label = T, label.size = 1.5)

```

continuation of supp figure 4 - violin plots for all the connexin
```{r, supp_fig4_vln, fig.height = 3, fig.width = 20}
library(ggplot2)
VlnPlot(all, features = c("gja1a")) + NoLegend() 
VlnPlot(all, features = c("gja1b")) + NoLegend() 
VlnPlot(all, features = c("gja2")) + NoLegend()
VlnPlot(all, features = c("gja3")) + NoLegend()
VlnPlot(all, features = c("gja4")) + NoLegend() 
VlnPlot(all, features = c("gja5a")) + NoLegend()
VlnPlot(all, features = c("gja5b")) + NoLegend()
VlnPlot(all, features = c("gja8a")) + NoLegend()
VlnPlot(all, features = c("gja8b")) + NoLegend()
VlnPlot(all, features = c("gja9a")) + NoLegend() 
VlnPlot(all, features = c("gja9b")) + NoLegend() 
VlnPlot(all, features = c("gja10a")) + NoLegend() 
VlnPlot(all, features = c("gja10b")) + NoLegend()
VlnPlot(all, features = c("gja11")) + NoLegend()
VlnPlot(all, features = c("gja12.1")) + NoLegend()
VlnPlot(all, features = c("gja12.2")) + NoLegend()
VlnPlot(all, features = c("gja13.1")) + NoLegend()
VlnPlot(all, features = c("gja13.2")) + NoLegend()
VlnPlot(all, features = c("gjb1a")) + NoLegend()
VlnPlot(all, features = c("gjb1b")) + NoLegend()
VlnPlot(all, features = c("gjb3")) + NoLegend()
VlnPlot(all, features = c("gjb7")) + NoLegend()
VlnPlot(all, features = c("gjb8")) + NoLegend()
VlnPlot(all, features = c("gjb9a")) + NoLegend()
VlnPlot(all, features = c("gjb9b")) + NoLegend()
VlnPlot(all, features = c("gjb10")) + NoLegend()
VlnPlot(all, features = c("gjc1")) + NoLegend()
VlnPlot(all, features = c("gjc2")) + NoLegend()
VlnPlot(all, features = c("gjc4a.1")) + NoLegend()
VlnPlot(all, features = c("gjc4a.2")) + NoLegend()
VlnPlot(all, features = c("gjc4b")) + NoLegend()
VlnPlot(all, features = c("gjd1a")) + NoLegend()
VlnPlot(all, features = c("gjd1b")) + NoLegend()
VlnPlot(all, features = c("gjd2a")) + NoLegend()
VlnPlot(all, features = c("gjd2b")) + NoLegend()
VlnPlot(all, features = c("gjd4")) + NoLegend()
VlnPlot(all, features = c("gjd5")) + NoLegend()
VlnPlot(all, features = c("gjd6")) + NoLegend()
VlnPlot(all, features = c("gje1a")) + NoLegend()
VlnPlot(all, features = c("gje1b")) + NoLegend()
VlnPlot(all, features = c("gjz1")) + NoLegend()
```

Supp figure 5 - tissue markers with connexin expression
```{r, supp_fig5}
lens <- subset(all, idents = c("237", "193", "112"))
FeaturePlot(lens, features = c("cryba2b"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(lens, features = c("crybb1"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(lens, features = c("gja8b"), order = T, label = F, cols =c("gray93", "#FE4A49"))
DimPlot(lens, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))

heart <- subset(all, idents = c("232"))
FeaturePlot(heart, features = c("myl7"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(heart, features = c("nppa"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(heart, features = c("gjd6"), order = T, label = F, cols =c("gray93", "#FE4A49"))
DimPlot(heart, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))


muscle <- subset(all, idents = c("234", "5", "4", "233", "11", "10", "9", "8", "7"))
FeaturePlot(muscle, features = c("smyhc1"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(muscle, features = c("myhz2"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(muscle, features = c("gja2"), order = T, label = F, cols =c("gray93", "#FE4A49"))
DimPlot(muscle, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))

horizontal <- subset(all, idents = c("203"))
FeaturePlot(horizontal, features = c("mdka"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(horizontal, features = c("gja9b"), order = T, label = F, cols =c("gray93", "#FE4A49"))
FeaturePlot(horizontal, features = c("gja10b"), order = T, label = F, cols =c("gray93", "#FE4A49"))
DimPlot(horizontal, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))

neuro <- subset(all, idents = c("178", "170", "139", "105", "196", "180", "179", "153", "149", "104", "103", "212", "204", "191", "205", "199", "125", "124",  "143", "146", "144", "198", "111", "73", "211", "209", "208", "155", "142", "141", "110", "108", "106", "107", "238", "202",  "109", "102", "220", "172", "171", "161", "190", "189", "182", "127", "197", "183", "195", "181", "176", "154", "152", "151", "150", "148", "207", "218", "188", "147", "186", "160", "184", "194", "192", "187", "175", "174", "173", "145"))

FeaturePlot(neuro, features = c("pcna"), order = F, label = F, label.size = 2, repel = F,  cols =c("gray93", "#FE4A49"))
FeaturePlot(neuro, features = c("snap25a"), order = F,  label = F, label.size = 2,  cols =c("gray93", "#FE4A49"))
FeaturePlot(neuro, features = c("gjc4b"), order = F,  label = F, label.size = 2, cols =c("gray93", "#FE4A49"))
DimPlot(neuro, group.by = "orig.ident",  cols = c("#82A9D9", "#E08D6D", "#EBDE93"), label = T)

```

Supp figure 6 - intestine over time 
```{r, supp_fig6}
intestine <- subset(all, idents = c("140", "67", "60", "3"))
intestine_group <- RenameIdents(intestine, "140" = "Intestine", "67" = "Intestine", "60" = "Intestine", "3" = "Intestine")
VlnPlot(intestine_group, features = c("cldnc", "fabp2", "foxa3"))
VlnPlot(intestine, features = c("gjc4b", "gja12.1", "gja13.1"), group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))
```

Supp fig 7 - PGC and connnexin expression 
```{r, supp_fig7}
pgc <- subset(all, idents = c("128"))
FeaturePlot(pgc, features = c("ddx4"), order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(pgc, features = c("nanos3"), order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(pgc, features = c("gjb8"), order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(pgc, features = c("gja9a"), order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(pgc, features = c("gjd1b"), order = T, cols =c("gray93", "#FE4A49"))
FeaturePlot(pgc, features = c("gjc4b"), order = T, cols =c("gray93", "#FE4A49"))
DimPlot(pgc, group.by = "orig.ident", cols = c("#82A9D9", "#E08D6D", "#EBDE93"))

```

Supp figure 8 - marker genes and connexin expression within integument cells 
```{r, supp_Fig8, fig.height = 11, fig.width = 8.5}
skin_sub_annotated <- RenameIdents(skin_sub, "222" = "Basal Cell", "223" = "Basal Cell", "224" = "Basal Cell", "23" = "Basal Cell", "25" = "Basal Cell", "27" = "Basal Cell", "24" = "Basal Cell", "26" = "Basal Cell", "28" = "Basal Cell", "29" = "Basal Cell", "30" = "Basal Cell", "31" = "Basal Cell", "32" = "Basal Cell", "51" = "Basal Cell", "52" = "Pigment Cell", "53" = "Pigment Cell", "54" = "Pigment Cell", "56" = "Pigment Cell", "57" = "Pigment Cell", "39" = "Ionocyte", "38" = "Ionocyte", "47" = "Ionocyte", "48" = "Ionocyte", "40" = "Periderm", "41" = "Periderm", "42" = "Periderm", "43" = "Periderm", "44" = "Periderm", "45" = "Periderm", "46" = "Periderm")
VlnPlot(skin_sub_annotated, features = c("ppl", "gjb9a","ppl", "ppl", "foxi3a", "gjb7", "gjb9b", "foxi3b", "sox10", "aox5", "gja5b", "gja4", "tp63", "gjc4a.2", "gjc4a.1", "gjb8", "gjb3", "gjb10", "gjc4b"), combine = TRUE)
```


```{r, groupings for dimplots}
play <- RenameIdents(all, "178" = "Differentiating Neuron", "170" = "Differentiating Neuron", "139" = "Differentiating Neuron", "105" = "Differentiating Neuron", "196" = "Differentiating Neuron", "180" = "Differentiating Neuron", "179" = "Differentiating Neuron", "153" = "Differentiating Neuron", "149" = "Differentiating Neuron", "104" = "Differentiating Neuron", "103" = "Differentiating Neuron", "212" = "Differentiating Neuron", "204" = "Differentiating Neuron", "191" = "Differentiating Neuron", "205" = "Differentiating Neuron", "199" = "Differentiating Neuron", "125" = "Differentiating Neuron", "124" = "Differentiating Neuron", "143" = "Differentiating Neuron", "146" = "Differentiating Neuron")

play <- RenameIdents(play, "156" = "Non-Specific", "138" = "Non-Specific", "136" = "Non-Specific", "135" = "Non-Specific", "134" = "Non-Specific", "133" = "Non-Specific", "132" = "Non-Specific", "131" = "Non-Specific")

play <- RenameIdents(play, "226" = "Endothelial", "225" = "Endothelial", "120" = "Endothelial")

play <- RenameIdents(play, "230" = "Blood", "229" = "Blood", "228" = "Blood", "227" = "Blood")

play <- RenameIdents(play, "20" = "Macrophage", "19" = "Macrophage", "17" = "Macrophage")

play <- RenameIdents(play, "72" = "Vasculature")

play <- RenameIdents(play, "232" = "Heart", "165" = "Heart", "101" = "Heart")

play <- RenameIdents(play, "144" = "Neuron", "198" = "Neuron", "111" = "Neuron", "73" = "Neuron", "211" = "Neuron", "209" = "Neuron", "208" = "Neuron", "155" = "Neuron", "142" = "Neuron", "141" = "Neuron", "110" = "Neuron", "108" = "Neuron", "106" = "Neuron", "107" = "Neuron", "238" = "Neuron", "202" = "Neuron", "109" = "Neuron", "102" = "Neuron", "220" = "Neuron")

play <- RenameIdents(play, "18" = "Neutrophil")

play <- RenameIdents(play, "128" = "Germ Cell")

play <- RenameIdents(play, "59" = "Kidney")

play <- RenameIdents(play, "21" = "Spleen", "22" = "Spleen")

play <- RenameIdents(play, "167" = "Pectoral Fin Bud", "66" = "Pectoral Fin Bud", "65" = "Pectoral Fin Bud", "64" = "Pectoral Fin Bud")

play <- RenameIdents(play, "234" = "Skeletal Muscle", "5" = "Skeletal Muscle", "4" = "Skeletal Muscle", "233" = "Skeletal Muscle", "11" = "Skeletal Muscle", "10" = "Skeletal Muscle", "9" = "Skeletal Muscle", "8" = "Skeletal Muscle", "7" = "Skeletal Muscle")

play <- RenameIdents(play, "224" = "Pharyngeal Endoderm", "223" = "Pharyngeal Endoderm", "222" = "Pharyngeal Endoderm")

play <- RenameIdents(play, "116" = "Hatching Gland", "115" = "Hatching Gland", "58" = "Hatching Gland")

play <- RenameIdents(play, "68" = "Mesenchyme-Related", "69" = "Mesenchyme-Related", "93" = "Mesenchyme-Related")

play <- RenameIdents(play, "82" = "Thymus", "81" = "Thymus")

play <- RenameIdents(play, "50" = "Placode", "100" = "Placode", "121" = "Placode")

play <- RenameIdents(play, "83" = "Paraxial Mesoderm", "119" = "Paraxial Mesoderm", "84" = "Paraxial Mesoderm", "123" = "Paraxial Mesoderm", "114" = "Paraxial Mesoderm", "113" = "Paraxial Mesoderm", "118" = "Paraxial Mesoderm", "91" = "Paraxial Mesoderm", "92" = "Paraxial Mesoderm", "88" = "Paraxial Mesoderm", "87" = "Paraxial Mesoderm")

play <- RenameIdents(play, "16" = "Pancreas", "219" = "Pancreas")

play <- RenameIdents(play, "1" = "Liver", "2" = "Liver")

play <- RenameIdents(play, "55" = "Retina", "14" = "Retina", "13" = "Retina", "126" = "Retina", "217" = "Retina", "216" = "Retina", "215" = "Retina", "214" = "Retina", "213" = "Retina", "201" = "Retina", "200" = "Retina", "80" = "Retina", "79" = "Retina", "78" = "Retina", "12" = "Retina", "210" = "Retina", "206" = "Retina", "203" = "Retina", "77" = "Retina", "185" = "Retina", "159" = "Retina", "157" = "Retina", "137" = "Retina", "158" = "Retina", "129" = "Retina")

play <- RenameIdents(play, "140" = "Intestine", "67" = "Intestine", "60" = "Intestine", "3" = "Intestine")

play <- RenameIdents(play, "75" = "Hypochord")

play <- RenameIdents(play, "15" = "Epiphysis")

play <- RenameIdents(play, "94" = "Hypophysis")

play <- RenameIdents(play, "172" = "Neuroblast", "171" = "Neuroblast", "161" = "Neuroblast", "190" = "Neuroblast", "189" = "Neuroblast", "182" = "Neuroblast", "127" = "Neuroblast", "197" = "Neuroblast", "183" = "Neuroblast", "195" = "Neuroblast", "181" = "Neuroblast", "176" = "Neuroblast", "154" = "Neuroblast", "152" = "Neuroblast", "151" = "Neuroblast", "150" = "Neuroblast", "148" = "Neuroblast", "207" = "Neuroblast", "218" = "Neuroblast", "188" = "Neuroblast", "147" = "Neuroblast", "186" = "Neuroblast", "160" = "Neuroblast", "184" = "Neuroblast", "194" = "Neuroblast", "192" = "Neuroblast", "187" = "Neuroblast", "175" = "Neuroblast", "174" = "Neuroblast")

play <- RenameIdents(play, "63" = "Osteoblast")

play <- RenameIdents(play, "168" = "Pharyngeal Arch", "89" = "Pharyngeal Arch", "70" = "Pharyngeal Arch")

play <- RenameIdents(play, "236" = "Notochord", "235" = "Notochord")

play <- RenameIdents(play, "231" = "Muscle", "71" = "Muscle", "6" = "Muscle")

play <- RenameIdents(play, "53" = "Neural Crest", "52" = "Neural Crest", "57" = "Neural Crest", "56" = "Neural Crest", "54" = "Neural Crest", "90" = "Neural Crest", "169" = "Neural Crest", "166" = "Neural Crest", "163" = "Neural Crest", "130" = "Neural Crest","122" = "Neural Crest", "99" = "Neural Crest", "98" = "Neural Crest", "97" = "Neural Crest", "96" = "Neural Crest", "95" = "Neural Crest", "62" = "Neural Crest", "61" = "Neural Crest", "86" = "Neural Crest", "85" = "Neural Crest", "164" = "Neural Crest", "162" = "Neural Crest", "76" = "Neural Crest")

play <- RenameIdents(play, "173" = "Oligodendrocyte", "145" = "Oligodendrocyte")
play <- RenameIdents(play, "117" = "Tailbud")
play <- RenameIdents(play, "237" = "Lens", "193" = "Lens", "112" = "Lens")

play <- RenameIdents(play, "74" = "FloorPlate", "177" = "FloorPlate")

play <- RenameIdents(play, "49" = "Integument", "46" = "Integument", "45" = "Integument", "44" = "Integument", "43" = "Integument", "42" = "Integument", "41" = "Integument", "40" = "Integument", "221" = "Integument", "33" = "Integument", "48" = "Integument", "39" = "Integument", "38" = "Integument", "47" = "Integument", "37" = "Integument", "36" = "Integument", "35" = "Integument", "34" = "Integument", "27" = "Integument", "25" = "Integument", "23" = "Integument", "51" = "Integument", "32" = "Integument", "31" = "Integument", "30" = "Integument", "29" = "Integument", "28" = "Integument", "26" = "Integument", "24" = "Integument")

DimPlot(play, label = T, label.size = 1.5)
```

Supp table 3 - pulling the % of cells that express given connnexin 
```{r}
genes <- c("gjc4b", "gja1b", "gje1a", "gja8a", "gja8b", "gja3", "gjd5", "gjd1b", "gjd2b",   "gja10b", "gja9b", "gjd2a", "gjd1a",  "gjb1a", "gjb1b", "gjc2", "gja1a", "gja4", "gja5b",  "gjb3", "gjb10", "gjb8", "gjb9a",   "gjc4a.1", "gjc4a.2", "gjb7", "gjb9b", "gja11", "gja12.1", "gja13.1", "gja2", "gje1b", "gjd4",  "gja5a", "gja10a", "gjc1", "gjd6", "gja13.2", "gja12.2", "gja9a", "gjz1")

PrctCellExpringGene <- function(object, genes, group.by = "all"){
    if(group.by == "all"){
        prct = unlist(lapply(genes,calc_helper, object=object))
        result = data.frame(Markers = genes, Cell_proportion = prct)
        return(result)
    }

    else{        
        list = SplitObject(object, group.by)
        factors = names(list)

        results = lapply(list, PrctCellExpringGene, genes=genes)
        for(i in 1:length(factors)){
        results[[i]]$Feature = factors[i]
        }
        combined = do.call("rbind", results)
        return(combined)
    }
}

calc_helper <- function(object,genes){
    counts = object[['RNA']]@counts
    ncells = ncol(counts)
    if(genes %in% row.names(counts)){
    sum(counts[genes,]>0)/ncells
    }else{return(NA)}
}

exprinfo <- PrctCellExpringGene(all, genes, group.by = "tree.ident")

#write.csv(exprinfo, 'CxExpressionList_updated.csv')
```

