---
title: "GranulomaFinalSoupX"
output: html_document
---

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


## Format SoupX 


```{r, }

FormatSoupXInput <- function(RawPrefix,WritePrefix,BarcodesName){
  BarcodesLoadName <- paste(BarcodesName,".Rdata",sep="")
  load(BarcodesLoadName)
  #FilteredName <- paste(FilteredPrefix,".umi.dge.txt.gz",sep="")
  #Filtered <- read.table(FilteredName,header=T,row.names=1)
 
  Load50KName <- paste(RawPrefix,"_50K.Rdata",sep="")

  loadRData <- function(fileName){
#loads an RData file, and returns it
    load(fileName)
    get(ls()[ls() != "fileName"])
}
Raw <- loadRData(Load50KName)

  Filtered <- Raw[,Barcodes]
  
  BarcodesFilt <- paste(Barcodes,"1",sep="-")
  BarcodesFilt <- as.data.frame(BarcodesFilt)
  colnames(BarcodesFilt) <- NULL
  write.table(BarcodesFilt,file=paste(WritePrefix,"/filtered_gene_bc_matrices/GRCh38/barcodes.tsv",sep=""),quote=FALSE,sep='\t',row.names=FALSE)
  NamesRaw <- colnames(Raw)
  BarcodesRaw <- paste(NamesRaw,"1",sep="-")
  BarcodesRaw <- as.data.frame(BarcodesRaw)
  colnames(BarcodesRaw) <- NULL
  write.table(BarcodesRaw,file=paste(WritePrefix,"/raw_gene_bc_matrices/GRCh38/barcodes.tsv",sep=""),quote=FALSE,sep='\t',row.names=FALSE)
  
  #Genes.tsv
  ## Column1 = ENSG_IDS ; Column2 = Gene Names
  nGenes <- dim(Filtered)[1]
  ENSGNames <- paste(rep("ENSG0000",nGenes),c(1:nGenes),sep="")
  GeneNames <- row.names(Filtered)
  GeneTSV <- as.data.frame(cbind(ENSGNames,GeneNames))
  colnames(GeneTSV) <- NULL
  row.names(GeneTSV) <- NULL
  
  write.table(GeneTSV,file=paste(WritePrefix,"/filtered_gene_bc_matrices/GRCh38/genes.tsv",sep=""),quote=FALSE,sep='\t',row.names = FALSE)
  write.table(GeneTSV,file=paste(WritePrefix,"/raw_gene_bc_matrices/GRCh38/genes.tsv",sep=""),quote=FALSE,sep='\t',row.names = FALSE)
  
  
  #matrix.mtx
  Raw <- Raw[GeneNames,]
  colnames(Raw) <- NULL
  row.names(Raw) <- NULL
  colnames(Filtered) <- NULL
  row.names(Filtered) <- NULL
  RawSparse <- Matrix(as.matrix(Raw), sparse = TRUE)
  RawSparse <- as(RawSparse, "dgTMatrix")
  FiltSparse <- Matrix(as.matrix(Filtered), sparse = TRUE)
  FiltSparse <- as(FiltSparse, "dgTMatrix")
  
  writeMM(FiltSparse,paste(WritePrefix,"/filtered_gene_bc_matrices/GRCh38/matrix.mtx",sep=""))
  writeMM(RawSparse,paste(WritePrefix,"/raw_gene_bc_matrices/GRCh38/matrix.mtx",sep=""))
}



FormatSoupXInput("Array1","Array1","Array1_Barcodes")

```

### Round 1

