function NormalizeData(X::SparseMatrixCSC; scale_factor::Int = 10000)::SparseMatrixCSC{Float32,Int32}
    total = sum(X, dims = 1) /scale_factor;
    colind = X.colptr;
    rown, coln = size(X)
    Y = Float32.(X)
    @inbounds @views Threads.@threads for i in 1:coln
            a = colind[i]
            b = colind[i + 1] - 1
            Y.nzval[a:b] = (Float32.(X.nzval[a:b] ./ total[i]))  
        end
    Y.nzval .= log1p.(Y.nzval)        
    return Y
end

function NormalizeData(X::Pissaro; scale_factor::Int = 10000, assay = X.default_assay)::Pissaro
    println("Normalizing ", assay, " with scale factor ", scale_factor)
    X.assays[assay].norm_count.array = NormalizeData(X.assays[assay].count.array, scale_factor = scale_factor) 
    return X 
end

function NormalizeData!(X::Pissaro; scale_factor::Int = 10000, assay = X.default_assay)
    println("Normalizing ", assay, " with scale factor ", scale_factor)
    X.assays[assay].norm_count.array = NormalizeData(X.assays[assay].count.array, scale_factor = scale_factor) 
    nothing 
end

function ScaleData(X::AbstractArray)::Array{Float32, 2}
    scale_data::Array{Float32, 2} = (Matrix(X) .- mean(X, dims = 2)) ./ std(X, dims = 2)
    return scale_data
end

function ScaleData(X::Pissaro; features::Vector{String}, assay::String = X.default_assay)::Pissaro
    X.assays[assay].scaled_count = ScaleData(X.assays[assay].norm_count.array)
    return(X)
end

function ScaleData!(X::Pissaro; assay::String = X.default_assay)
    X.assays[assay].scaled_count.array = ScaleData(X.assays[assay].norm_count.array)
    nothing
end

function PercentageFeatureSet(X; pattern::String, assay::String = X.default_assay)
    genes = names(X.assays[assay].count)[1]
    pat_genes = genes[occursin.(Regex(pattern), genes)]
    percent = sum(X.assays[assay].count[pat_genes,:], dims = 1) ./ sum(X.assays[assay].count, dims = 1)
    return DataFrame(cells = names(percent)[2], percentage = percent.array[1,:]) 
end

function PercentageFeatureSet(X; features::Vector{String}, assay::String = X.default_assay)
    genes = names(X.assays[assay].count)[1]
    if all(intersect(features, genes) .== features)
        percent = sum(X.assays[assay].count[features,:], dims = 1) ./ sum(X.assays[assay].count, dims = 1)
    else
        throw("nope")
    end
    return DataFrame(cells = names(percent)[2], percentage = percent.array[1,:]) 
end
function PercentageFeatureSet(X; features::String, assay::String = X.default_assay)
    X.assay[X.default_assay]
end

features()
["CD79A", "CD79B"]
X = pissaro.assays["RNA"].norm_count.array

@time ScaledData = (Matrix(X) .- mean(X, dims = 2)) ./ std(X, dims = 2)


@time PercentageFeatureSet(pissaro, pattern =  "^RP")
@time PercentageFeatureSet(pissaro, features =  ["CD79A", "CD79B", "DNTT"])

Matrix(X)
pattern = "^MT-"
X = pissaro

