---
title: "PUBLISHED DATA PREPARATION"
author: "Matthew J Shannon"
date: "12/08/2022"
output: html_document
---

# Load in all dependencies

```{r, include = FALSE}
#devtools::install_github  ('chris-mcginnis-ucsf/DoubletFinder')
#devtools::install_github  ('satijalab/seurat-wrappers')
#remotes::install_github   ("mojaveazure/seurat-disk")
#devtools::install_github  ("velocyto-team/velocyto.R")
#BiocManager::install      ("pcaMethods")
#BiocManager::install      ('MAST')
#BiocManager::install      ('org.Hs.eg.db')
#BiocManager::install      (version = "3.12")
#BiocManager::install      ('Seurat')
#BiocManager::install      ('readxl')
#BiocManager::install      ('modes')
#BiocManager::install      ('pheatmap')
#BiocManager::install      ('limma')
#BiocManager::install      ('clustree')
#BiocManager::install      ('clusterProfiler')
#BiocManager::install      ('EnhancedVolcano')
#install.packages          ('fields')
#install.packages          ('plotly')
#install.packages          ("VennDiagram")
 library                   (Seurat)
 library                   (Matrix)
 library                   (dplyr)
 library                   (ggplot2)
 library                   (cowplot)
 library                   (tibble)
 library                   (readxl)
 library                   (sctransform)
 library                   (fields)
 library                   (KernSmooth)
 library                   (ROCR)
 library                   (parallel)
 library                   (reshape2)
 library                   (pheatmap)
 library                   (DoubletFinder)
 library                   (limma)
 library                   (SeuratWrappers)
 library                   (SeuratDisk)
 library                   (plotly)
 library                   (clustree)
 # library                   (velocyto.R)
 library                   (MAST)
 library                   (EnhancedVolcano)
#library                   (clusterProfiler)
 library                   (AnnotationDbi)
 library                   (org.Hs.eg.db)
 library                   (VennDiagram)
```

## Read in cell cycle marker genes
 
```{r}
s.genes   <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
```

# Placental Villi Dataset Prepation, Quality Control, and Sample Integration

## Read in the Beristain single-cell Placenta data

```{r}
# Read in the data for placental villi sample one
getwd()
sample1 <- Read10X(data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P1")

# Read in the respective velocity data
#           ldat1            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/1_12weeks.loom")
# colnames (ldat1$spliced)   <- paste        (substring (colnames (ldat1$spliced),   15, 30), sep = "")
# colnames (ldat1$spliced)   <- paste                   (colnames (ldat1$spliced),   "-1",    sep = "")
# 
# colnames (ldat1$unspliced) <- paste        (substring (colnames (ldat1$unspliced), 15, 30), sep = "")
# colnames (ldat1$unspliced) <- paste                   (colnames (ldat1$unspliced), "-1",    sep = "")
# 
# colnames (ldat1$ambiguous) <- paste        (substring (colnames (ldat1$ambiguous), 15, 30), sep = "")
# colnames (ldat1$ambiguous) <- paste                   (colnames (ldat1$ambiguous), "-1",    sep = "")

# Convert to a seurat object
sample1 <- CreateSeuratObject (counts = sample1, min.cells = 1, min.features = 1, project = "B_P277")

# for (i in names (x = ldat1)) { sample1[[i]] <- CreateAssayObject(counts = ldat1[[i]]) }

# Add the sample metadata
sample1 <- AddMetaData        (object = sample1, metadata = "B_P277",   col.name = "ID"         ) # Sample identifier

sample1 <- AddMetaData        (object = sample1, metadata = "B_P277",   col.name = "rawID"      ) # Raw sample identifier (for samples with multiple files)

sample1 <- AddMetaData        (object = sample1, metadata = "1",        col.name = "nSex"       ) # Numeric biological sex of sample

sample1 <- AddMetaData        (object = sample1, metadata = "Male",     col.name = "Sex"        ) # Biological sex of sample

sample1 <- AddMetaData        (object = sample1, metadata = "6.4",      col.name = "rawGA"      ) # Raw gestational age of sample

sample1 <- AddMetaData        (object = sample1, metadata = "Early",    col.name = "binGA"      ) # Bin gestational age of sample

sample1 <- AddMetaData        (object = sample1, metadata = "1",        col.name = "nGA"        ) # Numeric gestational age of sample

sample1 <- AddMetaData        (object = sample1, metadata = "Early",    col.name = "GA"         ) # Gestational age of sample

sample1 <- AddMetaData        (object = sample1, metadata = "Placenta", col.name = "Tissue"     ) # Sample tissue type (either placenta, CD45+ Decidua, or CD45- Decidua)

sample1 <- AddMetaData        (object = sample1, metadata = "1",        col.name = "Batch"      ) # Sequencing batch

sample1 <- AddMetaData        (object = sample1, metadata = "1",        col.name = "Source"     ) # Data source (Beristain or Vento-Tormo)

sample1 <- AddMetaData        (object = sample1, metadata = "in vivo",  col.name = "Type"       ) # Data type (in vivo or in vitro)

sample1 <- AddMetaData        (object = sample1, metadata = "24",       col.name = "MA"         ) # Maternal Age (years)

sample1 <- AddMetaData        (object = sample1, metadata = "20.90",    col.name = "BMI"        ) # Patient Body Mass Index

sample1 <- AddMetaData        (object = sample1, metadata = "No",       col.name = "Smoke"      ) # Did the patient smoke during pregnancy?

sample1 <- AddMetaData        (object = sample1, metadata = "No",       col.name = "Diabetic"   ) # Is the patient living with diabetes (T1D or T2D)?

sample1 <- AddMetaData        (object = sample1, metadata = "None",     col.name = "Medications") # Did the patient consume anti-inflammatory and/or anti-hypertensive medications during pregnancy?

# Normalize the sample data and add cell cycle scores
sample1 <- NormalizeData      (sample1)

sample1 <- CellCycleScoring   (sample1, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample1)
dim   (sample1)
sample1@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample1 [["percent.mt"]]   <- PercentageFeatureSet (sample1, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample1), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample1, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample1, slot   = 'counts')              );

sample1 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample1, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample1, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample1, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample1@meta.data
```

```{r}
# Filter sample one
fltsample1 <- subset (sample1, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample1, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample1, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample1, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample one before and after filtering to observe the effect of filtering on this sample
   sample1
fltsample1
```

