gene plot ideas

gene_plot_ideas.R

# 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")

gene <- "FBgn0001098"

load_all("~/bioc/fishpond/github/fishpond")
plotAllelicGene(y, gene, db=txdb, genome="dm6")

gene_plot_mouse.R

library(SummarizedExperiment)

library(AnnotationHub)
library(ensembldb)
ah <- AnnotationHub()
#query(ah, c("EnsDb","102","Mus musculus"))
edb <- ah[["AH89211"]]

## finding genes...

library(Gviz)
library(plyranges)

load("../osteoblast-quant/data/tss_se_filtered.rda")
y <- gse[,gse$cross == "CASTxB6"]
devtools::load_all("~/bioc/fishpond/github/fishpond")
y <- labelKeep(y)
y <- swish(y, x="allele", pair="day")
hist(mcols(y)$pvalue)
y <- computeInfRV(y)
save(y, file="tss_results_se.rda")
load("tss_results_se.rda")

library(dplyr)
library(tidyr)
library(tibble)
tss <- mcols(y) %>%
  as.data.frame() %>%
  rownames_to_column("id") %>%
  tibble()

tss %>% filter(qvalue < .01) %>%
  arrange(pvalue, -abs(log2FC)) %>%
  select(-tx_id, -keep) %>%
  write.csv(file="tss_global.csv", row.names=FALSE)

load("gene_results.rda")

gene_dat <- gene %>%
  select(symbol, gene_id = id, gene_lfc = log2FC, gene_q = qvalue)

dat <- tss %>%
  select(id, gene_id, log10mean, tss_lfc = log2FC, tss_q = qvalue) %>%
  inner_join(gene_dat, by="gene_id") %>%
  filter(tss_q < .01 & gene_q < .01) %>%
  filter(abs(gene_lfc) < 7 & abs(tss_lfc) < 7)

dat2 <- dat %>%
  filter(sign(gene_lfc) != sign(tss_lfc) &
         abs(gene_lfc) > .5 & abs(tss_lfc) > 1 &
         log10mean > 1)
  
library(ggrepel)
dat %>%
  filter(gene_lfc > 0 & tss_lfc < 0) %>%
  ggplot(aes(gene_lfc, tss_lfc, color=log10mean)) + 
  geom_point() +
  scale_color_gradient(low="blue",high="yellow") +
  geom_label_repel(data=dat2[dat2$gene_lfc>0,],
                   aes(gene_lfc, tss_lfc, label=symbol)) +
  xlab("Gene-level aLFC") + ylab("TSS-level aLFC")

dat %>%
  filter(gene_lfc < 0 & tss_lfc > 0) %>%
  ggplot(aes(gene_lfc, tss_lfc, color=log10mean)) + 
  geom_point() +
  scale_color_gradient(low="blue",high="yellow") +
  geom_label_repel(data=dat2[dat2$gene_lfc<0,],
                   aes(gene_lfc, tss_lfc, label=symbol)) +
  xlab("Gene-level aLFC") + ylab("TSS-level aLFC")

write.csv(dat2, file="discordant_global.csv", row.names=FALSE)

levels(y$allele) <- c("B6","CAST")
gene <- "Nfatc3"
(dat3 <- dat %>% filter(symbol == gene) %>% select(id, log10mean, tss_lfc))
par(mfrow=c(3,1))
for (i in 1:nrow(dat3)) {
  plotInfReps(y, dat3$id[i], x="day", cov="allele", legend=TRUE, legendPos="topleft")
}

### moving plot to fishpond

## load("../osteoblast-quant/data/tss_se_filtered.rda")
## y <- gse[,gse$cross == "CASTxB6"]
## y <- labelKeep(y)
## y <- swish(y, x="allele", pair="day")
## hist(mcols(y)$pvalue)
## y <- computeInfRV(y)

load("tss_results_se.rda")
load_all("~/bioc/fishpond/github/fishpond")

gene <- "ENSMUSG00000026193"
gene <- "ENSMUSG00000031902"
plotAllelicGene(y, gene=gene, db=edb)
plotAllelicGene(y, gene=gene, db=edb, tpmFilter=10)
plotAllelicGene(y, gene=gene, db=edb, countFilter=500)

plotAllelicGene(y, symbol="Nfatc3", db=edb, tpmFilter=.01, countFilter=1, isoPropFilter=.01)
plotAllelicGene(y, symbol="Npr3", db=edb)

dat %>% filter(symbol == "Igfbp4") %>% select(id, log10mean, tss_lfc)
gene <- "ENSMUSG00000017493"
plotAllelicGene(y, gene=gene, db=edb)
plotInfReps(y, idx="ENSMUSG00000017493-99043597", x="day", cov="allele")
# this without tpmFilter needs to be de-bugged

dat %>% filter(symbol == "Fuca2") %>% select(id, log10mean, tss_lfc)
gene <- "ENSMUSG00000019810"
plotAllelicGene(y, gene=gene, db=edb)

sym <- "Zmat3"
sym <- "Csnk1a1"
sym <- "Eno1"
sym <- "Malat1"
sym <- "Il10rb"
gene <- mcols(y)$gene_id[match(sym, mcols(y)$symbol)]
plotAllelicGene(y, gene=gene, db=edb)

gene_plot_mouse_dynamic.R

library(SummarizedExperiment)
load("../osteoblast-quant/data/tss_se_filtered.rda")
y <- gse[,gse$cross == "CASTxB6"]
load_all("~/bioc/fishpond/github/fishpond")
y <- labelKeep(y)

