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"
   ]
  },
  {

mikelove

library(Biostrings)
dna <- readDNAStringSet("test_sequences.fasta")
ids <- names(dna)
library(tidyverse)
library(plyranges)
dat <- data.frame(ids)
into <- c("seqnames","start","end","pos","relpos",
          "ref","alt","allele","rsid","RE")
dat <- dat %>%
  separate(ids, sep="_", into=into, convert=TRUE)

mikelove

# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2095218
#url <- "ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2095nnn/GSM2095218/suppl/GSM2095218%5FDEX%5F3hr%5FGR%5FChIPseq%2ERep1%5Fpeaks%2Ebed%2Egz"
#download.file(url, destfile="GSM2095218_DEX_3hr_GR_ChIPseq.Rep1_peaks.bed.gz")

library(plyranges)
x <- read_bed("GSM2095218_DEX_3hr_GR_ChIPseq.Rep1_peaks.bed.gz")
genome(x) <- "hg19"

options(Biostrings.coloring=FALSE)

mikelove

x <- read.csv("GSE91061_BMS038109Sample.hg19KnownGene.raw.csv", row.names=1)
condition <- factor(sub(".+_(.+)_.+", "\\1", colnames(x)))

library(DESeq2)
dds <- DESeqDataSetFromMatrix(x, colData=data.frame(condition), ~condition)
vsd <- vst(dds, blind=FALSE)
plotPCA(vsd)

rv <- rowVars(assay(vsd))
pc <- prcomp(t(assay(vsd)[head(order(-rv),1000),]))

mikelove

library(SingleCellExperiment)
sce <- readRDS("sce_metaCells.rds")

dim(sce)
assayNames(sce)

a1 <- as.matrix( round(assay(sce, "dmel_counts")) )
a2 <- as.matrix( round(assay(sce, "counts") - a1) )
a1[is.na(a1)] <- 0
a2[is.na(a2)] <- 0

mikelove

# Variants and peaks coverage
# Michael Love
# March 11, 2022

library(readr)
library(dplyr)
library(tidyr)
library(plyranges)

# read_narrowpeaks -> just works

mikelove

# Gviz library uses TxDb to make gene models
# this is convenient bc we can build TxDb from any source
library(Gviz)
library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm6.ensGene

# load our TSS-summarized data, now has GRanges on rows
suppressPackageStartupMessages(library(SummarizedExperiment))
load("se_tss.rda")

mikelove

library(SingleCellExperiment)
sce <- readRDS("sce.rds")
labels <- readRDS("cellLabels.rds")

sce <- sce[,names(labels)]
colLabels(sce) <- labels

table(labels)
assayNames(sce)
total_count <- rowSums(assay(sce))

mikelove