```{r}
normsample1 <- NormalizeData        ( fltsample1)
normsample1 <- ScaleData            (normsample1)
normsample1 <- FindVariableFeatures (normsample1, selection.method = "vst", nfeatures = 2000)
normsample1 <- RunPCA               (normsample1)
normsample1 <- RunUMAP              (normsample1, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample1 <- paramSweep_v3    (normsample1, PCs = 1:50, sct = FALSE)
sweep.stats_sample1    <- summarizeSweep   (sweep.res.list_sample1,  GT  = FALSE)
bcmvn_sample1          <- find.pK          (sweep.stats_sample1)
annotations1           <- normsample1@meta.data$ClusteringResults
homotypic.prop1        <- modelHomotypic   (annotations1)
nExp_poi1              <- round            (0.023 * length (normsample1@meta.data$ID))  ## Assuming 2.3% doublet formation rate
nExp_poi.adj1          <- round            (nExp_poi1 *    (1 - homotypic.prop1)     )
normsample1            <- doubletFinder_v3 (normsample1,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi1,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample1            <- doubletFinder_v3 (normsample1,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj1,
                                            reuse.pANN = "pANN_0.25_0.09_66",
                                            sct        = FALSE)
```

```{r}
normsample1@meta.data
DimPlot (normsample1, group.by = "DF.classifications_0.25_0.09_66")
table   (normsample1$DF.classifications_0.25_0.09_66)
```

```{r}
# Subset to remove identified cell doublets
fltsample1 <- subset (normsample1, subset = DF.classifications_0.25_0.09_66 == 'Singlet')
fltsample1
table (fltsample1$DF.classifications_0.25_0.09_66)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample1)
rm (sweep.res.list_sample1)
rm (sweep.stats_sample1)
rm (bcmvn_sample1)
rm (annotations1)
rm (homotypic.prop1)
rm (nExp_poi1)
rm (nExp_poi.adj1)
rm (normsample1)
rm (ldat1)
```

```{r}
# Read in the data for placental villi sample two
sample2 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P2")

# Read in the respective velocity data
#           ldat2            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/2_12weeks.loom")
# colnames (ldat2$spliced)   <- paste        (substring (colnames (ldat2$spliced),   15, 30), sep = "")
# colnames (ldat2$spliced)   <- paste                   (colnames (ldat2$spliced),   "-1",    sep = "")
# 
# colnames (ldat2$unspliced) <- paste        (substring (colnames (ldat2$unspliced), 15, 30), sep = "")
# colnames (ldat2$unspliced) <- paste                   (colnames (ldat2$unspliced), "-1",    sep = "")
# 
# colnames (ldat2$ambiguous) <- paste        (substring (colnames (ldat2$ambiguous), 15, 30), sep = "")
# colnames (ldat2$ambiguous) <- paste                   (colnames (ldat2$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample2 <- CreateSeuratObject (counts = sample2, min.cells = 1, min.features = 1, project = "B_P89")

# for (i in names (x = ldat2)) { sample2[[i]] <- CreateAssayObject(counts = ldat2[[i]]) }

# Add the sample metadata
sample2 <- AddMetaData        (object = sample2, metadata = "B_P89",    col.name = "ID"         )

sample2 <- AddMetaData        (object = sample2, metadata = "B_P89",    col.name = "rawID"      )

sample2 <- AddMetaData        (object = sample2, metadata = "0",        col.name = "nSex"       )

sample2 <- AddMetaData        (object = sample2, metadata = "Female",   col.name = "Sex"        )

sample2 <- AddMetaData        (object = sample2, metadata = "6.4",      col.name = "rawGA"      )

sample2 <- AddMetaData        (object = sample2, metadata = "Early",    col.name = "binGA"      )

sample2 <- AddMetaData        (object = sample2, metadata = "1",        col.name = "nGA"        )

sample2 <- AddMetaData        (object = sample2, metadata = "Early",    col.name = "GA"         )

sample2 <- AddMetaData        (object = sample2, metadata = "Placenta", col.name = "Tissue"     )

sample2 <- AddMetaData        (object = sample2, metadata = "1",        col.name = "Batch"      )

sample2 <- AddMetaData        (object = sample2, metadata = "1",        col.name = "Source"     )

sample2 <- AddMetaData        (object = sample2, metadata = "in vivo",  col.name = "Type"       )

sample2 <- AddMetaData        (object = sample2, metadata = "22",       col.name = "MA"         )

sample2 <- AddMetaData        (object = sample2, metadata = "24.21",    col.name = "BMI"        )

sample2 <- AddMetaData        (object = sample2, metadata = "No",       col.name = "Smoke"      )

sample2 <- AddMetaData        (object = sample2, metadata = "No",       col.name = "Diabetic"   )

sample2 <- AddMetaData        (object = sample2, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample2 <- NormalizeData      (sample2)

sample2 <- CellCycleScoring   (sample2, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample2)
dim   (sample2)
sample2@meta.data
```

```{r} 
# Add sample quality control metrics (% MT and % RP data)
sample2 [["percent.mt"]] <- PercentageFeatureSet (sample2, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample2), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample2, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample2, slot   = 'counts')              );

sample2 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample2, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample2, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample2, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample2@meta.data
```

```{r}
# Filter sample two
fltsample2 <- subset (sample2, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample2, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample2, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample2, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample two before and after filtering to observe the effect of filtering on this sample
   sample2
fltsample2
```

```{r}
normsample2 <- NormalizeData        ( fltsample2)
normsample2 <- ScaleData            (normsample2)
normsample2 <- FindVariableFeatures (normsample2, selection.method = "vst", nfeatures = 2000)
normsample2 <- RunPCA               (normsample2)
normsample2 <- RunUMAP              (normsample2, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample2 <- paramSweep_v3    (normsample2, PCs = 1:50, sct = FALSE)
sweep.stats_sample2    <- summarizeSweep   (sweep.res.list_sample2,  GT  = FALSE)
bcmvn_sample2          <- find.pK          (sweep.stats_sample2)
annotations2           <- normsample2@meta.data$ClusteringResults
homotypic.prop2        <- modelHomotypic   (annotations2)
nExp_poi2              <- round            (0.023 * length (normsample2@meta.data$ID))  ## Assuming 2.3% doublet formation rate
nExp_poi.adj2          <- round            (nExp_poi2 *    (1 - homotypic.prop2)     )
normsample2            <- doubletFinder_v3 (normsample2,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi2,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample2            <- doubletFinder_v3 (normsample2,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj2,
                                            reuse.pANN = "pANN_0.25_0.09_63",
                                            sct        = FALSE)
```

