---
title: "echibnoScElifeCodeSample"
author: "Annie"
date: '2023-02-10'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#session set up
##load packages
```{r}
library(Seurat)
library(dplyr)
library(patchwork)

```

-------------------------------------------------
This file contains an outline of the pipelne used to create the P. miniata nuclear atlas, the S. purpuratus whole cell atlas, and the integrated atlas. 

#Pminiata atlas creation 
  
  
**NOTE: for 1:1 orthologs in P. miniata data, the gene symbol was replaced with the symbol used for that gene in S. pupuratus**
  
## open raw data files
read in raw cellranger output
```{r}
pm16<-readRDS(file="./pm16hpf.rds")
pm26<-readRDS(file="./pm26hpf.rds")
pm40<-readRDS(file="./pm40hpf.rds")
```
  

##quality control 
1. quantify the percentage of transcripts in each file that are rRNAs
  -rRNAs previousy identified using gff file
```{r}
rRNA<-c("LOC590670",'LOC119727991','LOC581711','LOC119719604','LOC119719678',"LOC594000","LOC752098","LOC105445410","LOC119729016",'LOC119730276',"LOC582829","LOC119733178","LOC105447583","LOC582250","LOC592405","LOC589889","LOC119738685","LOC119741035",'LOC587799','LOC119719424') 

pm16[["percent.rRNA"]] <- PercentageFeatureSet(pm16, features = intersect(rownames(pm16),rRNA))
pm26[["percent.rRNA"]] <- PercentageFeatureSet(pm26, features = intersect(rownames(pm26),rRNA))
pm40[["percent.rRNA"]] <- PercentageFeatureSet(pm40, features = intersect(rownames(pm40),rRNA))
```
  
visualize distribution of raw data
```{r}
VlnPlot(pm16,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
VlnPlot(pm26,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
VlnPlot(pm40,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
```
  
2. subset data to apply determined cutoffs for nFeature_RNA, nCount_RNA, and percenrt.rRNA
```{r}
pm16 <-subset(pm16, nFeature_RNA < 2000 &  nFeature_RNA > 200 & nCount_RNA < 2600 & nCount_RNA > 200)
pm26 <-subset(pm26,nFeature_RNA > 200 & nFeature_RNA < 1000 & nCount_RNA < 1500 & nCount_RNA > 200) 
pm40 <-subset(pm40,nFeature_RNA > 200 & nFeature_RNA < 2000 & nCount_RNA > 200 & nCount_RNA < 2500)
```
  
visualize distribution of filtered data 
```{r}
VlnPlot(pm16,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
VlnPlot(pm26,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
VlnPlot(pm40,features=c("nFeature_RNA","nCount_RNA","percent.rRNA"),ncol=3)
```
  
3. merge individual data sets into one item and split into list of individual experiments
  -cell IDs (orig.ident attribute) assigned based on time point of sampling 
  -each sample separated into item in list based on these cell IDs 
```{r}
allPm.combo<-merge(pm40,y=c(pm16,pm26),add.cell.ids=c("pm40","pm16","pm26"),project="PmNuclei")
allPm.list<-SplitObject(allPm.combo, split.by="orig.ident")
allPm.list
```
  
##Normalization and integration of Samples 
1. normalize each dataset using the SCTransform function, regressing out percent.RNA
```{r}
for (i in 1:length(allPm.list)) {
    allPm.list[[i]] <- SCTransform(allPm.list[[i]], vars.to.regress = c("percent.rRNA"), verbose = TRUE)
}
```
  
2. Integrate the different timepoints into shared space 
  -2a. prep for integration 
```{r}
features <- SelectIntegrationFeatures(object.list = allPm.list, nfeatures = 3000)
allPm.list <- PrepSCTIntegration(object.list = allPm.list, anchor.features = features)
```
  -2b find anchors
```{r}
allPm.anchors <- FindIntegrationAnchors(object.list = allPm.list, normalization.method = "SCT", anchor.features = features)
```

  -2c. integrate datasets
```{r}
allPm.combined.sct <- IntegrateData(anchorset = allPm.anchors, normalization.method = "SCT")
```

##Dimensional reduction and clustering
1. run PCA
```{r}
allPm.combined.sct<-RunPCA(allPm.combined.sct)
```
2. visualize PCA
```{r}
ElbowPlot(allPm.combined.sct,ndims = 50)
```
  -from this, 20 dimensions were selected for further analysis

3. find knn neighbors 
```{r}
allPm.combined.sct <- FindNeighbors(allPm.combined.sct, dims = 1:20)
```

4. find clusters (louvian) 
  -resolution of 2 used 
```{r}
allPm.combined.sct <- FindClusters(allPm.combined.sct, resolution = 2)
```

