### Install SpiecEasi using devtools in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Installs the SpiecEasi package from BiocManager in R. This is the recommended method for stable releases. Ensure BiocManager is installed first.

```r
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")

BiocManager::install("SpiecEasi")
```

--------------------------------

### Install gfortran for OSX Users

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Installs the gfortran package on OSX systems, which is a prerequisite for compiling SpiecEasi source code. This command should be run in the terminal.

```sh
xcode-select --install
```

--------------------------------

### Install SpiecEasi Development Version using remotes

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Installs the development version of SpiecEasi from GitHub using the remotes package. This is useful for accessing the latest features or bug fixes. Ensure remotes is installed first.

```r
if (!require("remotes", quietly = TRUE)) install.packages("remotes")

remotes::install_github("zdk123/SpiecEasi")
```

--------------------------------

### Load SpiecEasi Package in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Loads the SpiecEasi package into the current R session, making its functions available for use. This command should be run after successful installation.

```r
library(SpiecEasi)
```

--------------------------------

### Install SpiecEasi using Conda

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Installs the SpiecEasi package via conda, a popular package and environment manager. This command ensures that SpiecEasi is up-to-date with the main branch of the repository.

```sh
conda install -c bioconda r-spieceasi
```

--------------------------------

### Network Inference with SPIEC-EASI (MB Method)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Performs basic network inference using the Meinshausen-Bühlmann neighborhood selection (MB) method. It loads example data, runs the spiec.easi function, extracts the optimal network, and calculates network statistics. Dependencies include the SpiecEasi package.

```r
library(SpiecEasi)

# Load example American Gut Project data
data(amgut1.filt)

# Basic network inference using MB method
se <- spiec.easi(amgut1.filt, method='mb',
                 lambda.min.ratio=1e-2,
                 nlambda=20,
                 pulsar.params=list(rep.num=50))

# Extract the optimal network adjacency matrix
network <- getRefit(se)

# Get network statistics
num_edges <- sum(network)/2
num_nodes <- ncol(network)
cat(sprintf("Network has %d nodes and %d edges\n", num_nodes, num_edges))

# Get the optimal lambda value and stability score
optimal_lambda <- getOptLambda(se)
stability <- getStability(se)
cat(sprintf("Optimal lambda: %.4f, Stability: %.4f\n", optimal_lambda, stability))
```

--------------------------------

### Working with phyloseq Objects in SpiecEasi

Source: https://context7.com/zdk123/spieceasi/llms.txt

Demonstrates seamless integration of SpiecEasi with phyloseq objects for microbiome network inference. This example shows running spiec.easi directly on a phyloseq object and creating an igraph visualization with taxa names from the phyloseq object. Requires SpiecEasi, phyloseq, and igraph packages.

```r
library(phyloseq)

# Load phyloseq data
data('amgut2.filt.phy')

# Run SPIEC-EASI directly on phyloseq object
se.phylo <- spiec.easi(amgut2.filt.phy,
                       method='mb',
                       lambda.min.ratio=1e-2,
                       nlambda=20,
                       pulsar.params=list(rep.num=50))

# Create igraph with taxa names from phyloseq
ig <- adj2igraph(getRefit(se.phylo),
                 vertex.attr=list(name=taxa_names(amgut2.filt.phy)))

# Use phyloseq plotting functions
plot_network(ig, amgut2.filt.phy, type='taxa', color="Rank3")
```

--------------------------------

### Get Network and Model Parameters from spieceasi Object

Source: https://context7.com/zdk123/spieceasi/llms.txt

Extracts key network structures and optimal parameters from a fitted spieceasi object (se). This includes the adjacency matrix, lambda index and value, stability scores, merge matrix, and method-specific matrices like beta for MB or covariance/inverse covariance for glasso. Dependencies include the spieceasi object itself.

```r
# Get optimal network (adjacency matrix)
adj_matrix <- getRefit(se)

# Get optimal lambda index and value
opt_idx <- getOptInd(se)
opt_lambda <- getOptLambda(se)

# Get stability score at optimal lambda
stability <- getStability(se)

# Get edge-wise stability matrix
merge_matrix <- getOptMerge(se)

# For MB method: get coefficient matrix
if (se$est$method == "mb") {
  beta_matrix <- getOptBeta(se)

  # Symmetrize beta matrix
  sym_beta_ave <- symBeta(beta_matrix, mode='ave')
  sym_beta_max <- symBeta(beta_matrix, mode='maxabs')
}

# For glasso method: get covariance matrices
if (se$est$method == "glasso") {
  icov <- getOptiCov(se)  # Inverse covariance
  cov <- getOptCov(se)    # Covariance
}
```