```{r}
normsample2@meta.data
DimPlot (normsample2, group.by = "DF.classifications_0.25_0.09_63")
table   (normsample2$DF.classifications_0.25_0.09_63)
```

```{r}
# Subset to remove identified cell doublets
fltsample2 <- subset (normsample2, subset = DF.classifications_0.25_0.09_63 == 'Singlet')
fltsample2
table (fltsample2$DF.classifications_0.25_0.09_63)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample2)
rm (sweep.res.list_sample2)
rm (sweep.stats_sample2)
rm (bcmvn_sample2)
rm (annotations2)
rm (homotypic.prop2)
rm (nExp_poi2)
rm (nExp_poi.adj2)
rm (normsample2)
rm (ldat2)
```

```{r}
# Read in the data for placental villi sample three
sample3 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P3")

# Read in the respective velocity data
#           ldat3            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/3_6weeks.loom")
# colnames (ldat3$spliced)   <- paste        (substring (colnames (ldat3$spliced),   15, 30), sep = "")
# colnames (ldat3$spliced)   <- paste                   (colnames (ldat3$spliced),   "-1",    sep = "")
# 
# colnames (ldat3$unspliced) <- paste        (substring (colnames (ldat3$unspliced), 15, 30), sep = "")
# colnames (ldat3$unspliced) <- paste                   (colnames (ldat3$unspliced), "-1",    sep = "")
# 
# colnames (ldat3$ambiguous) <- paste        (substring (colnames (ldat3$ambiguous), 15, 30), sep = "")
# colnames (ldat3$ambiguous) <- paste                   (colnames (ldat3$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample3 <- CreateSeuratObject (counts = sample3, min.cells = 1, min.features = 1, project = "B_P21")
# 
# for (i in names (x = ldat3)) { sample3[[i]] <- CreateAssayObject(counts = ldat3[[i]]) }

# Add the sample metadata
sample3 <- AddMetaData        (object = sample3, metadata = "B_P21",    col.name = "ID"         )

sample3 <- AddMetaData        (object = sample3, metadata = "B_P21",    col.name = "rawID"      )

sample3 <- AddMetaData        (object = sample3, metadata = "0",        col.name = "nSex"       )

sample3 <- AddMetaData        (object = sample3, metadata = "Female",   col.name = "Sex"        )

sample3 <- AddMetaData        (object = sample3, metadata = "11.2",     col.name = "rawGA"      )

sample3 <- AddMetaData        (object = sample3, metadata = "Late",     col.name = "binGA"      )

sample3 <- AddMetaData        (object = sample3, metadata = "2",        col.name = "nGA"        )

sample3 <- AddMetaData        (object = sample3, metadata = "Late",     col.name = "GA"         )

sample3 <- AddMetaData        (object = sample3, metadata = "Placenta", col.name = "Tissue"     )

sample3 <- AddMetaData        (object = sample3, metadata = "1",        col.name = "Batch"      )

sample3 <- AddMetaData        (object = sample3, metadata = "1",        col.name = "Source"     )

sample3 <- AddMetaData        (object = sample3, metadata = "in vivo",  col.name = "Type"       )

sample3 <- AddMetaData        (object = sample3, metadata = "30",       col.name = "MA"         )

sample3 <- AddMetaData        (object = sample3, metadata = "21.16",    col.name = "BMI"        )

sample3 <- AddMetaData        (object = sample3, metadata = "No",       col.name = "Smoke"      )

sample3 <- AddMetaData        (object = sample3, metadata = "No",       col.name = "Diabetic"   )

sample3 <- AddMetaData        (object = sample3, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample3 <- NormalizeData      (sample3)

sample3 <- CellCycleScoring   (sample3,s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample3)
dim   (sample3)
sample3@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample3 [["percent.mt"]] <- PercentageFeatureSet (sample3, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample3), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample3, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample3, slot   = 'counts')              );

sample3 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample3, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample3, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample3, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample3@meta.data
```

```{r}
# Filter sample three
fltsample3 <- subset (sample3, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample3, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample3, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample3, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample three before and after filtering to observe the effect of filtering on this sample
   sample3
fltsample3
```

```{r}
normsample3 <- NormalizeData        ( fltsample3)
normsample3 <- ScaleData            (normsample3)
normsample3 <- FindVariableFeatures (normsample3, selection.method = "vst", nfeatures = 2000)
normsample3 <- RunPCA               (normsample3)
normsample3 <- RunUMAP              (normsample3, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample3 <- paramSweep_v3    (normsample3, PCs = 1:50, sct = FALSE)
sweep.stats_sample3    <- summarizeSweep   (sweep.res.list_sample3,  GT  = FALSE)
bcmvn_sample3          <- find.pK          (sweep.stats_sample3)
annotations3           <- normsample3@meta.data$ClusteringResults
homotypic.prop3        <- modelHomotypic   (annotations3)
nExp_poi3             <- round             (0.015 * length (normsample3@meta.data$ID) )  ## Assuming 1.5% doublet formation rate
nExp_poi.adj3          <- round            (nExp_poi3 * (1 - homotypic.prop3)        )
normsample3            <- doubletFinder_v3 (normsample3,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi3,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample3            <- doubletFinder_v3 (normsample3,
                                            PCs        = 1:50, 
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj3,
                                            reuse.pANN = "pANN_0.25_0.09_25",
                                            sct        = FALSE)
```

```{r}
normsample3@meta.data
DimPlot (normsample3, group.by = "DF.classifications_0.25_0.09_25")
table   (normsample3$DF.classifications_0.25_0.09_25)
```

```{r}
# Subset to remove identified cell doublets
fltsample3 <- subset (normsample3, subset = DF.classifications_0.25_0.09_25 == 'Singlet')
fltsample3
table (fltsample3$DF.classifications_0.25_0.09_25)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample3)
rm (sweep.res.list_sample3)
rm (sweep.stats_sample3)
rm (bcmvn_sample3)
rm (annotations3)
rm (homotypic.prop3)
rm (nExp_poi3)
rm (nExp_poi.adj3)
rm (normsample3)
rm (ldat3)
```