```{r, }
library(SoupX)
devtools::install_github("constantAmateur/SoupX")

StirThePotRound1 <- function(dataDirs,UMICutoff,Name,Animal){
  scl = load10X(dataDirs,soupRange = c(0, UMICutoff))
  scl = inferNonExpressedGenes(scl)
  tstGenes = rownames(scl$channels$Channel1$nonExpressedGenes)[seq(20)]
  gg = plotMarkerDistribution(scl, "Channel1", tstGenes)
  MarkerDistributionPlotName <- paste(Animal,"/",Name,"_MarkerDistributionPlot.pdf",sep="")
  pdf(file=MarkerDistributionPlotName,height=5,width=10)
  print(plot(gg))
  SoupMatrix <- as.data.frame(scl$soupMatrix)
  SoupMatrix$Genes <- row.names(SoupMatrix)
  SoupMatrix <- SoupMatrix[order(SoupMatrix[,1],decreasing = TRUE),]
  SoupMatrixName <- paste(Animal,"/",Name,"_SoupMatrix.Rdata",sep="")
  save(file=SoupMatrixName,SoupMatrix)
  dev.off()
}

StirThePotRound1("Array1/",100,"Array1","Animal")
```



### Stir the Pot Round 2

```{r,}


StirThePotRound2 <- function(dataDirs,UMICutoff,Name,CorrectList,Animal){
  scl = load10X(dataDirs,soupRange = c(0, UMICutoff))
  scl = inferNonExpressedGenes(scl)
  tstGenes = rownames(scl$channels$Channel1$nonExpressedGenes)[seq(20)]
  gg = plotMarkerDistribution(scl, "Channel1", tstGenes)
  MarkerDistributionPlotName <- paste(Name,"MarkerDistributionPlot.pdf",sep="_")
  pdf(file=MarkerDistributionPlotName,height=5,width=10)
  print(plot(gg))
  dev.off()
  SoupMatrix <- as.data.frame(scl$soupMatrix)
  SoupMatrix$Genes <- row.names(SoupMatrix)
  SoupMatrix <- SoupMatrix[order(SoupMatrix[,1],decreasing = TRUE),]
  SoupMatrixName <- paste(Name,"SoupMatrix.Rdata",sep="_")
  save(file=SoupMatrixName,SoupMatrix)
  scl = calculateContaminationFraction(scl, "Channel1", CorrectList)
  gg = plotChannelContamination(scl, "Channel1")
  ContaminationFractionName <- paste(Name,"ContaminationFraction.pdf",sep="_")
  pdf(file=ContaminationFractionName,height=5,width=10)
  plot(gg)
  dev.off()
  scl = interpolateCellContamination(scl, "Channel1", useGlobal = TRUE)
  InterpolatedContaminationFractionName <- paste(Name,"InterpolatedContaminationFraction.pdf",sep="_")
  gg = plotChannelContamination(scl, "Channel1")
  pdf(file=InterpolatedContaminationFractionName,height=5,width=10)
  plot(gg)
  dev.off()
  head(scl$channels$Channel1$rhos)
  scl = strainCells(scl)
  scl = adjustCounts(scl)
  atoc <- as.matrix(scl$atoc)
  atoc <- as.data.frame(atoc)
  atoc_names <- colnames(atoc)
  
  atocName <- paste(Animal,"/",Name,"_atoc.Rdata",sep="")
  save(file=atocName,atoc)
  SoupXName <- paste(Animal,"/",Name,"_SoupX.Rdata",sep="")
  save(file=SoupXName,scl)
  #SeuratObject = createCleanedSeurat(scl)
  #SeuratName <- paste(Animal,"/",Name,"_Seurat.Rdata",sep="")
  #save(file=SeuratName,SeuratObject)
}


StirThePotRound2("Array1/",100,"Array1",c(T2P=c("SFTPC","SFTPB"),Plasma=c("JCHAIN"),Mast=c("CPA3")),"Animal")


### Doublet Finder Final 
```{r, }
devtools::install_github('chris-mcginnis-ucsf/DoubletFinder')

library(DoubletFinder)