--------------------------------

### Sparse Plus Low-Rank Network Inference with SpiecEasi

Source: https://context7.com/zdk123/spieceasi/llms.txt

Learns networks with latent variable effects using the sparse plus low-rank (slr) decomposition method in SpiecEasi. This example fits the SLR model, extracts the low-rank residual component, and performs robust PCA to identify latent factors. Dependencies include SpiecEasi.

```r
# Fit sparse + low rank model with 10 latent components
se.slr <- spiec.easi(amgut1.filt,
                     method='slr',
                     r=10,
                     lambda.min.ratio=1e-2,
                     nlambda=20,
                     pulsar.params=list(rep.num=20, ncores=4))

# Extract the low-rank residual component
X <- se.slr$est$data
L <- se.slr$est$resid[[getOptInd(se.slr)]]

# Perform robust PCA to get latent factors
pca <- robustPCA(X, L)

# Examine latent factor scores (samples x factors)
dim(pca$scores)
head(pca$scores)

# Check correlation with metadata
# age_cor <- cor(age_vector, pca$scores, use='pairwise')
```

--------------------------------

### Troubleshooting Empty Networks with StARS (Initial)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Illustrates a common issue where StARS selection results in an empty network. It shows how to diagnose this by examining the optimal lambda index and refitted network size, often indicating that the lambda path needs adjustment.

```r
pargs <- list(seed=10010)
se <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=5e-1, nlambda=10, pulsar.params=pargs)
# Warning in .starsind(est$summary, stars.thresh): Optimal lambda may be outside the supplied values
getOptInd(se)
# [1] 10
sum(getRefit(se))/2
# [1] 58
```

--------------------------------

### Troubleshooting Empty Networks with StARS (Further Refinement)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Further refines the troubleshooting process by increasing 'nlambda' after adjusting 'lambda.min.ratio'. This allows for a more granular sampling of the lambda path, leading to a more accurate network stability estimate and a denser network.

```r
se <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-1, nlambda=100, pulsar.params=pargs)
getStability(se)
# [1] 0.04946882
sum(getRefit(se))/2
# [1] 210
```

--------------------------------

### Run SpiecEasi Pipeline (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Executes the core SpiecEasi pipeline for sparse inverse covariance estimation and model selection. This involves applying data transformations and fitting the model using the 'mb' method with a specified lambda minimum ratio and number of lambda values.

```r
se <- spiec.easi(X, method='mb', lambda.min.ratio=1e-2, nlambda=15)
# Applying data transformations...
# Selecting model with pulsar using stars...
# Fitting final estimate with mb...
# done
```

--------------------------------

### Troubleshooting Empty Networks with StARS (Adjusted Lambda Min Ratio)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Provides a solution to the empty network problem by lowering the 'lambda.min.ratio' to explore denser networks. This adjustment helps ensure that the stability selection process finds a suitable lambda value.

```r
se <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-1, nlambda=10, pulsar.params=pargs)
```

--------------------------------

### Visualize Ecological Networks with igraph (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Loads the 'igraph' library, which is typically used for visualizing and analyzing network structures. This snippet assumes that network objects (like `ig.mb`, `ig.gl`, `ig.sparcc` from previous steps) have already been created and are ready for plotting.

```r
library(igraph)
```

--------------------------------

### Comparing Network Inference Methods

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Compares the results of different network inference methods (serial vs. parallel, StARS vs. bstars) by checking the identity of refitted networks and comparing computation times. This helps in validating the efficiency and consistency of the methods.

```r
## serial vs parallel
identical(getRefit(se1), getRefit(se2))
t1[3] > t2[3]
## stars vs bstars
identical(getRefit(se1), getRefit(se3))
t1[3] > t3[3]
identical(getRefit(se2), getRefit(se4))
t2[3] > t4[3]
```

--------------------------------

### Parallel Network Inference with StARS

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Demonstrates parallel execution of SpiecEasi using the StARS selection criterion with multiple cores. This speeds up the computation by distributing the task across available processors.

```r
pargs2 <- list(rep.num=50, seed=10010, ncores=4)
t2 <- system.time(
se2 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
              sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs2)
)
```

--------------------------------

### Compare Degree Distributions of Inferred Networks (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code snippet calculates and plots the degree distributions for three inferred networks: glasso, MB, and sparcc. It computes the degree distribution for each network and then plots them on the same graph for comparison, with a legend indicating which color corresponds to which method. The y-axis represents frequency and the x-axis represents degree.

