GabrielHoffman

# Run code from 
# https://diseaseneurogenomics.github.io/zenith/reference/zenithPR_gsa.html
# 
# Then run DESeq() and zenithPR_gsa()

# object construction
dds <- DESeqDataSetFromMatrix(geneCounts[keep,], df_metadata, ~ gender)

# standard analysis
dds <- DESeq(dds)

GabrielHoffman

library(dreamlet)
library(muscat)
library(mashr)
library(SingleCellExperiment)

data(example_sce)

# create pseudobulk for each sample and cell cluster
pb <- aggregateToPseudoBulk(example_sce[1:100, ],

GabrielHoffman

# Plot within and between correlations for a continuous property
plot_within_btw = function( data, info){

	library(tidyverse)
	library(corrr)

	if( ! all(c("ID", "Property") %in% colnames(info)) ){
		stop("info must have column names ID and Property")
	}

GabrielHoffman

# Example code for analysis of 
# 1) Trait specificity 
# 2) Cell-type specificity   


# 1) Trait specificity    
######################

library(arrow)

GabrielHoffman

#' Remove variance partition component and recompute fractions
#' 
#' @param vp result of \code{fitVarPart()}
#' @param exclude array of column names to exclude 
removeVarPartComponent = function(vp, exclude){

	# get column index of exclude
	j = match(exclude, colnames(vp)[-c(1:2)])

	# compute new variance fractions

GabrielHoffman

# x is days, y is arbitrary scale
n = 100
x = 10 + 1:n
y = scale(x) + 10*scale(x^2) + rnorm(n)

# relationship is non-linear
# but we have stoing knowledge that it *must* be monotonic 
plot(x, y)

GabrielHoffman

library(zellkonverter)
library(SingleCellExperiment)

file = "240118_PsychAD_FULL_1494_subclass_EN_NF_hvg3k_umap.h5ad"
sce = readH5AD(file, use_hdf5=TRUE, verbose=TRUE, raw=TRUE)

# normalized counts
assay(sce, "X")[1:3, 1:3]

# see "alternative expression", i.e. raw counts

GabrielHoffman

Gene expression Y
Genotype matrix X

TODO
1) Compute sd(Y[i,]) for gene i
2) Compute var(X[,j]) for variant j
3) beta = z / sqrt(v*(n+z^2))
	sd = 1 / sqrt(v*(n+z^2))
	varbeta = sd^2
	where v = var(X[,j])

GabrielHoffman

library(mashr)
library(corrplot)
library(remaCor)
library(ggplot2)

load("/sc/arion/projects/roussp01a/sanan/230706_NPSAD_filePrep/mashr/ckpts/mashr_ckpt2.RData")

# how many genes are non-zero in each tissue?
ngenes = apply(betaTab, 2, function(x) sum(x!=0))

GabrielHoffman

Age = 1:100

dsgn = ns(Age,df=2)
rownames(dsgn) = as.character(Age)
dsgn['66',]

dsgn = ns(Age[-c(1:5)],df=2)
rownames(dsgn) = as.character(Age[-c(1:5)])
dsgn['66',]