Spectral Inference

Install

type:

julia-repl ] dev git@github.com:aramanlab/SpectralInference.jl.git into the julia command line

Examples

calculate Spectral distances ```julia using SpectralInference

M = rand(100, 100) usv = svd(M) # from linear algebra spectralpartitions = getintervals(usv.S, alpha=1.0, q=0.5) Dij = spectraldistances(usv.U, usv.S, spectralpartitions) make a tree from the distances using basic hierarchical clusting julia

shortcut of hclust(spimtx; linkage=:average, branchorder=:optimal) from Clustering package

spitree = UPGMA_tree(Dij) ```

or with neighbor joining julia using NeighborJoining # https://github.com/BenjaminDoran/NeighborJoining.jl njclusts = regNJ(Dij) # or fastNJ(Dij) ids = string.(1:100); nwstring = NeighborJoining.newickstring(njclusts, ids)

example script for running SpectralInference on a gene expression matrix in a .csv file to generate spectral tree.

```julia

Example SpectralInference code for gene expression matrix in .csv format (genes = rows, cells = columns)

Ben Doran / Noah Gamble

2023/08/07

Load packages

using SpectralInference using NeighborJoining using CSV, DataFrames using LinearAlgebra using Clustering using Muon using JLD2

Load matrix as data frame from .csv

df = CSV.read("/path/to/geneexpressionmatrix.csv", DataFrame)

Copy to matrix

m = Matrix(df[:,2:end]) #

filter out low abundance genes

filtthresh = 10 genecounts = mapslices(r -> sum(r.>0), m, dims=2) |> vec idxkeep = genecounts .> filtthresh m = m[idxkeep,:]

Organize into annotated data object

m' is transpose of m (make samples rows)

adata = AnnData(X = m') adata.obsnames .= names(df)[2:end] adata.varnames .= string.(df.Column1[idx_keep]) @info size(adata.X)

Perform SVD

@info "Peforming SVD..." @time usv = svd(adata.X)

Perform SpectralInference

@info "Calculating SpectralInference matrix..." @time dij = spectraldistances(usv.U, usv.S, getintervals(usv.S))

Cluster to tree

@info "Clustering SpectralInference matrix..."

@time spitree = UPGMA_tree(dij)

@time nwtreestring = SpectralInference.newickstring(spitree, adata.obs_names.vals)

@time spitree = regNJ(dij) @time nwtreestring = NeighborJoining.newickstring(spitree, adata.obs_names.vals)

Add SpectralInference matrix to adata object

@info "Adding SpectralInference matrix to ADATA..." adata.obsp["spectraldistances"] = dij adata.obs[:,"order"] = spitree.order adata.uns["cellmerges"] = spitree.merges adata.uns["heights"] = spitree.heights adata.uns["treelinkage"] = string(spitree.linkage) adata.uns["newicktreestring"] = nwtreestring

Write annotated data to and Newick tree to disk

@info "Writing results to disk..." @time begin writeh5ad("spitest.h5ad", adata) open("testtree.nw.txt", "w") do io println(io, nwtreestring) end jldsave("test_tree.jld2"; spitree) end

```

example script for creating bootstrap spectral tree and calculating support values

```julia using SpectralInference using NeighborJoining using Muon

spitst = readh5ad("spitest.h5ad") mtx = spitst.X[:, :] NBOOT = 100 boottrees = map(1:NBOOT) do i ncols = size(mtx, 2) tmpmtx = mtx[:, rand(1:ncols, ncols)] usv = svd(tmpmtx) dij = spectraldistances(usv.U, usv.S, getintervals(usv.S)) nwtreestring = NeighborJoining.newickstring(regNJ(dij)) end open("testtree_bootstraps.nw.txt", "w") do io for t in boottrees println(io, t) end end

calculate support values

using GoTreejll reftreefile = joinpath(pwd(), "testtree.nw.txt") boottreefile = joinpath(pwd(), "testtreebootstraps.nw.txt") supporttreefile = joinpath(pwd(), "testtreebootstraps.nw.txt") run($(gotree()) -i $reftreefile -b $boottreefile -o $supporttreefile) ```

script for calculating MI between metavariables and tree depths

```julia using SpectralInference using NewickTree using Muon, DataFrames

spitst = readh5ad("spitest.h5ad") spitree = readnw(readline(joinpath(pwd(), "testtreebootstraps.nw.txt"))) spitree50pct = aspolytomy(n->NewickTree.support(n)<.5, spitree) leafnames = getleafnames(spitree50pct) mi, treedepths = pairedMIacrosstreedepth(eachcolumn(spitst.obs), leafnames, spitree_50pct; ncuts=100)

NBOOTS = 50 results = map(1:NBOOTS) do i pairedMIacrosstreedepth(eachcolumn(spitst.obs), leafnames, spitree_50pct; bootstrap=true, ncuts=100) end

rowmask = spitst.obs.count .> 15 results = map(1:NBOOTS) do i pairedMIacrosstreedepth(eachcolumn(spitst.obs), leafnames, spitree50pct; mask=rowmask, bootstrap=true, ncuts=100) end allmiresults = first.(results) alltreedepths_results = last.(results)

```

Citing

See CITATION.bib for the relevant reference(s).