```r
dd.gl     <- degree.distribution(ig.gl)
dd.mb     <- degree.distribution(ig.mb)
dd.sparcc <- degree.distribution(ig.sparcc)

plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
      ylab="Frequency", xlab="Degree", main="Degree Distributions")
points(0:(length(dd.gl)-1), dd.gl, col="red" , type='b')
points(0:(length(dd.mb)-1), dd.mb, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
        col=c("forestgreen", "red", "black"), pch=1, lty=1)
```

--------------------------------

### Perform parallel model selection with pulsar in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code illustrates how to use the pulsar package as a backend for SpiecEasi's model selection, specifically using the StARS criterion. It shows how to set `pulsar.select=TRUE` and pass parallel computation parameters (like `rep.num` and `ncores`) via `pulsar.params`. Setting a seed ensures reproducible results. The `system.time` function is used to measure the execution time.

```r
## Default settings ##
pargs1 <- list(rep.num=50, seed=10010)
t1 <- system.time(
se1 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
              sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs1)
)
```

--------------------------------

### Batch Network Inference with batchtools

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Configures and runs SpiecEasi in batch mode using the batchtools package for distributed computing. This allows independent execution of subsamples on an HPC cluster.

```r
## bargs <- list(rep.num=50, seed=10010, conffile="path/to/conf.txt")
bargs <- list(rep.num=50, seed=10010, conffile="parallel")
## See the config file stores:
pulsar::findConfFile('parallel')

## uncomment line below to turn off batchtools reporting
# options(batchtools.verbose=FALSE)
se5 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
            sel.criterion='stars', pulsar.select='batch', pulsar.params=bargs)
```

--------------------------------

### Network Inference with Bounded-StARS

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Implements network inference using the bounded-StARS (bstars) method, which offers computational savings by restricting the lambda path search. This is particularly useful for large datasets.

```r
t3 <- system.time(
se3 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
               sel.criterion='bstars', pulsar.select=TRUE, pulsar.params=pargs1)
)
t4 <- system.time(
se4 <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-3, nlambda=30,
               sel.criterion='bstars', pulsar.select=TRUE, pulsar.params=pargs2)
)
```

--------------------------------

### Fit cross-domain interaction networks in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code demonstrates fitting a network model using SpiecEasi for cross-domain interaction data, such as combining 16S rRNA and Proteomics data. It takes a list of phyloseq objects as input and uses the 'mb' method. The `pulsar.params` argument can be used to control parallel processing. The resulting network can be visualized, coloring vertices based on their data domain.

```r
library(phyloseq)
data(hmp2)
se.hmp2 <- spiec.easi(list(hmp216S, hmp2prot), method='mb', nlambda=40,
              lambda.min.ratio=1e-2, pulsar.params = list(thresh = 0.05))

dtype <- c(rep(1,ntaxa(hmp216S)), rep(2,ntaxa(hmp2prot)))
plot(adj2igraph(getRefit(se.hmp2)), vertex.color=dtype+1, vertex.size=9)
```

--------------------------------

### Evaluate and Compare Network Edge Weights (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This code evaluates and compares the edge weights of networks inferred by SparCC, SpiecEasi MB, and glasso. It calculates correlations from covariance matrices, extracts beta coefficients, summarizes edge lists, and plots histograms of the edge weights for visual comparison. Dependencies include the Matrix package and SpiecEasi functions.

```r
library(Matrix)
secor  <- cov2cor(getOptCov(se.gl.amgut))
sebeta <- symBeta(getOptBeta(se.mb.amgut), mode='maxabs')
elist.gl     <- summary(triu(secor*getRefit(se.gl.amgut), k=1))
elist.mb     <- summary(sebeta)
elist.sparcc <- summary(sparcc.graph*sparcc.amgut$Cor)

hist(elist.sparcc[,3], main='', xlab='edge weights')
hist(elist.mb[,3], add=TRUE, col='forestgreen')
hist(elist.gl[,3], add=TRUE, col='red')
```

--------------------------------

### Visualize Network with Vertex Size Proportional to Mean CLR (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This code snippet visualizes three networks (MB, glasso, sparcc) using the Fruchterman-Reingold layout. The size of each vertex is set to be proportional to the mean centered log-ratio (clr) of the corresponding row, plus an offset of 6. Vertex labels are suppressed for clarity.