```{r}
# Read in the data for placental villi sample four
sample4 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P4")

# Read in the respective velocity data
#           ldat4            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/4_6weeks.loom")
# colnames (ldat4$spliced)   <- paste        (substring (colnames (ldat4$spliced),   15, 30), sep = "")
# colnames (ldat4$spliced)   <- paste                   (colnames (ldat4$spliced),   "-1",    sep = "")
# 
# colnames (ldat4$unspliced) <- paste        (substring (colnames (ldat4$unspliced), 15, 30), sep = "")
# colnames (ldat4$unspliced) <- paste                   (colnames (ldat4$unspliced), "-1",    sep = "")
# 
# colnames (ldat4$ambiguous) <- paste        (substring (colnames (ldat4$ambiguous), 15, 30), sep = "")
# colnames (ldat4$ambiguous) <- paste                   (colnames (ldat4$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample4 <- CreateSeuratObject (counts = sample4, min.cells = 1, min.features = 1, project = "B_P130")

# for (i in names (x = ldat4)) { sample4[[i]] <- CreateAssayObject(counts = ldat4[[i]]) }

# Add the sample metadata
sample4 <- AddMetaData        (object = sample4, metadata = "B_P130",   col.name = "ID"         )

sample4 <- AddMetaData        (object = sample4, metadata = "B_P130",   col.name = "rawID"      )

sample4 <- AddMetaData        (object = sample4, metadata = "1",        col.name = "nSex"       )

sample4 <- AddMetaData        (object = sample4, metadata = "Male",     col.name = "Sex"        )

sample4 <- AddMetaData        (object = sample4, metadata = "11.4",     col.name = "rawGA"      )

sample4 <- AddMetaData        (object = sample4, metadata = "Late",     col.name = "binGA"      )

sample4 <- AddMetaData        (object = sample4, metadata = "2",        col.name = "nGA"        )

sample4 <- AddMetaData        (object = sample4, metadata = "Late",     col.name = "GA"         )

sample4 <- AddMetaData        (object = sample4, metadata = "Placenta", col.name = "Tissue"     )

sample4 <- AddMetaData        (object = sample4, metadata = "1",        col.name = "Batch"      )

sample4 <- AddMetaData        (object = sample4, metadata = "1",        col.name = "Source"     )

sample4 <- AddMetaData        (object = sample4, metadata = "in vivo",  col.name = "Type"       )

sample4 <- AddMetaData        (object = sample4, metadata = "32",       col.name = "MA"         )

sample4 <- AddMetaData        (object = sample4, metadata = "22.86",    col.name = "BMI"        )

sample4 <- AddMetaData        (object = sample4, metadata = "No",       col.name = "Smoke"      )

sample4 <- AddMetaData        (object = sample4, metadata = "No",       col.name = "Diabetic"   )

sample4 <- AddMetaData        (object = sample4, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample4 <- NormalizeData      (sample4)

sample4 <- CellCycleScoring   (sample4, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample4)
dim   (sample4)
sample4@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample4 [["percent.mt"]] <- PercentageFeatureSet (sample4, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample4), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample4, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample4, slot   = 'counts')              );

sample4 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample4, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample4, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample4, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample4@meta.data
```

```{r}
# Filter sample four
fltsample4 <- subset (sample4, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot and decide threshold for % mtRNA
VlnPlot (fltsample4, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample4, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample4, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample four before and after filtering to observe the effect of filtering on this sample
   sample4
fltsample4
```

```{r}
normsample4 <- NormalizeData        ( fltsample4)
normsample4 <- ScaleData            (normsample4)
normsample4 <- FindVariableFeatures (normsample4, selection.method = "vst", nfeatures = 2000)
normsample4 <- RunPCA               (normsample4)
normsample4 <- RunUMAP              (normsample4, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample4 <- paramSweep_v3    (normsample4, PCs = 1:50, sct = FALSE)
sweep.stats_sample4    <- summarizeSweep   (sweep.res.list_sample4,  GT  = FALSE)
bcmvn_sample4          <- find.pK          (sweep.stats_sample4)
annotations4           <- normsample4@meta.data$ClusteringResults
homotypic.prop4        <- modelHomotypic   (annotations4)
nExp_poi4              <- round            (0.023 * length (normsample4@meta.data$ID))  ## Assuming 2.3% doublet formation rate
nExp_poi.adj4          <- round            (nExp_poi4 *    (1 - homotypic.prop4)     )
normsample4            <- doubletFinder_v3 (normsample4,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi4,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample4            <- doubletFinder_v3 (normsample4,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj4,
                                            reuse.pANN = "pANN_0.25_0.09_67",
                                            sct        = FALSE)
```

```{r}
normsample4@meta.data
DimPlot (normsample4, group.by = "DF.classifications_0.25_0.09_67")
table   (normsample4$DF.classifications_0.25_0.09_67)
```

```{r}
# Subset to remove identified cell doublets
fltsample4 <- subset (normsample4, subset = DF.classifications_0.25_0.09_67 == 'Singlet')
fltsample4
table (fltsample4$DF.classifications_0.25_0.09_67)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample4)
rm (sweep.res.list_sample4)
rm (sweep.stats_sample4)
rm (bcmvn_sample4)
rm (annotations4)
rm (homotypic.prop4)
rm (nExp_poi4)
rm (nExp_poi.adj4)
rm (normsample4)
rm (ldat4)
```

