## ----global_options, include=FALSE---------------------------------------
library(knitr)
opts_chunk$set(tidy=TRUE)

## ----quick_start, eval=FALSE---------------------------------------------
#  install.packages('SoupX')
#  library(SoupX)
#  #Load data and estimate soup profile
#  sc = load10X('Path/to/cellranger/outs/folder/')
#  #Estimate rho
#  sc = autoEstCont(sc)
#  #Clean the data
#  out = adjustCounts(sc)

## ----install_CRAN,eval=FALSE---------------------------------------------
#  install.packages('SoupX')

## ----install, eval=FALSE-------------------------------------------------
#  devtools::install_github("constantAmateur/SoupX",ref='devel')

## ----load----------------------------------------------------------------
library(SoupX)

## ----download,eval=FALSE-------------------------------------------------
#  tmpDir = tempdir(check=TRUE)
#  download.file('http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz',destfile=file.path(tmpDir,'tod.tar.gz'))
#  download.file('http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc4k/pbmc4k_filtered_gene_bc_matrices.tar.gz',destfile=file.path(tmpDir,'toc.tar.gz'))
#  untar(file.path(tmpDir,'tod.tar.gz'),exdir=tmpDir)
#  untar(file.path(tmpDir,'toc.tar.gz'),exdir=tmpDir)

## ----load_data,eval=FALSE------------------------------------------------
#  sc = load10X(tmpDir)

## ----load_data_manual,eval=FALSE-----------------------------------------
#  toc = Seurat::Read10X(file.path(tmpDir,'filtered_gene_bc_matrices','GRCh38'))
#  tod = Seurat::Read10X(file.path(tmpDir,'raw_gene_bc_matrices','GRCh38'))
#  sc = SoupChannel(tod,toc)

## ----load_saved----------------------------------------------------------
data(PBMC_sc)
sc = PBMC_sc
sc

## ----estimateSoup, eval=FALSE--------------------------------------------
#  sc = SoupChannel(tod,toc,calcSoupProfile=FALSE)
#  sc = estimateSoup(sc)

## ----estimateNoDrops, eval=TRUE------------------------------------------
library(Matrix)
toc = sc$toc
scNoDrops = SoupChannel(toc,toc,calcSoupProfile=FALSE)
#Calculate soup profile
soupProf = data.frame(row.names = rownames(toc),
                      est = rowSums(toc)/sum(toc),
                      counts = rowSums(toc))
scNoDrops = setSoupProfile(scNoDrops,soupProf)

## ----set_clustering------------------------------------------------------
data(PBMC_metaData)
sc = setClusters(sc,setNames(PBMC_metaData$Cluster,rownames(PBMC_metaData)))

## ----add_DR--------------------------------------------------------------
sc = setDR(sc,PBMC_metaData[colnames(sc$toc),c('RD1','RD2')])

## ----plot_annot----------------------------------------------------------
library(ggplot2)
dd = PBMC_metaData[colnames(sc$toc),]
mids = aggregate(cbind(RD1,RD2) ~ Annotation,data=dd,FUN=mean)
gg = ggplot(dd,aes(RD1,RD2)) + 
  geom_point(aes(colour=Annotation),size=0.2) +
  geom_label(data=mids,aes(label=Annotation)) +
  ggtitle('PBMC 4k Annotation') +
  guides(colour = guide_legend(override.aes = list(size=1)))
plot(gg)

## ----plot_IGKC-----------------------------------------------------------
dd$IGKC = sc$toc['IGKC',]
gg = ggplot(dd,aes(RD1,RD2)) +
  geom_point(aes(colour=IGKC>0))
plot(gg)

## ----sanity_check--------------------------------------------------------
gg = plotMarkerMap(sc,'IGKC')
plot(gg)

## ----set_rho-------------------------------------------------------------
sc = setContaminationFraction(sc,0.2)

## ----genes1--------------------------------------------------------------
nonExpressedGeneList = list(HB=c('HBB','HBA2'),IG = c('IGKC'))

## ----genes2--------------------------------------------------------------
nonExpressedGeneList = list(HB=c('HBB','HBA2'),IG = c('IGKC','IGHG1','IGHG3'))

## ----auto_est------------------------------------------------------------
sc = autoEstCont(sc)

## ----auto_est_unif_prior-------------------------------------------------
sc = autoEstCont(sc,priorRhoStdDev=0.3)

## ----topSoupGenes--------------------------------------------------------
head(sc$soupProfile[order(sc$soupProfile$est,decreasing=TRUE),],n=20)

## ----inferNonExpressed---------------------------------------------------
plotMarkerDistribution(sc)

## ----igGenes-------------------------------------------------------------
igGenes = c('IGHA1','IGHA2','IGHG1','IGHG2','IGHG3','IGHG4','IGHD','IGHE','IGHM',
            'IGLC1','IGLC2','IGLC3','IGLC4','IGLC5','IGLC6','IGLC7',
            'IGKC')

## ----calculateNullMatrix-------------------------------------------------
useToEst = estimateNonExpressingCells(sc,nonExpressedGeneList = list(IG=igGenes),clusters=FALSE)

## ----visNullMatrix-------------------------------------------------------
plotMarkerMap(sc,geneSet=igGenes,useToEst=useToEst)

## ----calcNullMatrixWithClustering----------------------------------------
useToEst = estimateNonExpressingCells(sc,nonExpressedGeneList = list(IG=igGenes))
plotMarkerMap(sc,geneSet=igGenes,useToEst=useToEst)

## ----calcContamination---------------------------------------------------
sc = calculateContaminationFraction(sc,list(IG=igGenes),useToEst=useToEst)

## ----viewCont------------------------------------------------------------
head(sc$metaData)

## ----decontaminate-------------------------------------------------------
out = adjustCounts(sc)

## ----mostZeroed----------------------------------------------------------
cntSoggy = rowSums(sc$toc>0)
cntStrained = rowSums(out>0)
mostZeroed = tail(sort((cntSoggy-cntStrained)/cntSoggy),n=10)
mostZeroed

## ----mostReduced---------------------------------------------------------
tail(sort(rowSums(sc$toc>out)/rowSums(sc$toc>0)),n=20)

## ----IGKC_change---------------------------------------------------------
plotChangeMap(sc,out,'IGKC')

## ----change_plots,eval=FALSE---------------------------------------------
#  plotChangeMap(sc,out,'LYZ')
#  plotChangeMap(sc,out,'CD74')
#  plotChangeMap(sc,out,'IL32')
#  plotChangeMap(sc,out,'TRAC')
#  plotChangeMap(sc,out,'S100A9')
#  plotChangeMap(sc,out,'NKG7')
#  plotChangeMap(sc,out,'GNLY')
#  plotChangeMap(sc,out,'CD4')
#  plotChangeMap(sc,out,'CD8A')

## ----writeOut,eval=FALSE-------------------------------------------------
#  DropletUtils:::write10xCounts('./strainedCounts',out)