5. UMAP reduction 
```{r}
allPm.combined.sct <- RunUMAP(allPm.combined.sct, dims = 1:20)
```
  
visualize clustering 
```{r}
DimPlot(allPm.combined.sct, reduction = "umap") + DimPlot(allPm.combined.sct, reduction = "umap", group.by="orig.ident") 
```

##Marker gene determination
identified positive markrs expressed in at least 5% of nuclei in a cluster wiht at least a 0.25 log fold changre in expression relative to rest of data
```{r}
allPm.combined.sct.markers <-FindAllMarkers(allPm.combined.sct, only.pos = TRUE, min.pct = 0.05, logfc.threshold = 0.25)
allPm.combined.sct.markers 
```

##load file of final P. miniata atlas with named clusters
loading in file provided in supplement of paper 
```{r}
pm<-readRDS("../pminiataDevAtlas.rds")
DimPlot(pm)
```
**this pminiata file will be used in the intgration 


-------------------------------------------------
#S.purpuratus atlas creation 


## open raw data files
read in raw cellranger output (raw datea from Foster et al 2020)
```{r}
hatchBlast <- readRDS(file ="Sp_hatchedBlastula.rds")
hatchBlast
mesBlast <- readRDS(file ="Sp_mesenchymeBlastula.rds")
mesBlast
earlyGast <- readRDS(file ="Sp_earlyGastrula.rds")
earlyGast
```
  

##quality control 

1. quantify the percentage of transcripts in each file that are rRNAs
```{r}
rRNA<-c("LOC590670",'LOC119727991','LOC581711','LOC119719604','LOC119719678',"LOC594000","LOC752098","LOC105445410","LOC119729016",'LOC119730276',"LOC582829","LOC119733178","LOC105447583","LOC582250","LOC592405","LOC589889","LOC119738685","LOC119741035",'LOC587799','LOC119719424') 

hatchBlast[["percent.rRNA"]] <- PercentageFeatureSet(hatchBlast, features = intersect(rownames(hatchBlast),rRNA))
mesBlast[["percent.rRNA"]] <- PercentageFeatureSet(mesBlast, features = intersect(rownames(mesBlast),rRNA))
earlyGast[["percent.rRNA"]] <- PercentageFeatureSet(earlyGast, features = intersect(rownames(earlyGast),rRNA))
```

2. quantify the percentae of transcripts in each file that originate form the mitochondria
```{r}
mtgenes <- c("CYTB","ND6","ND5", "ND4","ND3","ATP6","ATP8","COX2","ND4L","COX1","ND2","ND1")

hatchBlast[["percent.mt"]] <- PercentageFeatureSet(hatchBlast, features = intersect(rownames(hatchBlast),mtgenes))
mesBlast[["percent.mt"]] <- PercentageFeatureSet(mesBlast, features = intersect(rownames(mesBlast),mtgenes))
earlyGast[["percent.mt"]] <- PercentageFeatureSet(earlyGast, features = intersect(rownames(earlyGast),mtgenes))
```

3. subset data to apply determined cutoffs for nFeature_RNA, nCount_RNA, and percenrt.rRNA
```{r}
hatchBlast <-subset(hatchBlast, nFeature_RNA < 2500 &  nFeature_RNA > 400 & percent.mt < 5)
mesBlast <-subset(mesBlast,nFeature_RNA > 200 & nFeature_RNA < 1000 & percent.mt < 5) 
earlyGast <-subset(earlyGast,nFeature_RNA > 200 & nFeature_RNA < 2000 & percent.mt < 5 & nCount_RNA < 3000)
```
  
visualize distribution of filtered data 
```{r}
VlnPlot(hatchBlast,features=c("nFeature_RNA","nCount_RNA","percent.rRNA","percent.mt"),ncol=4)
VlnPlot(mesBlast,features=c("nFeature_RNA","nCount_RNA","percent.rRNA","percent.mt"),ncol=4)
VlnPlot(earlyGast,features=c("nFeature_RNA","nCount_RNA","percent.rRNA","percent.mt"),ncol=4)
```

combine into list 
```{r}
allSp.combo <-merge(hatchBlast,y=c(earlyGast,mesBlast))
allSp.list <-SplitObject(allSp.combo,split.by="orig.ident")
allSp.list
```
##normalization and integration of samples
1. normalize each dataset using the SCTransform function, regressing out percent.RNA and percent.mt
```{r}
for (i in 1:length(allSp.list)) {
    allSp.list[[i]] <- SCTransform(allSp.list[[i]], vars.to.regress = c("percent.mt","percent.rRNA"), verbose = TRUE)
}
```

  
2. Integrate the different timepoints into shared space 
  -2a. prep for integration 