```{r}
# Read in the data for placental villi sample five
sample5 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P5")

# Read in the respective velocity data
#           ldat5            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/sc244_6w.loom")
# colnames (ldat5$spliced)   <- paste        (substring (colnames (ldat5$spliced),   18, 33), sep = "")
# colnames (ldat5$spliced)   <- paste                   (colnames (ldat5$spliced),   "-1",    sep = "")
# 
# colnames (ldat5$unspliced) <- paste        (substring (colnames (ldat5$unspliced), 18, 33), sep = "")
# colnames (ldat5$unspliced) <- paste                   (colnames (ldat5$unspliced), "-1",    sep = "")
# 
# colnames (ldat5$ambiguous) <- paste        (substring (colnames (ldat5$ambiguous), 18, 33), sep = "")
# colnames (ldat5$ambiguous) <- paste                   (colnames (ldat5$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample5 <- CreateSeuratObject (counts = sample5, min.cells = 1, min.features = 1, project = "B_P244")
# 
# for (i in names (x = ldat5)) { sample5[[i]] <- CreateAssayObject(counts = ldat5[[i]]) }

# Add the sample metadata
sample5 <- AddMetaData        (object = sample5, metadata = "B_P244",   col.name = "ID"         )

sample5 <- AddMetaData        (object = sample5, metadata = "B_P244",   col.name = "rawID"      )

sample5 <- AddMetaData        (object = sample5, metadata = "0",        col.name = "nSex"       )

sample5 <- AddMetaData        (object = sample5, metadata = "Female",   col.name = "Sex"        )

sample5 <- AddMetaData        (object = sample5, metadata = "6.2",      col.name = "rawGA"      )

sample5 <- AddMetaData        (object = sample5, metadata = "Early",    col.name = "binGA"      )

sample5 <- AddMetaData        (object = sample5, metadata = "1",        col.name = "nGA"        )

sample5 <- AddMetaData        (object = sample5, metadata = "Early",    col.name = "GA"         )

sample5 <- AddMetaData        (object = sample5, metadata = "Placenta", col.name = "Tissue"     )

sample5 <- AddMetaData        (object = sample5, metadata = "2",        col.name = "Batch"      )

sample5 <- AddMetaData        (object = sample5, metadata = "1",        col.name = "Source"     )

sample5 <- AddMetaData        (object = sample5, metadata = "in vivo",  col.name = "Type"       )

sample5 <- AddMetaData        (object = sample5, metadata = "23",       col.name = "MA"         )

sample5 <- AddMetaData        (object = sample5, metadata = "23.53",    col.name = "BMI"        )

sample5 <- AddMetaData        (object = sample5, metadata = "No",       col.name = "Smoke"      )

sample5 <- AddMetaData        (object = sample5, metadata = "No",       col.name = "Diabetic"   )

sample5 <- AddMetaData        (object = sample5, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample5 <- NormalizeData      (sample5)

sample5 <- CellCycleScoring   (sample5, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample5)
dim   (sample5)
sample5@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample5 [["percent.mt"]] <- PercentageFeatureSet (sample5, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample5), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample5, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample5, slot   = 'counts')              );

sample5 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample5, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample5, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample5, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample5@meta.data
```

```{r}
# Filter sample five
fltsample5 <- subset (sample5, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample5, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample5, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample5, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample five before and after filtering to observe the effect of filtering on this sample
   sample5
fltsample5
```

```{r}
normsample5 <- NormalizeData        ( fltsample5)
normsample5 <- ScaleData            (normsample5)
normsample5 <- FindVariableFeatures (normsample5, selection.method = "vst", nfeatures = 2000)
normsample5 <- RunPCA               (normsample5)
normsample5 <- RunUMAP              (normsample5, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample5 <- paramSweep_v3    (normsample5, PCs = 1:50, sct = FALSE)
sweep.stats_sample5    <- summarizeSweep   (sweep.res.list_sample5,  GT  = FALSE)
bcmvn_sample5          <- find.pK          (sweep.stats_sample5)
annotations5           <- normsample5@meta.data$ClusteringResults
homotypic.prop5        <- modelHomotypic   (annotations5)
nExp_poi5              <- round            (0.004 * length (normsample5@meta.data$ID))  ## Assuming 0.4% doublet formation rate
nExp_poi.adj5          <- round            (nExp_poi5 *    (1 - homotypic.prop5)     )
normsample5            <- doubletFinder_v3 (normsample5,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi5,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample5            <- doubletFinder_v3 (normsample5,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj5,
                                            reuse.pANN = "pANN_0.25_0.09_2",
                                            sct        = FALSE)
```

```{r}
normsample5@meta.data
DimPlot (normsample5, group.by = "DF.classifications_0.25_0.09_2")
table   (normsample5$DF.classifications_0.25_0.09_2)
```

```{r}
# Subset to remove identified cell doublets
fltsample5 <- subset (normsample5, subset = DF.classifications_0.25_0.09_2 == 'Singlet')
fltsample5
table (fltsample5$DF.classifications_0.25_0.09_2)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample5)
rm (sweep.res.list_sample5)
rm (sweep.stats_sample5)
rm (bcmvn_sample5)
rm (annotations5)
rm (homotypic.prop5)
rm (nExp_poi5)
rm (nExp_poi.adj5)
rm (normsample5)
rm (ldat5)
```

```{r}
# Read in the data for placental villi sample six
sample6 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P6")

# Read in the respective velocity data
#           ldat6            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/sc295_6w.loom")
# colnames (ldat6$spliced)   <- paste        (substring (colnames (ldat6$spliced),   18, 33), sep = "")
# colnames (ldat6$spliced)   <- paste                   (colnames (ldat6$spliced),   "-1",    sep = "")
# 
# colnames (ldat6$unspliced) <- paste        (substring (colnames (ldat6$unspliced), 18, 33), sep = "")
# colnames (ldat6$unspliced) <- paste                   (colnames (ldat6$unspliced), "-1",    sep = "")
# 
# colnames (ldat6$ambiguous) <- paste        (substring (colnames (ldat6$ambiguous), 18, 33), sep = "")
# colnames (ldat6$ambiguous) <- paste                   (colnames (ldat6$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample6 <- CreateSeuratObject (counts = sample6, min.cells = 1, min.features = 1, project = "B_P295")

# for (i in names (x = ldat6)) { sample6[[i]] <- CreateAssayObject(counts = ldat6[[i]]) }

# Add the sample metadata
sample6 <- AddMetaData        (object = sample6, metadata = "B_P295",   col.name = "ID"         )

sample6 <- AddMetaData        (object = sample6, metadata = "B_P295",   col.name = "rawID"      )

sample6 <- AddMetaData        (object = sample6, metadata = "1",        col.name = "nSex"       )

sample6 <- AddMetaData        (object = sample6, metadata = "Male",     col.name = "Sex"        )

sample6 <- AddMetaData        (object = sample6, metadata = "6.5",      col.name = "rawGA"      )

sample6 <- AddMetaData        (object = sample6, metadata = "Early",    col.name = "binGA"      )

sample6 <- AddMetaData        (object = sample6, metadata = "1",        col.name = "nGA"        )

sample6 <- AddMetaData        (object = sample6, metadata = "Early",    col.name = "GA"         )

sample6 <- AddMetaData        (object = sample6, metadata = "Placenta", col.name = "Tissue"     )

sample6 <- AddMetaData        (object = sample6, metadata = "2",        col.name = "Batch"      )

sample6 <- AddMetaData        (object = sample6, metadata = "1",        col.name = "Source"     )

sample6 <- AddMetaData        (object = sample6, metadata = "in vivo",  col.name = "Type"       )

sample6 <- AddMetaData        (object = sample6, metadata = "30",       col.name = "MA"         )

sample6 <- AddMetaData        (object = sample6, metadata = "23.51",    col.name = "BMI"        )

sample6 <- AddMetaData        (object = sample6, metadata = "No",       col.name = "Smoke"      )

sample6 <- AddMetaData        (object = sample6, metadata = "No",       col.name = "Diabetic"   )

sample6 <- AddMetaData        (object = sample6, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample6 <- NormalizeData      (sample6)

sample6 <- CellCycleScoring   (sample6, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample6)
dim   (sample6)
sample6@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample6 [["percent.mt"]] <- PercentageFeatureSet (sample6, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample6), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample6, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample6, slot   = 'counts')              );

sample6 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample6, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample6, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample6, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample6@meta.data
```

