All genes, ‘EM’ samples only +EM, Minus E3 (40 pTE samples per stage)

library(destiny)

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

library(RColorBrewer)
library(scater)

## Loading required package: SingleCellExperiment

## Loading required package: SummarizedExperiment

## Loading required package: MatrixGenerics

## Loading required package: matrixStats

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: BiocGenerics

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:destiny':
## 
##     distance

## The following object is masked from 'package:grDevices':
## 
##     windows

## Loading required package: GenomeInfoDb

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

## Loading required package: scuttle

## Loading required package: ggplot2

library(SingleCellExperiment)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:IRanges':
## 
##     slice

## The following object is masked from 'package:S4Vectors':
## 
##     rename

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

setwd("C:\\Users\\Jenni\\OneDrive\\PhD\\HumanEmbryoImplantation\\Run1\\EBremoved\\trajectory")
#read in a counts table and metadata file to make your own seurat object
counts<-read.csv("DEA_TEMcounts.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_mTE_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning: 'as(<dsCMatrix>, "dsTMatrix")' is deprecated.
## Use 'as(., "TsparseMatrix")' instead.
## See help("Deprecated") and help("Matrix-deprecated").

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, 
  color = ~Timepoint, colors = c("skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers(size=10) %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="mTE Diffusion Plot")

p

counts<-read.csv("DEA_TEPcounts.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_pTE_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="pTE Diffusion Plot")

p

counts<-read.csv("DEA_EPIcounts_LOG2FC1.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_EPI_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 1688 genes. Consider passing e.g.
## n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="EPI Diffusion Plot")

p

counts<-read.csv("DEA_PEcounts_LOG2FC1.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_PE_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 545 genes. Consider passing e.g.
## n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="PE Diffusion Plot")

p

counts<-read.csv("Pet-counts-ALL.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_ALL.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 17511 genes. Consider passing
## e.g. n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("green", "darkgreen", "tan", "skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="ALL Diffusion Plot, Including E3")

p

counts<-read.csv("Pet-counts-ALL-minusE3.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_ALL_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 17511 genes. Consider passing
## e.g. n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("green", "darkgreen", "skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="ALL Diffusion Plot, Minus E3")

p

counts<-read.csv("Pet-counts-ALLgenes-minusE3-onlyTEsamples.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_TE_AB_minusE3.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 17511 genes. Consider passing
## e.g. n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("green", "darkgreen", "skyblue", "deepskyblue", "blue", "navy")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="ALL genes: Diffusion Plot, Minus E3, TE samples only +AB")

p

counts<-read.csv("Pet-counts-ALLgenes-minusE3-onlyEMsamples-40pTEperstage.csv", header = T)
counts$X<-as.character(counts$X)
rownames(counts)<-make.unique(counts$X)
counts<-counts[,-1]

labels=read.table("Pet_meta_ICMpTE_EM_minusE3-40pTEperstage.txt",header=TRUE)

# create a SingleCellExperiment object
sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(counts),
    logcounts = log2(as.matrix(counts) + 1)
  ), 
  colData = labels
)

sce_log <- logcounts(sce)
dm <- DiffusionMap(t(sce_log))

## Warning in DiffusionMap(t(sce_log)): You have 17511 genes. Consider passing
## e.g. n_pcs = 50 to speed up computation.

tmp <- data.frame(DC1 = eigenvectors(dm)[,1],
                  DC2 = eigenvectors(dm)[,2],
                  DC3 = eigenvectors(dm)[,3],
                  Timepoint = sce$stage)


p <- plot_ly(
  tmp, x = ~DC1, y = ~DC2, z = ~DC3, size=5,
  color = ~Timepoint, colors = c("skyblue", "deepskyblue", "blue", "navy", "red2", "orange")
  ) %>%
  add_markers() %>%
  layout(
    scene = list(xaxis = list(title = 'DC1'),
        yaxis = list(title = 'DC2'),
        zaxis = list(title = 'DC3'))
        ) %>%
  layout(title="ALL genes: Diffusion Plot, Minus E3, EPI/PE/pTE samples only +EM (40 pTE samples per stage)")

p