```r
vsize    <- rowMeans(clr(amgut1.filt, 1))+6
am.coord <- layout.fruchterman.reingold(ig.mb)

par(mfrow=c(1,3))
plot(ig.mb, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="MB")
plot(ig.gl, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="sparcc")
```

--------------------------------

### Parallel and Batch Processing for spieceasi

Source: https://context7.com/zdk123/spieceasi/llms.txt

Illustrates how to accelerate the spieceasi inference process using parallel processing, bounded StARS criterion, and batch mode for High-Performance Computing (HPC) clusters. It compares serial processing with multicore parallel execution and shows how to use 'bstars' for faster selection and 'batch' mode for cluster environments. Key parameters include 'ncores', 'sel.criterion', and 'pulsar.select'.

```r
# Serial processing (default)
# Note: Ensure amgut1.filt is defined before this code.
pargs1 <- list(rep.num=50, seed=10010)
se_serial <- spiec.easi(amgut1.filt, method='mb',
                        lambda.min.ratio=1e-3, nlambda=30,
                        pulsar.params=pargs1)

# Parallel multicore processing
pargs2 <- list(rep.num=50, seed=10010, ncores=4)
se_parallel <- spiec.easi(amgut1.filt, method='mb',
                          lambda.min.ratio=1e-3, nlambda=30,
                          pulsar.params=pargs2)

# Bounded StARS for faster inference
se_bstars <- spiec.easi(amgut1.filt, method='mb',
                        lambda.min.ratio=1e-3, nlambda=30,
                        sel.criterion='bstars',
                        pulsar.params=pargs2)

# Batch mode for HPC clusters
bargs <- list(rep.num=50, seed=10010, conffile="parallel")
se_batch <- spiec.easi(amgut1.filt, method='mb',
                       lambda.min.ratio=1e-3, nlambda=30,
                       pulsar.select='batch',
                       pulsar.params=bargs)
```

--------------------------------

### Load and Normalize Data for SpiecEasi (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Loads the American Gut dataset, calculates sequencing depths, and normalizes counts to a common scale. This prepares the data for subsequent steps in the SpiecEasi pipeline, such as fitting count margins and generating synthetic networks.

```r
data(amgut1.filt)
depths <- rowSums(amgut1.filt)
amgut1.filt.n  <- t(apply(amgut1.filt, 1, norm_to_total))
amgut1.filt.cs <- round(amgut1.filt.n * min(depths))

d <- ncol(amgut1.filt.cs)
n <- nrow(amgut1.filt.cs)
e <- d
```

--------------------------------

### Diversity Metrics Calculation in R

Source: https://context7.com/zdk123/spieceasi/llms.txt

Demonstrates the calculation of ecological diversity indices from sample count data using functions from the spieceasi package (or related ecological packages). It covers total sum normalization, Shannon entropy calculation, and the computation of the effective number of species (exp(Shannon)). Requires a vector of sample counts as input.

```r
# Sample counts
counts <- c(10, 20, 15, 5, 30)

# Total sum normalization
normalized <- norm_to_total(counts)

# Shannon entropy
shannon_val <- shannon(counts)

# Effective number of species (exponential of Shannon)
n_eff <- neff(counts)
```

--------------------------------

### Generate Diverse Network Topologies (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Generates various network topologies including Erdos-Renyi random graphs, scale-free graphs, hub graphs, cluster graphs, band graphs, and block diagonal graphs. It also demonstrates converting a graph to a precision matrix with specified limits for positive and negative theta values and a target condition number. Requires the 'igraph' package.

```r
# Erdos-Renyi random graph
g_random <- make_graph("erdos_renyi", D=50, e=75)

# Scale-free graph
g_scalefree <- make_graph("scale_free", D=50, e=75)

# Hub graph with 3 hub nodes
g_hub <- make_graph("hub", D=50, e=75, numHubs=3)

# Cluster graph with community structure
g_cluster <- make_graph("cluster", D=50, e=75)

# Band (chain-like) graph
g_band <- make_graph("band", D=50, e=75)

# Block diagonal graph
g_block <- make_graph("block", D=50, e=75, numHubs=3)

# Convert any graph to precision matrix
prec_random <- graph2prec(g_random,
                          posThetaLims=c(0.2, 0.5),
                          negThetaLims=c(-0.5, -0.2),
                          targetCondition=100)

```

--------------------------------