```{r}
features <- SelectIntegrationFeatures(object.list = allSp.list, nfeatures = 3000)
allSp.list <- PrepSCTIntegration(object.list = allSp.list, anchor.features = features)
```
  -2b find anchors
```{r}
allSp.anchors <- FindIntegrationAnchors(object.list = allSp.list, normalization.method = "SCT", anchor.features = features)
```

  -2c. integrate datasets
```{r}
allSp.combined.sct <- IntegrateData(anchorset = allSp.anchors, normalization.method = "SCT")
```

##dimensional reduction and clustering

1. run PCA
```{r}
allSp.combined.sct <- RunPCA(allSp.combined.sct, verbose = FALSE)
```
2. visualize PCA
```{r}
ElbowPlot(allSp.combined.sct,ndims=50)
```
  -from this, 20 dimensions were selected for further analysis

3. find knn neighbors 
```{r}
allSp.combined.sct <- FindNeighbors(allSp.combined.sct, dims = 1:20)
```

4. find clusters (louvian) 
  -resolution of 1 used 
```{r}
allSp.combined.sct <- FindClusters(allSp.combined.sct, resolution = 1)
```

5. UMAP reduction 
```{r}
allSp.combined.sct <- RunUMAP(allSp.combined.sct, reduction = "pca", dims = 1:20)
```
  
visualize clustering 
```{r}
DimPlot(allSp.combined.sct, reduction = "umap",label=TRUE) + NoLegend()+  DimPlot(allSp.combined.sct, reduction = "umap", group.by = "orig.ident")
```


##marker gene determination

identified positive markrs expressed in at least 5% of nuclei in a cluster wiht at least a 0.25 log fold changre in expression relative to rest of data
```{r}
allSp.combined.sct.markers <-FindAllMarkers(allSp.combined.sct, only.pos = TRUE, min.pct = 0.05, logfc.threshold = 0.25)
allSp.combined.sct.markers 
```


-------------------------------------------------
#Integration of S.purpuratus and P. miniata atlases 

##load file the named atlases to be compared 
loading in file provided in supplement of paper 
**note the gene symbols of 1:1 orthologs were renamed to the gene symbol used by S. pupuratus in the P. miniata atlas
```{r}
sp<-readRDS("../spDevAtlas.rds")
pm<-readRDS("../pmDevAtlas.rds")
```


##format atlases
add original cluster name attribute and species attribute
```{r}
sp$orig.clust<-Idents(sp)
sp$species<-"Sp"
pm$orig.clust <- Idents(pm)
pm$species<-"Pm"
```


merge files and clear previous normalization (will re-run normalization with SCtransform)
```{r}
pm.sp<-merge(pm,sp)

DefaultAssay(pm.sp)<-"RNA"
pm.sp[["integrated"]]<-NULL
pm.sp[["SCT"]]<-NULL

pm.sp.list <- SplitObject(pm.sp, split.by="orig.ident")
pm.sp.list
```


##Normalization and integration of Samples 
1. normalize each dataset using the SCTransform function, regressing out percent.RNA
```{r}
for (i in 1:length(pm.sp)) {
    pm.sp[[i]] <- SCTransform(pm.sp[[i]], vars.to.regress = c("percent.rRNA","percent.mt"), verbose = TRUE)
}
```
  
2. Integrate the different timepoints into shared space 
  -2a. prep for integration 
```{r}
features<-SelectIntegrationFeatures(object.list = pm.sp.list, nfeatures=3000)
pm.sp.list<-PrepSCTIntegration(object.list=pm.sp.list, anchor.features = features)
```
  -2b find anchors
```{r}
pm.sp.anchors<-FindIntegrationAnchors(object.list = pm.sp.list,normalization.method = "SCT",anchor.features = features)
```

  -2c. integrate datasets
```{r}
pm.sp.int <-IntegrateData(anchorset = pm.sp.anchors,normalization.method = "SCT")
```

##dimensional reductino and clustering

1. run PCA
```{r}
pm.sp.int<-RunPCA(pm.sp.int)
```
2. visualize PCA
```{r}
ElbowPlot(pm.sp.int,ndims=50)
```
  -from this, 10 dimensions were selected for further analysis

3. find knn neighbors 
```{r}
pm.sp.int<-FindNeighbors(pm.sp.int, dims=1:10)
```

4. find clusters (louvian) 
  -resolution of 1 used 
```{r}
pm.sp.int<-FindClusters(pm.sp.int, resolution=1)
```

5. UMAP reduction 
```{r}
pm.sp.int<-RunUMAP(pm.sp.int,dims= 1:10)
```
  
visualize clustering 
```{r}
DimPlot(pm.sp.int,label=TRUE) + DimPlot(pm.sp.int,group.by="species") + DimPlot(pm.sp.int, group.by="orig.clust",label=TRUE)+NoLegend()
```