```{r}
# Filter sample six
fltsample6 <- subset (sample6, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample6, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample6, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample6, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample six before and after filtering to observe the effect of filtering on this sample
   sample6
fltsample6
```

```{r}
normsample6 <- NormalizeData        ( fltsample6)
normsample6 <- ScaleData            (normsample6)
normsample6 <- FindVariableFeatures (normsample6, selection.method = "vst", nfeatures = 2000)
normsample6 <- RunPCA               (normsample6)
normsample6 <- RunUMAP              (normsample6, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample6 <- paramSweep_v3    (normsample6, PCs = 1:50, sct = FALSE)
sweep.stats_sample6    <- summarizeSweep   (sweep.res.list_sample6,  GT  = FALSE)
bcmvn_sample6          <- find.pK          (sweep.stats_sample6)
annotations6           <- normsample6@meta.data$ClusteringResults
homotypic.prop6        <- modelHomotypic   (annotations6)
nExp_poi6              <- round            (0.008 * length (normsample6@meta.data$ID))  ## Assuming 0.8% doublet formation rate
nExp_poi.adj6          <- round            (nExp_poi6 *    (1 - homotypic.prop6)     )
normsample6            <- doubletFinder_v3 (normsample6,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi6,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample6            <- doubletFinder_v3 (normsample6,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj6,
                                            reuse.pANN = "pANN_0.25_0.09_6",
                                            sct        = FALSE)
```

```{r}
normsample6@meta.data
DimPlot (normsample6, group.by = "DF.classifications_0.25_0.09_6")
table   (normsample6$DF.classifications_0.25_0.09_6)
```

```{r}
# Subset to remove identified cell doublets
fltsample6 <- subset (normsample6, subset = DF.classifications_0.25_0.09_6 == 'Singlet')
fltsample6
table (fltsample6$DF.classifications_0.25_0.09_6)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample6)
rm (sweep.res.list_sample6)
rm (sweep.stats_sample6)
rm (bcmvn_sample6)
rm (annotations6)
rm (homotypic.prop6)
rm (nExp_poi6)
rm (nExp_poi.adj6)
rm (normsample6)
rm (ldat6)
```

```{r}
# Read in the data for placental villi sample seven
sample7 <- Read10X            (data.dir = "/Users/young-junekim/Desktop/CAREER/Graduate Studies in Physiology/Year 1/Labwork/Meta-analysis Work/Datasets/Individuals/Shannon et al. (2022)/P7")

# Read in the respective velocity data
#           ldat7            <- ReadVelocity ("/Volumes/OneTouch/R/5_scVelo/sc899_12w.loom")
# colnames (ldat7$spliced)   <- paste        (substring (colnames (ldat7$spliced),   18, 33), sep = "")
# colnames (ldat7$spliced)   <- paste                   (colnames (ldat7$spliced),   "-1",    sep = "")
# 
# colnames (ldat7$unspliced) <- paste        (substring (colnames (ldat7$unspliced), 18, 33), sep = "")
# colnames (ldat7$unspliced) <- paste                   (colnames (ldat7$unspliced), "-1",    sep = "")
# 
# colnames (ldat7$ambiguous) <- paste        (substring (colnames (ldat7$ambiguous), 18, 33), sep = "")
# colnames (ldat7$ambiguous) <- paste                   (colnames (ldat7$ambiguous), "-1",    sep = "")

# Convert to a seurat object for analysis
sample7 <- CreateSeuratObject (counts = sample7, min.cells = 1, min.features = 1, project = "B_P899")
# 
# for (i in names (x = ldat7)) { sample7[[i]] <- CreateAssayObject(counts = ldat7[[i]]) }

# Add the sample metadata
sample7 <- AddMetaData        (object = sample7, metadata = "B_P899",   col.name = "ID"         )

sample7 <- AddMetaData        (object = sample7, metadata = "B_P899",   col.name = "rawID"      )

sample7 <- AddMetaData        (object = sample7, metadata = "0",        col.name = "nSex"       )

sample7 <- AddMetaData        (object = sample7, metadata = "Female",   col.name = "Sex"        )

sample7 <- AddMetaData        (object = sample7, metadata = "11.4",     col.name = "rawGA"      )

sample7 <- AddMetaData        (object = sample7, metadata = "Late",     col.name = "binGA"      )

sample7 <- AddMetaData        (object = sample7, metadata = "2",        col.name = "nGA"        )

sample7 <- AddMetaData        (object = sample7, metadata = "Late",     col.name = "GA"         )

sample7 <- AddMetaData        (object = sample7, metadata = "Placenta", col.name = "Tissue"     )

sample7 <- AddMetaData        (object = sample7, metadata = "2",        col.name = "Batch"      )

sample7 <- AddMetaData        (object = sample7, metadata = "1",        col.name = "Source"     )

sample7 <- AddMetaData        (object = sample7, metadata = "in vivo",  col.name = "Type"       )

sample7 <- AddMetaData        (object = sample7, metadata = "20",       col.name = "MA"         )

sample7 <- AddMetaData        (object = sample7, metadata = "22.08",    col.name = "BMI"        )

sample7 <- AddMetaData        (object = sample7, metadata = "No",       col.name = "Smoke"      )

sample7 <- AddMetaData        (object = sample7, metadata = "No",       col.name = "Diabetic"   )

sample7 <- AddMetaData        (object = sample7, metadata = "None",     col.name = "Medications")

# Normalize the sample data and add cell cycle scores
sample7 <- NormalizeData      (sample7)

sample7 <- CellCycleScoring   (sample7, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# Observe the sample data for quality control
class (sample7)
dim   (sample7)
sample7@meta.data
```