# filter
infRep1 <- assay(y,"infRep1")
a1 <- y$allele == "a1"
mcols(y)$someInfo <- rowSums(abs(infRep1[,!a1] - infRep1[,a1]) < 1) < ncol(y)/2
table(mcols(y)$someInfo)
y <- y[mcols(y)$someInfo,]

y <- swish(y, x="allele", pair="day", cov="day", cor="pearson")
hist(mcols(y)$pvalue)
y <- computeInfRV(y)
save(y, file="tss_results_se_dynamic.rda")
load("tss_results_se.rda")

library(dplyr)
library(tidyr)
library(tibble)
tss <- mcols(y) %>%
  as.data.frame() %>%
  rownames_to_column("id") %>%
  tibble()

tss %>% filter(qvalue < .1) %>%
  arrange(pvalue, -abs(log2FC)) %>%
  select(-tx_id, -keep) %>%
  write.csv(file="tss_dynamic.csv", row.names=FALSE)

### find discordant

load("gene_results_dynamic.rda")

gene_dat <- gene %>%
  select(symbol, gene_id = id, gene_stat = stat, gene_q = qvalue)

dat <- tss %>%
  select(id, gene_id, log10mean, tss_stat = stat, tss_q = qvalue) %>%
  inner_join(gene_dat, by="gene_id") %>%
  filter(tss_q < .1)

dat2 <- dat %>%
  filter(abs(gene_stat - tss_stat) > .25 &
         log10mean > 1)
  
library(ggrepel)
dat %>%
  ggplot(aes(gene_stat, tss_stat, color=log10mean)) + 
  geom_point() +
  scale_color_gradient(low="blue",high="yellow") +
  geom_label_repel(data=dat2,
                   aes(gene_stat, tss_stat, label=symbol))

write.csv(dat2, file="discordant_dynamic.csv", row.names=FALSE)

sparc.R

library(SummarizedExperiment)
library(ensembldb)
library(AnnotationHub)

ah <- AnnotationHub()
#query(ah, c("EnsDb","102","Mus musculus"))
edb <- ah[["AH89211"]]

# or:
edbfile <- ensDbFromGtf("~/Downloads/Mus_musculus.GRCm38.102.gtf.gz")
edb <- EnsDb(edbfile)

tbg <- transcriptsBy(edb, by="gene")
tbg[["ENSMUSG00000018593"]]

load("tss_results_se.rda")

pdf("~/Desktop/sparc_gene_model.pdf")
r <- GRanges("11", IRanges(start=55.418e6, end=55.422e6))
plotAllelicGene(y, gene="ENSMUSG00000018593", db=edb, region=r,
                labels=list(a1="CAST", a2="B6"),
                tpmFilter=5, qvalue=FALSE, transcriptAnnotation="transcript")
dev.off()

pdf("~/Desktop/sparc_gene_model_2.pdf")
plotAllelicGene(y, gene="ENSMUSG00000018593", db=edb, 
                labels=list(a1="CAST", a2="B6"),
                tpmFilter=5, qvalue=FALSE, transcriptAnnotation="transcript")
dev.off()


library(tibble)
library(dplyr)
mcols(y)[mcols(y)$gene_id == "ENSMUSG00000018593",] %>%
  as_tibble() %>%
  mutate(tx_id = sapply(tx_id, paste, collapse=",")) %>%
  select(tx_id, group_id, log10mean, log2FC, qvalue)

rowMeans(assay(y, "abundance")[mcols(y)$gene_id == "ENSMUSG00000018593",])
rowMeans(assay(y, "counts")[mcols(y)$gene_id == "ENSMUSG00000018593",])

test_gene.R

library(SummarizedExperiment)
load("../osteoblast-quant/data/gse_filtered.rda")
y <- gse[,gse$cross == "CASTxB6"]
load_all("~/bioc/fishpond/github/fishpond")
y <- labelKeep(y)
y <- swish(y, x="allele", pair="day")
hist(mcols(y)$pvalue)

library(tidyr)
library(tibble)
gene <- mcols(y) %>%
  as.data.frame() %>%
  rownames_to_column("id") %>%
  tibble()
save(gene, file="gene_results.rda")

library(dplyr)
gene %>% filter(qvalue < .01) %>%
  arrange(pvalue, -abs(log2FC)) %>%
  select(-keep) %>%
  write.csv(file="gene_global.csv", row.names=FALSE)

### dynamic

library(SummarizedExperiment)
load("../osteoblast-quant/data/gse_filtered.rda")
y <- gse[,gse$cross == "CASTxB6"]
load_all("~/bioc/fishpond/github/fishpond")
y <- labelKeep(y)

# filter
infRep1 <- assay(y,"infRep1")
a1 <- y$allele == "a1"
mcols(y)$someInfo <- rowSums(abs(infRep1[,!a1] - infRep1[,a1]) < 1) < ncol(y)/2
table(mcols(y)$someInfo)
y <- y[mcols(y)$someInfo,]

y <- swish(y, x="allele", pair="day", cov="day", cor="pearson")
hist(mcols(y)$pvalue)

library(tidyr)
library(tibble)
gene <- mcols(y) %>%
  as.data.frame() %>%
  rownames_to_column("id") %>%
  tibble()
save(gene, file="gene_results_dynamic.rda")

library(dplyr)
gene %>% filter(qvalue < .1) %>%
  arrange(pvalue, -abs(log2FC)) %>%
  select(-keep, -someInfo) %>%
  write.csv(file="gene_dynamic.csv", row.names=FALSE)