0
$\begingroup$

I am running a public 10x dataset through SCONE in which one of the normalization techniques is from SCARN

#10x scone
suppressPackageStartupMessages({
  library(scater) # BioConductor
  library(SingleCellExperiment) # BioConductor
  library(scRNAseq) #BioConductor
  library(DropletUtils) # BioConductor
  library(tidyverse) # CRAN
  library(here) # CRAN
  library(DT) # CRAN
  library(edgeR) #BioConductor
  library(pheatmap) # CRAN
})

#Load Data
data_dir <- "filtered_feature_bc_matrix"
sce <- read10xCounts(data_dir,col.names=TRUE)
rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ID, rowData(sce)$Symbol)

#Add QC for Mitochondrial genes
is.mito <- grep("^MT-", rownames(sce))
sce <- addPerCellQC(sce, subsets=list(Mito=is.mito))

keep=c("sum","detected",          
         "percent_top_50","percent_top_100","percent_top_200","percent_top_500",      
         "subsets_Mito_sum","subsets_Mito_detected","subsets_Mito_percent","total")

# Adaptive Filter using MAD method
qc.lib2 <- isOutlier(sce$sum, log=TRUE, type="lower") #lib sizes
#Gene Filtering step 1
qc.nexprs2 <- isOutlier(sce$detected, log=TRUE, type="lower") #expressed genes
qc.mito2 <- isOutlier(sce$subsets_Mito_percent, type="higher")
attr(qc.lib2, "thresholds") #check cutoffs
discard2 <- qc.lib2 | qc.nexprs2 | qc.mito2
df2<-DataFrame(LibSize=sum(qc.lib2), NExprs=sum(qc.nexprs2),
               MitoProp=sum(qc.mito2), Total=sum(discard2))
> df2
DataFrame with 1 row and 4 columns
    LibSize    NExprs  MitoProp     Total
  <integer> <integer> <integer> <integer>
1       395       344       412       482

plot(sce$sum, sce$subsets_Mito_percent, log="x", xlab="Total count", ylab='Mitochondrial %')
abline(h=attr(qc.mito2, "thresholds")["higher"], col="red")
filtered <- sce[,!discard2]
is.mito <-which(qc.mito2)
lost <- calculateAverage(counts(sce)[,discard2])
kept <- calculateAverage(counts(sce)[,!discard2])
logged <- cpm(cbind(lost, kept), log=TRUE, prior.count=2)
logFC <- logged[,1] - logged[,2]
abundance <- rowMeans(logged)
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
points(abundance[is.mito], logFC[is.mito], col="dodgerblue", pch=16)

library(scone)
data(housekeeping)

metadata(filtered)$which_qc <- colData(filtered)[keep]
sce.matrix <- as.matrix(counts(filtered))
rownames(sce.matrix) <- rowData(filtered)$Symbol
dim(sce.matrix)
> dim(sce.matrix)
[1] 33538  3514

#Gene Filtering step 2
sce.matrix <- sce.matrix[rowSums(sce.matrix)>0,]
> dim(sce.matrix)
[1] 21195  3514

num_reads <- quantile(sce.matrix[sce.matrix > 0])[4]
num_cells = 0.25*ncol(sce.matrix)
is_common = rowSums(sce.matrix >= num_reads ) >= num_cells
> table(is_common)
is_common
FALSE  TRUE 
20405   790

hk = intersect(housekeeping$V1,rownames(sce.matrix))
> length(hk)
[1] 503

# Metric-based Filtering
#ralign = colData(sce)$detected,
mfilt = metric_sample_filter(sce.matrix,
                             nreads = colSums(sce.matrix),
                             gene_filter = is_common,
                             pos_controls = hk,
                             zcut = 3, mixture = FALSE,
                             plot = TRUE,
                             hard_nreads=2000)

filter_cell <- !apply(simplify2array(mfilt[!is.na(mfilt)]),1,any)

# Final Gene Filtering: Highly expressed in at least 5 cells
num_reads <- quantile(sce.matrix[sce.matrix > 0])[4]
num_cells = 5
is_quality = rowSums(sce.matrix >= num_reads ) >= num_cells

#Gene Filtering step 3
filtered <- sce.matrix[is_quality, filter_cell]
> dim(filtered)
[1] 11105  3434

qc <- colData(sce)[colnames(filtered),]

posi<- c("CD8B",
         "LEF1",
         "OXNAD1",
         "TRABD2A",
         "CCR7",
         "MAL",
         "TCF7",
         "PIK3IP1")

#Making the positive and negative controls
poscon = intersect(rownames(filtered),posi)
negcon = intersect(rownames(filtered),hk)
qc=qc[keep]
ppq = scale(qc[,apply(qc,2,sd) > 0],center = TRUE, scale = TRUE)

my_scone <- SconeExperiment(filtered,
                            qc=ppq, 
                            negcon_ruv = rownames(filtered) %in% negcon,
                            poscon = rownames(filtered) %in% poscon
)

EFF_FN = function (ei)
{
  sums = colSums(ei > 0)
  eo = t(t(ei)*sums/mean(sums))
  return(eo)
}

SCRAN_FN3 = function (ei) 
{
  if (!requireNamespace("scran", quietly = TRUE)) {
    stop("scran package needed for SCRAN_FN()")
  }
  scales = scran::calculateSumFactors(ei) #, sizes = ceiling(sqrt(ncol(ei))))
  eo = t(t(ei) * mean(scales)/scales)
  return(eo)
}

## ----- Scaling Argument -----

scaling1=list(none=identity, # Identity - do nothing
             uq = UQ_FN,
             scran = SCRAN_FN3,
             sum = SUM_FN) #, # SCONE library wrappers...
             #tmm = TMM_FN, 
             #uq = UQ_FN,
             #fq = FQT_FN,
             #deseq = DESEQ_FN
             #eff = EFF_FN, # User-defined function)

my_scone1 <- scone(my_scone,
                  scaling=scaling1,
                  k_qc=3, k_ruv = 3,
                  adjust_bio="no",
                  run=FALSE)

BiocParallel::register(
  BiocParallel::SerialParam()
) # Register BiocParallel Serial Execution

my_scone <- scone(my_scone1,
                  scaling=scaling1,
                  run=TRUE,
                  eval_kclust = 2:6,
                  return_norm = "in_memory",
                  zero = "postadjust")


pc_obj = prcomp(apply(t(get_scores(my_scone)),1,rank),
                center = TRUE,scale = FALSE)
bp_obj = biplot_color(pc_obj,y = -get_score_ranks(my_scone),expand = .6)

I get are following warning messages:

Warning message:
1: In .local(x, ...) :
  none_uq returned at least one NA value. Consider removing it from the comparison.
 In the meantime, failed methods have been filtered from the output.
2: In FUN(...) : encountered negative size factor estimates
3: In FUN(...) : encountered negative size factor estimates

While I do not know why upper quartile technique fails, I have found some useful information regarding negative size factor estimates relating the SCRAN normalization strategy. According to SCRAN's documentation if one has too many negative size factor estimates, then low abundant genes need to be filtered out more stringently. However, I carry out 3 gene filtering steps (See #1 Gene Filtering step, #2 Gene Filtering step and #3 Gene Filtering step in code), so I am not sure how to increase stringency. I have also increased the min.mean value in Scran method, which is another technique they suggest but I still get the same warnings.

Does anyone have insights into this?

$\endgroup$
1
  • $\begingroup$ Warnings are not errors. Warnings do not mean anything failed. $\endgroup$ – swbarnes2 Jul 12 '20 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.