FindDoublets <- function(Data,Name,doubletRate,Animal){
  load(Data)
  NameUse <- paste(Name,"atoc",sep="_")
  SeuratObject <- CreateSeuratObject(counts = eval(parse(text=NameUse))) 
  SeuratObject <- NormalizeData(object = SeuratObject)
  SeuratObject <- ScaleData(object = SeuratObject)
  SeuratObject <- FindVariableFeatures(object = SeuratObject, do.plot = FALSE)
  SeuratObject = RunPCA(SeuratObject,pc.genes=SeuratObject@var.genes)
  SeuratObject = RunTSNE(object=SeuratObject,dims.use=1:20,do.fast=TRUE,perplexity=40)
  
  nCells <- dim(SeuratObject@data)[2]
  nDoublets <- floor(nCells*doubletRate)
  nDoublet5 <- floor(nCells*0.05)
  nDoublet10 <- floor(nCells*0.10)
  nDoublet15 <- floor(nCells*0.15)
  nDoublet20 <- floor(nCells*0.20)
  nDoublet25 <- floor(nCells*0.25)
  
  SeuratObject <- doubletFinder(SeuratObject, expected.doublets = nDoublets, proportion.artificial = 0.25, proportion.NN = 0.01)
  SeuratObject@meta.data$Doublet[SeuratObject@meta.data$pANNPredictions == "Singlet"] <- 0
  SeuratObject@meta.data$Doublet[SeuratObject@meta.data$pANNPredictions == "Doublet"] <- 1
  HistogramName <- paste(Name,"pANN_Histogram.pdf",sep="_")
  pdf(file=HistogramName,height=7,width=8)
  hist(SeuratObject@meta.data$pANN,breaks=50)
  dev.off()
  DoubletList <- as.data.frame(SeuratObject@meta.data)
  DoubletSubset <- DoubletList[DoubletList$Doublet == 1,]
  pANNCutoff <- min(DoubletSubset$pANN)
  SeuratObject@meta.data$Doublet[SeuratObject@meta.data$pANN >= pANNCutoff] <- 1
  DoubletList <- as.data.frame(SeuratObject@meta.data)
  
  DoubletList <- DoubletList[order(-DoubletList$pANN),]
  pANNCutoff5 <- DoubletList$pANN[floor(nCells*0.05)]
  pANNCutoff10 <- DoubletList$pANN[floor(nCells*0.10)]
  pANNCutoff15 <- DoubletList$pANN[floor(nCells*0.15)]
  pANNCutoff20 <- DoubletList$pANN[floor(nCells*0.20)]
  pANNCutoff25 <-   DoubletList$pANN[floor(nCells*0.25)]
    
  SeuratObject@meta.data$Doublet5 <- 0
  SeuratObject@meta.data$Doublet10 <- 0
  SeuratObject@meta.data$Doublet15 <- 0
  SeuratObject@meta.data$Doublet20 <- 0
  SeuratObject@meta.data$Doublet25 <- 0
  
  SeuratObject@meta.data$Doublet5[SeuratObject@meta.data$pANN >= pANNCutoff5] <- 1
  SeuratObject@meta.data$Doublet10[SeuratObject@meta.data$pANN >= pANNCutoff10] <- 1
  SeuratObject@meta.data$Doublet15[SeuratObject@meta.data$pANN >= pANNCutoff15] <- 1
  SeuratObject@meta.data$Doublet20[SeuratObject@meta.data$pANN >= pANNCutoff20] <- 1
  SeuratObject@meta.data$Doublet25[SeuratObject@meta.data$pANN >= pANNCutoff25] <- 1

  #SeuratName <- paste(Animal,Name,"DoubletSeurat.Rdata",sep="_")
  #save(file=SeuratName,SeuratObject)
  DoubletListName <- paste(Animal,"/",Name,"_DoubletList.Rdata",sep="")
  save(file=DoubletListName,DoubletList)
  DoubletPlotName <- paste(Animal,"/",Name,"_DoubletPlot.pdf",sep="")
  pdf(file=DoubletPlotName,width=11,height=5)
  FeaturePlot(SeuratObject,c('pANN','Doublet'))
  dev.off()
  
  DoubletProgressionPlotName <- paste(Animal,"/",Name,"_DoubletProgresionPlot.pdf",sep="")
  pdf(file=DoubletProgressionPlotName,width=10,height=15)
  FeaturePlot(SeuratObject,c('pANN','Doublet5','Doublet10','Doublet15','Doublet20','Doublet25'))
  dev.off()
}

#### Run doubletFinder 4017 #### 
FindDoublets("Array1_atoc.Rdata","Array1",0.0139,"Animal") #15000 cells
```