```{r}
# Add sample quality control metrics (% MT and % RP data)
sample7 [["percent.mt"]] <- PercentageFeatureSet (sample7, pattern = "^MT-")

ribo.genes                 <- grep (pattern = "^RP[SL][[:digit:]]", x = rownames (sample7), value = TRUE);

percent.ribo               <- Matrix::colSums (x = GetAssayData (sample7, slot   = 'counts')[ribo.genes, ]) / 
                              Matrix::colSums (x = GetAssayData (sample7, slot   = 'counts')              );

sample7 [['percent.ribo']] <- percent.ribo

# Visualize QC metrics as a violin plot, cells with >20% mtRNA are taken as poor quality and removed 
VlnPlot (sample7, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (sample7, features = c("percent.mt"),   group.by = "ID")
VlnPlot (sample7, features = c("percent.ribo"), group.by = "ID")

# View sample metadata
sample7@meta.data
```

```{r}
# Filter sample seven
fltsample7 <- subset (sample7, subset = nFeature_RNA > 1000 & percent.mt < 20)

# Visualize QC metrics as a violin plot again (confirm filtering success)
VlnPlot (fltsample7, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)
VlnPlot (fltsample7, features = c("percent.mt"),   group.by = "ID")
VlnPlot (fltsample7, features = c("percent.ribo"), group.by = "ID")

# View the structure of sample seven before and after filtering to observe the effect of filtering on this sample
   sample7
fltsample7
```

```{r}
normsample7 <- NormalizeData        ( fltsample7)
normsample7 <- ScaleData            (normsample7)
normsample7 <- FindVariableFeatures (normsample7, selection.method = "vst", nfeatures = 2000)
normsample7 <- RunPCA               (normsample7)
normsample7 <- RunUMAP              (normsample7, dims = 1:50)
```

```{r}
# Remove doublets from this sample
sweep.res.list_sample7 <- paramSweep_v3    (normsample7, PCs = 1:50, sct = FALSE)
sweep.stats_sample7    <- summarizeSweep   (sweep.res.list_sample7,  GT  = FALSE)
bcmvn_sample7          <- find.pK          (sweep.stats_sample7)
annotations7           <- normsample7@meta.data$ClusteringResults
homotypic.prop7        <- modelHomotypic   (annotations7)
nExp_poi7              <- round            (0.016 * length (normsample7@meta.data$ID))  ## Assuming 1.6% doublet formation rate
nExp_poi.adj7          <- round            (nExp_poi7 *    (1 - homotypic.prop7)     )
normsample7            <- doubletFinder_v3 (normsample7,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi7,
                                            reuse.pANN = FALSE,
                                            sct        = FALSE)
normsample7            <- doubletFinder_v3 (normsample7,
                                            PCs        = 1:50,
                                            pN         = 0.25,
                                            pK         = 0.09,
                                            nExp       = nExp_poi.adj7,
                                            reuse.pANN = "pANN_0.25_0.09_26",
                                            sct        = FALSE)
```

```{r}
normsample7@meta.data
DimPlot (normsample7, group.by = "DF.classifications_0.25_0.09_26")
table   (normsample7$DF.classifications_0.25_0.09_26)
```

```{r}
# Subset to remove identified cell doublets
fltsample7 <- subset (normsample7, subset = DF.classifications_0.25_0.09_26 == 'Singlet')
fltsample7
table (fltsample7$DF.classifications_0.25_0.09_26)
```

```{r}
# Remove files no longer needed before proceeding (conserves computational space)
rm (sample7)
rm (sweep.res.list_sample7)
rm (sweep.stats_sample7)
rm (bcmvn_sample7)
rm (annotations7)
rm (homotypic.prop7)
rm (nExp_poi7)
rm (nExp_poi.adj7)
rm (normsample7)
rm (ldat7)
```

## Merge the Beristain and Vento-Tormo Placenta datasets

```{r # Run in Terminal}
# Merge the twenty-one placental villi seurat objects into one large data file
vel_interface <- merge (fltsample1, y = c(fltsample2 ,
                                          fltsample3 ,
                                          fltsample4 ,
                                          fltsample5 ,
                                          fltsample6 ,
                                          fltsample7),
                        add.cell.ids  = c("1", "2", "3", "4", "5", "6", "7"), # Beristain Placenta Datasets
                        merge.data    = TRUE,
                        project       = "10X_Integrated_Interface")

# View the merged placental villi seurat object for visual inspection before downstream analyses
vel_interface
```

```{r}
# View the merged sample data stratified by sample ID
table (vel_interface$ID)
```

```{r}
# View the dimensionality of the merged dataset
dim (vel_interface)
```

```{r}
# Visualize the distribution of cell cycle markers across the merged dataset
RidgePlot (vel_interface, features = c("PCNA", "TOP2A", "MCM6", "CDK1", "MKI67", "CCNA2"), ncol = 3)
```

```{r}
# Add metadata for the cell cycle difference scores, this will allow us to normalize for any variation present within a cell cycle stage but not for variation that exists between cell cycle stages (allows us to still use cell proliferation as an indicator for cell clustering without confounding variance of cell cycle gene expression changes that exist within distinct cell cycle states)
vel_interface$CC.Difference <- vel_interface$S.Score - vel_interface$G2M.Score
```

## Visualize metadata and mitochondrial DNA information

```{r}
# View the merged metadata
vel_interface@meta.data
```

```{r}
# Visualize final dataset QC metrics using violoin plots
VlnPlot(object   = vel_interface,
        features = c("nCount_RNA", "nFeature_RNA", "percent.mt", "percent.ribo"),
        group.by =   "ID",
        ncol     = 2)
```

