mikelove

library(dplyr)
library(tidyr)
library(purrr)
library(broom)
library(tibble)

n <- 1000
dat <- matrix(runif(6*n),nrow=n)
colnames(dat) <- paste0("parent", rep(c("x","y"),c(3,3)), c(1:3,1:3))
dat <- as_tibble(dat)

mikelove

ngene <- 1000
sig_ratio <- .1
niso <- rpois(ngene, lambda=5)
niso <- pmax(niso, 2)

set.seed(1)

gene <- factor(rep(1:ngene,niso))
pval <- runif(length(gene))
idx <- sapply(split(seq_along(gene), gene)[1:(ngene*sig_ratio)], `[`, 1)

mikelove

data {
  int<lower=0> n;
  int y[n];
  vector[n] x;
}
parameters {
  real<lower=0> mu;
  real<lower=0> delta;
}
transformed parameters {

mikelove

# edited from this function: https://github.com/seandavi/BiocPkgTools/blob/HEAD/R/pkgDependencyMetrics.R#L232
pkgDependencyGain <- function(pkg, depdf, depsToRemove) {

  ## fetch dependency graph
  g <- buildPkgDependencyIgraph(depdf)
  
  ## exclude 'R', 'base' and 'methods'
  excludedpkgs <- c("R", "base", "methods")
  g <- igraph::induced_subgraph(g, setdiff(names(V(g)), excludedpkgs))
  

mikelove

library(samr)
library(limma)
library(DESeq2)

dds <- makeExampleDESeqDataSet(m=100, betaSD=rep(c(0,.5),c(900,100)))

# remove batch effect and re-exponentiate
design_matrix <- model.matrix(design(dds), colData(dds))
batch <- runif(100)
mat <- log2(counts(dds)+1)

mikelove

meanSdPlot <- function(x, ranks = TRUE, xlab = ifelse(ranks, "rank(mean)", "mean"), 
                       ylab = "sd", pch, plot = TRUE, bins = 50, ...) {
      
      stopifnot(is.logical(ranks), length(ranks) == 1, !is.na(ranks))

      n = nrow(x)
      if (n == 0L) {
        warning("In 'meanSdPlot': input matrix 'x' has 0 rows. There is nothing to be done.")
        return()
      }

mikelove

library(mpra)
E <- 1e3 # Number of elements
B <- 3 # Number of barcodes
s <- 4 # Samples per tissue
nt <- 2 # Number of tissues
set.seed(434)
rna <- matrix(rpois(E*B*s*nt, lambda = 2 * 30), nrow = E*B, ncol = s*nt)
dna <- matrix(rpois(E*B*s*nt, lambda = 30), nrow = E*B, ncol = s*nt)
rn <- as.character(outer(paste0("barcode_", seq_len(B), "_"),
                         paste0("elem_", seq_len(E)), FUN = "paste0"))

mikelove

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
g <- genes(txdb)
library(org.Hs.eg.db)
g$symbol <- mapIds(org.Hs.eg.db, g$gene_id, "SYMBOL", "ENTREZID")

library(plyranges)
query <- data.frame(
  seqnames=c("chr1","chr2","chr3"),
  start=50e6, end=50e6) %>%

mikelove

library(magrittr)
library(purrr)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb %>%
  { list(e=exons(.), g=genes(.)) } %>%
  map(GenomicRanges::reduce) %>%
  map(~sum(width(.x))) %>%
  { diff(unname(unlist(.))) / sum(seqlengths(txdb)) }
## [1] 0.4007753

mikelove

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Searching modulated *br* TFBS on DM3 genome using GRAFIMO"
   ]
  },
  {