library(Seurat)

DS_NOV2017_rep1.data <- read.table("./DGEs/12_2017NOV_Sample1_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep1.data <- DS_NOV2017_rep1.data[,colSums(DS_NOV2017_rep1.data) >= 800]
tmp.rep1 <- CreateSeuratObject(DS_NOV2017_rep1.data, min.cells = 1, min.features = 400, project = "rep1")

DS_NOV2017_rep2.data <- read.table("./DGEs/12_2017NOV_Sample2_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep2.data <- DS_NOV2017_rep2.data[,colSums(DS_NOV2017_rep2.data) >= 800]
tmp.rep2 <- CreateSeuratObject(DS_NOV2017_rep2.data, min.cells = 1, min.features = 400, project = "rep2")

DS_NOV2017_rep3.data <- read.table("./DGEs/12_2017NOV_Sample3_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep3.data <- DS_NOV2017_rep3.data[,colSums(DS_NOV2017_rep3.data) >= 800]
tmp.rep3 <- CreateSeuratObject(DS_NOV2017_rep3.data, min.cells = 1, min.features = 400, project = "rep3")

DS_NOV2017_rep4.data <- read.table("./DGEs/12_2017NOV_Sample4_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep4.data <- DS_NOV2017_rep4.data[,colSums(DS_NOV2017_rep4.data) >= 800]
tmp.rep4 <- CreateSeuratObject(DS_NOV2017_rep4.data, min.cells = 1, min.features = 400, project = "rep4")

DS_NOV2017_rep5.data <- read.table("./DGEs/12_2017NOV_Sample5_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep5.data <- DS_NOV2017_rep5.data[,colSums(DS_NOV2017_rep5.data) >= 800]
tmp.rep5 <- CreateSeuratObject(DS_NOV2017_rep5.data, min.cells = 1, min.features = 400, project = "rep5")

DS_NOV2017_rep6.data <- read.table("./DGEs/12_2017NOV_Sample6_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep6.data <- DS_NOV2017_rep6.data[,colSums(DS_NOV2017_rep6.data) >= 800]
tmp.rep6 <- CreateSeuratObject(DS_NOV2017_rep6.data, min.cells = 1, min.features = 400, project = "rep6")

DS_NOV2017_rep7.data <- read.table("./DGEs/12_2017NOV_Sample7_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep7.data <- DS_NOV2017_rep7.data[,colSums(DS_NOV2017_rep7.data) >= 800]
tmp.rep7 <- CreateSeuratObject(DS_NOV2017_rep7.data, min.cells = 1, min.features = 400, project = "rep7")

DS_NOV2017_rep8.data <- read.table("./DGEs/12_2017NOV_Sample8_dmel625_v02_merged_gene_exon_tagged_top10000cells.dge.txt.gz", header = TRUE, sep = "\t", row.names = 1)
DS_NOV2017_rep8.data <- DS_NOV2017_rep8.data[,colSums(DS_NOV2017_rep8.data) >= 800]
tmp.rep8 <- CreateSeuratObject(DS_NOV2017_rep8.data, min.cells = 1, min.features = 400, project = "rep8")

tmp1.merged <- merge(tmp.rep1, tmp.rep2)
tmp2.merged <- merge(tmp1.merged, tmp.rep3)
tmp3.merged <- merge(tmp2.merged, tmp.rep4)
tmp4.merged <- merge(tmp3.merged, tmp.rep5)
tmp5.merged <- merge(tmp4.merged, tmp.rep6)
tmp6.merged <- merge(tmp5.merged, tmp.rep7)
NOV2017_DS_allreps <- merge(tmp6.merged, tmp.rep8)

NOV2017_DS_allreps <- subset(NOV2017_DS_allreps, subset = nFeature_RNA < 30000)

all.features <- as.vector(NOV2017_DS_allreps@assays$RNA@counts@Dimnames[1])
TE.features <- mapply(grep,"^TE-", all.features, value = TRUE)
GENE.features <- mapply(grep, "^TE-", all.features, value = TRUE, invert = TRUE)
NOV2017_DS_allreps.TE.data <- NOV2017_DS_allreps@assays$RNA@counts[TE.features[1:length(TE.features)],]
NOV2017_DS_allreps.GENE.data <- NOV2017_DS_allreps@assays$RNA@counts[GENE.features[1:length(GENE.features)],]

NOV2017_DS_allreps[["GENES"]] <- CreateAssayObject(counts = NOV2017_DS_allreps.GENE.data)
NOV2017_DS_allreps[["TES"]] <- CreateAssayObject(counts = NOV2017_DS_allreps.TE.data)

NOV2017_DS_allreps <- NormalizeData(NOV2017_DS_allreps, assay = "GENES")
NOV2017_DS_allreps <- FindVariableFeatures(NOV2017_DS_allreps, assay = "GENES")
NOV2017_DS_allreps <- ScaleData(NOV2017_DS_allreps, assay = "GENES")
NOV2017_DS_allreps <- RunPCA(NOV2017_DS_allreps, assay = "GENES", npcs = 80, reduction.name = "gene_PCA", reduction.key = "gene_PCA_")
NOV2017_DS_allreps <- RunTSNE(NOV2017_DS_allreps, assay = "GENES", dims = 1:80, reduction = "gene_PCA", reduction.name = "gene_tsne", reduction.key = "gene_tsne_", check_duplicates = F)
NOV2017_DS_allreps <- FindNeighbors(NOV2017_DS_allreps, reduction = "gene_PCA", dims = 1:80, assay = "GENES")
NOV2017_DS_allreps <- FindClusters(NOV2017_DS_allreps, resolution = 3.5, assay = "GENES", graph.name = "GENES_snn")
NOV2017_DS_allreps <- RunUMAP(NOV2017_DS_allreps, reduction = "gene_PCA", dims = 1:80, assay = "GENES", reduction.name = "gene_umap", reduction.key = "gene_umap_")

NOV2017_DS_allreps <- NormalizeData(NOV2017_DS_allreps, assay = "TES")
NOV2017_DS_allreps <- FindVariableFeatures(NOV2017_DS_allreps, assay = "TES")
NOV2017_DS_allreps <- ScaleData(NOV2017_DS_allreps, assay = "TES")
NOV2017_DS_allreps <- RunPCA(NOV2017_DS_allreps, assay = "TES", npcs = 80, reduction.name = "te_PCA", reduction.key = "te_PCA_")
NOV2017_DS_allreps <- RunTSNE(NOV2017_DS_allreps, assay = "TES", dims = 1:80, reduction = "te_PCA", reduction.name = "te_tsne", reduction.key = "te_tsne_", check_duplicates = F)
NOV2017_DS_allreps <- FindNeighbors(NOV2017_DS_allreps, reduction = "te_PCA", dims = 1:80, assay = "TES")
NOV2017_DS_allreps <- FindClusters(NOV2017_DS_allreps, resolution = 3.5, assay = "TES", graph.name = "TES_snn")
NOV2017_DS_allreps <- RunUMAP(NOV2017_DS_allreps, reduction = "te_PCA", dims = 1:80, assay = "TES", reduction.name = "te_umap", reduction.key = "te_umap_")