### Analyze American Gut Data with SpiecEasi (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Applies SpiecEasi directly to the American Gut data, comparing results from the 'mb' (Meinshausen-Bühlmann) and 'glasso' (graphical lasso) methods, as well as SparCC (correlation-based). Normalization is handled internally by `spiec.easi`. It also demonstrates passing parameters to the internal StARS selection function and utilizing multiple cores.

```r
se.mb.amgut <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=50))
se.gl.amgut <- spiec.easi(amgut1.filt, method='glasso', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=50))
sparcc.amgut <- sparcc(amgut1.filt)
## Define arbitrary threshold for SparCC correlation matrix for the graph
sparcc.graph <- abs(sparcc.amgut$Cor) >= 0.3
diag(sparcc.graph) <- 0
library(Matrix)
sparcc.graph <- Matrix(sparcc.graph, sparse=TRUE)
## Create igraph objects
ig.mb     <- adj2igraph(getRefit(se.mb.amgut))
ig.gl     <- adj2igraph(getRefit(se.gl.amgut))
ig.sparcc <- adj2igraph(sparcc.graph)
```

--------------------------------

### Network Analysis with Phyloseq Object (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code demonstrates using SpiecEasi with a phyloseq object. It loads the 'amgut2.filt.phy' dataset, runs SpiecEasi using the 'mb' method with specified lambda parameters, converts the resulting adjacency matrix to an igraph object, and then plots the network colored by 'Rank3' taxa. This requires the phyloseq package.

```r
library(phyloseq)
## Load round 2 of American gut project
data('amgut2.filt.phy')
se.mb.amgut2 <- spiec.easi(amgut2.filt.phy, method='mb', lambda.min.ratio=1e-2,
                           nlambda=20, pulsar.params=list(rep.num=50))
ig2.mb <- adj2igraph(getRefit(se.mb.amgut2),  vertex.attr=list(name=taxa_names(amgut2.filt.phy)))
plot_network(ig2.mb, amgut2.filt.phy, type='taxa', color="Rank3")
```

--------------------------------

### Network Visualization with igraph in R

Source: https://context7.com/zdk123/spieceasi/llms.txt

Demonstrates how to convert a network inferred by spiec.easi into an igraph object for visualization and analysis. It includes fitting the model, converting the adjacency matrix to an igraph object, calculating node properties, setting a layout, plotting the network, and computing basic network statistics like density and average degree. Requires the 'igraph' and 'spieceasi' packages.

```r
library(igraph)

# Fit model
se <- spiec.easi(amgut1.filt, method='mb',
                 lambda.min.ratio=1e-2, nlambda=20,
                 pulsar.params=list(rep.num=50))

# Convert to igraph
ig <- adj2igraph(getRefit(se),
                 vertex.attr=list(name=colnames(amgut1.filt)))

# Calculate node properties
vsize <- rowMeans(clr(amgut1.filt+1, 1)) + 6
degree_vals <- degree(ig)

# Layout
coords <- layout.fruchterman.reingold(ig)

# Plot
plot(ig, layout=coords,
     vertex.size=vsize,
     vertex.label=NA,
     vertex.color="lightblue",
     edge.width=1,
     edge.color="gray50")

# Network statistics
cat(sprintf("Nodes: %d, Edges: %d\n", vcount(ig), ecount(ig)))
cat(sprintf("Density: %.4f\n", edge_density(ig)))
cat(sprintf("Average degree: %.2f\n", mean(degree(ig))))
```

--------------------------------

### Compute SparCC Correlation and Build Network (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Computes SparCC correlations and covariance from microbial abundance data, then constructs a network by thresholding the absolute correlation values. The resulting network is converted to a sparse matrix and an igraph object for further analysis. It also compares the edge count with SPIEC-EASI. Requires the 'Matrix' library.

```r
# Compute SparCC correlations
sparcc_result <- sparcc(amgut1.filt, iter=20, inner_iter=10, th=0.1)

# Extract correlation and covariance matrices
cor_matrix <- sparcc_result$Cor
cov_matrix <- sparcc_result$Cov

# Create network by thresholding correlations
threshold <- 0.3
sparcc_network <- abs(cor_matrix) >= threshold
diag(sparcc_network) <- 0

# Convert to sparse matrix and igraph
library(Matrix)
sparcc_network <- Matrix(sparcc_network, sparse=TRUE)
ig_sparcc <- adj2igraph(sparcc_network)

# Compare edge count with SPIEC-EASI
cat(sprintf("SparCC edges: %d\n", sum(sparcc_network)/2))

```

--------------------------------