```{r}
# FeatureScatter to visualize feature-feature relationships
P1 <- FeatureScatter (object = vel_interface, feature1 = "nCount_RNA", feature2 = "percent.mt"  ) + NoLegend()
P2 <- FeatureScatter (object = vel_interface, feature1 = "nCount_RNA", feature2 = "percent.ribo") + NoLegend()
P3 <- FeatureScatter (object = vel_interface, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") + NoLegend()

CombinePlots (plots = list (P1, P2, P3) )
rm           (P1)
rm           (P2)
rm           (P3)
```

## Integrate All Datasets into One Seurat File

```{r}
# Integrate the data based on the "Batch" variable
vel_interface.list <- SplitObject (vel_interface, split.by = "Batch")

vel_interface.list <- vel_interface.list [ c("1", "2") ]

```

```{r}
# Normalize the data and regress out variation resulting from the % MT and CC Difference variables
for (i in 1:length(vel_interface.list)) {

vel_interface.list[[i]] <- SCTransform (vel_interface.list[[i]],
                                        verbose         = FALSE,
                                        vars.to.regress = c("percent.mt",
                                                            "CC.Difference" ))
                                         }
```

```{r}
# Increase global max size to sustain heavy computational loads
options (future.globals.maxSize = 10000000000000000000000000000000)
```

```{r}
# Integrate the data
vel_interface.features <- SelectIntegrationFeatures (object.list     = vel_interface.list, nfeatures = 3000)

vel_interface.list     <- PrepSCTIntegration        (object.list     = vel_interface.list,
                                                     anchor.features = vel_interface.features,
                                                     verbose         = FALSE)
``` 

```{r}
# Integrate the data (continued)
vel_interface.anchors <- FindIntegrationAnchors (object.list          = vel_interface.list,
                                                 normalization.method = "SCT",
                                                 anchor.features      = vel_interface.features,
                                                 verbose              = FALSE)
rm (vel_interface.list)
rm (vel_interface.features)
```

```{r}
# Integrate the data (continued)
Interface <- IntegrateData (anchorset            = vel_interface.anchors,
                            normalization.method = "SCT",
                            verbose              = FALSE)
rm (vel_interface.anchors)
```

```{r}
# Run principal component analysis on integrated data
Interface <- RunPCA (object         = Interface,
                     do.print       = TRUE,
                     pcs.print      = 1:5,
                     features.print = 5)
```

```{r}
# Visualize explained variation as a function of the number of principal components
ElbowPlot (Interface, ndims = 50)
```

```{r}
# Principal component analysis visualiztion, colour-coded by selected metadata variables
PCAPlot (Interface, group.by = "ID")
PCAPlot (Interface, group.by = "GA")
PCAPlot (Interface, group.by = "Sex")
PCAPlot (Interface, group.by = "Phase")
```

```{r}
# Additional data visualization
VizDimLoadings (Interface, dims = 1:4)
```

```{r}
# Cluster the maternal fetal interface dataset
resolution.range     <- seq (from = 0, to = 1, by = 0.1)

Interface <- FindNeighbors (Interface)
Interface <- FindClusters  (Interface, resolution = resolution.range)
```

```{r}
# Select maternal fetal interface organoid cluster resolution
clustree (Interface, prefix = "integrated_snn_res.")
```

```{r}
Interface <- FindClusters (Interface, resolution = 1.0)
# Change the resolution to get 7 or 8 cell types
# Reduce maternal fetal interface data dimensionality
Interface <- RunUMAP      (Interface, dims       = 1:50)
Interface <- RunTSNE      (Interface, dims       = 1:50)
```

```{r}
# Visualize the maternal fetal interface data by UMAP
UMAPPlot    (Interface, label = TRUE,  pt.size = 2)                            + NoLegend ()
UMAPPlot    (Interface, label = TRUE,  pt.size = 2, split.by = 'GA')           + NoLegend ()
UMAPPlot    (Interface, label = TRUE,  pt.size = 2, split.by = 'Sex')          + NoLegend ()
UMAPPlot    (Interface, label = TRUE,  pt.size = 2, split.by = 'Tissue')       + NoLegend ()
UMAPPlot    (Interface, label = FALSE, pt.size = 2, split.by = 'ID', ncol = 3) + NoLegend ()

FeaturePlot (Interface, features = c("nFeature_RNA"), min.cutoff = 0, pt.size = 2, order = TRUE, cols = c("grey", "red"))
FeaturePlot (Interface, features = c("nCount_RNA"),   min.cutoff = 0, pt.size = 2, order = TRUE, cols = c("grey", "red"))
FeaturePlot (Interface, features = c("percent.mt"),   min.cutoff = 0, pt.size = 2, order = TRUE, cols = c("grey", "red"))
FeaturePlot (Interface, features = c("percent.ribo"), min.cutoff = 0, pt.size = 2, order = TRUE, cols = c("grey", "red"))

save (Interface, file = "interface.Rdata")
```

Let's extract the cluster-specific markers at this point.

```{r}
getwd()
Merged <- load("/home/youkim/Meta-analysis/Shannon et al. (2021)/interface.RData")
Markers <- FindAllMarkers(Interface, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(Markers)

cat(tapply(Markers$gene,list(Markers$cluster),function(i) paste0(i[1:10],collapse = ',')),sep='\n')

head(Idents(Interface))
```

```{r}
# Label all maternal fetal interface clusters by cell type
new.cluster.ids       <- c('Syncytiotrophoblast',
                           'Endothelial cells',
                           'Proliferating trophoblasts',
                           'Pericytes',
                           'Hofbauer cells',
                           'Mesenchymal stem cells',
                           'Extravillous trophoblasts',
                           'Cytotrophoblasts','Fibroblasts',
                           'Smooth muscle cells',
                           'Antigen-presenting cells',
                           'Villous cytotrophoblasts',
                           'Vascular endothelial cells',
                           'Decidual cells',
                           'Proliferating cells',
                           'Macrophages',
                           'Immune cells',
                           'Myofibroblasts',
                           'Stromal cells')

names (new.cluster.ids)      <- levels       (Interface)
Interface$seurat_clusters <- Idents       (Interface)
Interface                 <- RenameIdents (Interface, new.cluster.ids)

saveRDS(Interface, "shSeurat_GPT4.rds")
shdataset_umap <- UMAPPlot(Interface) + NoLegend()
LabelClusters(plot = shdataset_umap, id = "ident")
```

```{r}
# Save the maternal fetal interface dataset

```