### Network Extraction from Fitted Model (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Fits a model using the SPIEC-EASI function, likely for network inference from microbial abundance data. It uses the 'mb' (Meinshausen-Bühlmann) method and specifies parameters for lambda selection, including lambda minimum ratio and the number of lambdas. It also includes parameters for pulsar, suggesting ensemble methods. The estimated network is stored in `se$est$path`.

```r
# Fit a model
se <- spiec.easi(amgut1.filt, method='mb',
                 lambda.min.ratio=1e-2, nlambda=20,
                 pulsar.params=list(rep.num=50))

```

--------------------------------

### Synthesize Ecological Network Data (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Generates a synthetic ecological network and associated data using SpiecEasi. It creates a graph, calculates the precision and correlation matrices, and then synthesizes community counts from a zero-inflated negative binomial distribution based on the generated correlation structure.

```r
set.seed(10010)
graph <- SpiecEasi::make_graph('cluster', d, e)
Prec  <- graph2prec(graph)
Cor   <- cov2cor(prec2cov(Prec))

X <- synth_comm_from_counts(amgut1.filt.cs, mar=2, distr='zinegbin', Sigma=Cor, n=n)
```

--------------------------------

### Construct robust PCA from SpiecEasi output in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code constructs a robust Principal Component Analysis (PCA) model using the estimated data and residuals from a SPIEC-EASI SLR network. It extracts the data (X) and residuals (L) from the `spiec.easi` object and then applies the `robustPCA` function. The resulting `pca` object contains scores, loadings, and other PCA components.

```r
X <- se2.slr.amgut$est$data
L <- se2.slr.amgut$est$resid[[getOptInd(se2.slr.amgut)]]
pca <- robustPCA(X, L)
```

--------------------------------

### Apply Compositional Data Transformations (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Demonstrates various transformations for compositional microbiome data, including Centered Log-Ratio (CLR), Total Sum Normalization, and Additive Log-Ratio (ALR). It shows how to apply these transformations to both matrix and data frame objects, specifying the margin (rows or columns) for transformations. Requires custom functions like 'clr', 'norm_to_total', 'norm_pseudo', and 'alr'.

```r
# Sample count data
counts <- matrix(rpois(50, lambda=10), nrow=10, ncol=5)

# Centered log-ratio (CLR) transformation
clr_data <- clr(counts, mar=2)  # mar=2 for columns

# Total sum normalization
normalized <- t(apply(counts, 1, norm_to_total))

# Add pseudocount before normalization
pseudo_norm <- t(apply(counts, 1, norm_pseudo))

# Additive log-ratio (ALR) transformation
alr_data <- alr(counts, mar=2, divcomp=1)

# For data frames
counts_df <- as.data.frame(counts)
clr_df <- clr(counts_df)
alr_df <- alr(counts_df)

```

--------------------------------

### Multi-Domain Network Inference with SpiecEasi

Source: https://context7.com/zdk123/spieceasi/llms.txt

Infers joint microbial association networks across multiple data types (domains) using SpiecEasi. This snippet demonstrates combining 16S and proteomics data, running spiec.easi for inference, and visualizing the resulting cross-domain network. Requires SpiecEasi and igraph packages.

```r
# Load multi-domain data (16S and proteomics)
data(hmp2)

# Combine 16S and proteomics data in a list
datalist <- list(hmp216S, hmp2prot)

# Infer joint network across domains
se.multi <- spiec.easi(datalist,
                       method='mb',
                       nlambda=40,
                       lambda.min.ratio=1e-2,
                       pulsar.params=list(thresh=0.05, rep.num=50))

# Extract network and identify cross-domain edges
network <- getRefit(se.multi)
n_16S <- ntaxa(hmp216S)
n_prot <- ntaxa(hmp2prot)

# Visualize with domain labels
library(igraph)
dtype <- c(rep(1, n_16S), rep(2, n_prot))
ig <- adj2igraph(network)
plot(ig, vertex.color=dtype+1, vertex.size=9,
     vertex.label=NA, edge.width=0.5)
```

--------------------------------

### Matrix Conversion Utilities in R

Source: https://context7.com/zdk123/spieceasi/llms.txt

Provides functions to convert between different types of matrices: precision (inverse covariance), covariance, and correlation. It demonstrates how to transform a precision matrix to covariance, covariance back to precision, covariance to correlation, and correlation to covariance using specified standard deviations. This is useful for various statistical modeling and data analysis tasks.

```r
# Create a precision (inverse covariance) matrix
prec <- matrix(c(2.0, -0.5, 0.0,
                 -0.5, 2.0, -0.5,
                 0.0, -0.5, 2.0), nrow=3)

# Convert precision to covariance
cov_mat <- prec2cov(prec)

# Convert covariance back to precision
prec_recovered <- cov2prec(cov_mat)

# Convert covariance to correlation
cor_mat <- cov2cor(cov_mat)

# Convert correlation to covariance with specific SDs
sds <- c(2, 3, 4)
cov_from_cor <- cor2cov(cor_mat, sds)

# Verify calculations
all.equal(prec, prec_recovered)
```

--------------------------------

### Bootstrap SparCC for P-values and Significance (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Performs bootstrap inference on microbial abundance data to calculate p-values for SparCC correlations. It allows for 'both' sided testing and extracts correlations with p-values below a specified threshold (e.g., 0.05). Requires the 'sparccboot' and 'pval.sparccboot' functions.

```r
# Perform bootstrap inference
set.seed(123)
boot_result <- sparccboot(amgut1.filt,
                          R=100,
                          ncpus=4)

# Calculate p-values
pval_result <- pval.sparccboot(boot_result, sided='both')

# Extract significant correlations
significant_cors <- pval_result$cors[pval_result$pvals < 0.05]
significant_pvals <- pval_result$pvals[pval_result$pvals < 0.05]

# Display results
length(significant_cors)
summary(significant_cors)

```

--------------------------------

### Network Inference with glasso Method in SpiecEasi

Source: https://context7.com/zdk123/spieceasi/llms.txt

Estimates microbial association networks using the graphical LASSO (glasso) method. This snippet demonstrates running spiec.easi with the 'glasso' method, extracting covariance and inverse covariance matrices, and deriving edge weights from the correlation structure. Requires the SpiecEasi and Matrix packages.

```r
# Run glasso method for network inference
se.gl <- spiec.easi(amgut1.filt,
                    method='glasso',
                    lambda.min.ratio=1e-2,
                    nlambda=20,
                    pulsar.params=list(rep.num=50, ncores=4))

# Extract covariance and inverse covariance matrices
icov_matrix <- getOptiCov(se.gl)   # Inverse covariance (precision)
cov_matrix <- getOptCov(se.gl)     # Covariance matrix

# Get edge weights from the correlation structure
library(Matrix)
secor <- cov2cor(cov_matrix)
network <- getRefit(se.gl)
edge_list <- summary(triu(secor * network, k=1))

# Display edge information
head(edge_list)
```

--------------------------------

### Generate Synthetic Microbiome Data with SPIEC-EASI (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Generates synthetic microbiome count data using a known network structure, fitting marginal distributions from real data. It normalizes counts, creates a graph topology, converts it to a precision matrix, and then generates synthetic data with a Zero-Inflated Negative Binomial (ZINB) distribution. Finally, it runs SPIEC-EASI on the synthetic data and evaluates reconstruction accuracy using 'huge'. Requires 'Matrix', 'igraph', and 'huge' libraries.

```r
# Load real data to fit marginal distributions
data(amgut1.filt)

# Normalize to common depth
depths <- rowSums(amgut1.filt)
amgut1.filt.n <- t(apply(amgut1.filt, 1, norm_to_total))
amgut1.filt.cs <- round(amgut1.filt.n * min(depths))

# Create synthetic network structure
set.seed(10010)
d <- ncol(amgut1.filt.cs)
n <- nrow(amgut1.filt.cs)
e <- d  # number of edges

# Generate cluster graph topology
graph <- make_graph('cluster', d, e)

# Convert graph to precision matrix
Prec <- graph2prec(graph)

# Convert precision to correlation
Cor <- cov2cor(prec2cov(Prec))

# Generate synthetic community with ZINB distribution
X_synth <- synth_comm_from_counts(amgut1.filt.cs,
                                  mar=2,
                                  distr='zinegbin',
                                  Sigma=Cor,
                                  n=n)

# Verify dimensions
dim(X_synth)

# Run SPIEC-EASI on synthetic data
se_synth <- spiec.easi(X_synth, method='mb',
                       lambda.min.ratio=1e-2,
                       nlambda=15)

# Evaluate reconstruction accuracy
library(huge)
roc_result <- huge::huge.roc(se_synth$est$path, graph, verbose=FALSE)

```

--------------------------------

### Generate Multivariate Count Data (R)

Source: https://context7.com/zdk123/spieceasi/llms.txt

Generates multivariate count data from various distributions including Zero-inflated Poisson (ZIP), Zero-inflated Negative Binomial (ZINB), Poisson (Pois), and Negative Binomial (NB), while maintaining a specified correlation structure (Sigma). It also generates multivariate normal data. The function `cor()` is used to check the empirical correlations of the generated data. Requires functions like 'rmvzipois', 'rmvzinegbin', 'rmvpois', 'rmvnegbin', and 'rmvnorm'.

```r
# Define correlation structure
Sigma <- matrix(c(1.0, 0.5, 0.3,
                  0.5, 1.0, 0.4,
                  0.3, 0.4, 1.0), nrow=3)

# Zero-inflated Poisson
mu_zip <- c(5, 10, 15)
data_zip <- rmvzipois(n=100, mu=mu_zip, Sigma=Sigma)

# Zero-inflated Negative Binomial
mu_zinb <- c(5, 10, 15)
data_zinb <- rmvzinegbin(n=100, mu=mu_zinb, Sigma=Sigma)

# Poisson
mu_pois <- c(5, 10, 15)
data_pois <- rmvpois(n=100, mu=mu_pois, Sigma=Sigma)

# Negative Binomial
mu_nb <- c(5, 10, 15)
data_nb <- rmvnegbin(n=100, mu=mu_nb, Sigma=Sigma)

# Multivariate Normal
mu_norm <- c(0, 0, 0)
data_norm <- rmvnorm(n=100, mu=mu_norm, Sigma=Sigma)

# Check empirical correlations
cor(data_zinb)

```

--------------------------------

### Evaluate SpiecEasi Model Performance (R)

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

Evaluates the performance of the fitted SpiecEasi model by examining the Receiver Operating Characteristic (ROC) curve over the lambda path and the Precision-Recall (PR) curve for the selected graph using the StARS criterion. It uses functions from the 'huge' and 'SpiecEasi' packages.

```r
huge::huge.roc(se$est$path, graph, verbose=FALSE) stars.pr(getOptMerge(se), graph, verbose=FALSE)
# stars selected final network under: getRefit(se)
```

--------------------------------

### Fit networks with and without latent effects (SPIEC-EASI slr) in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code fits two types of networks using the spiec.easi function: 'mb' for standard model-based networks and 'slr' for sparse + low-rank decomposition to remove latent effects. It demonstrates how to specify parameters like lambda.min.ratio, nlambda, and pulsar.params for parallel processing. The output can be used to compare network consistency.

```r
se1.mb.amgut <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=20, ncores=4))
se2.mb.amgut <- spiec.easi(amgut2.filt.phy, method='mb', lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=20, ncores=4))


se1.slr.amgut <- spiec.easi(amgut1.filt, method='slr', r=10, lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=20, ncores=4))
se2.slr.amgut <- spiec.easi(amgut2.filt.phy, method='slr', r=10, lambda.min.ratio=1e-2,
                          nlambda=20, pulsar.params=list(rep.num=20, ncores=4))
```

--------------------------------

### Correlate metadata with latent factors in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code calculates the correlation between ecological metadata (age, BMI, depth) and the latent factors derived from robust PCA. It extracts relevant metadata from a phyloseq object and computes pairwise correlations using the `cor` function with the `use='pairwise'` option to handle missing values. This helps understand the biological relevance of the identified latent variables.

```r
age <- as.numeric(as.character(sample_data(amgut2.filt.phy)$AGE))mi <- as.numeric(as.character(sample_data(amgut2.filt.phy)$BMI))
depth <- colSums(otu_table(amgut2.filt.phy))

cor(age, pca$scores, use='pairwise')
cor(bmi, pca$scores, use='pairwise')
cor(depth, pca$scores, use='pairwise')
```

--------------------------------

### Compare edge set consistency using Jaccard index in R

Source: https://github.com/zdk123/spieceasi/blob/master/README.md

This R code snippet calculates the Jaccard index to compare the consistency of edge sets between two fitted networks. It uses the edge.diss function, which requires the refitted network objects obtained from spiec.easi, the desired metric ('jaccard'), and the taxon names. This is useful for evaluating the improvement gained by removing latent effects.

```r
otu1 <- colnames(amgut1.filt)
otu2 <- taxa_names(amgut2.filt.phy)
edge.diss(getRefit(se1.mb.amgut), getRefit(se2.mb.amgut), 'jaccard', otu1, otu2)
edge.diss(getRefit(se1.slr.amgut), getRefit(se2.slr.amgut), 'jaccard', otu1, otu2)
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.