mikelove

---
title: "Toy tree example for collapsing"
author: "Michael Love"
---
  
Example data with 20 inferential replicates, here we just have 1
sample per condition and we calculate the LFC at each level of the
tree.

From the below simulation setup (see first chunk), the true DE signal

mikelove

library(plyranges)
library(readr)

bindata <- read_tsv("bindata.40000.hg19.tsv.gz")

fire <- read_bed("fire-adult-hg19.txt")

# not necessary, but nice to have
si <- Seqinfo(genome="hg19")
si <- keepStandardChromosomes(si)

mikelove

# minimap2 map reads:
minimap2 -t {threads} -a -x map-ont -p 1.0 -N 100 {input.index} {input.fastq} | samtools view -Sb > {output.bam}
# salmon Quant:
salmon quant --ont --noLengthCorrection -p 8 -t {input.transcriptome} -l U -a {input.bam} -o {output}

mikelove

Index CHR,Index SNP,BP Left,BP Right,# Conditional independent,Index BP,Index P value,Merged Indexes,2nd SNP,2nd BP,2nd P value,3rd SNP,3rd BP,3rd P value,4th SNP,4th BP,4th P value,5th SNP,5th BP,5th P value,,,,,,,,
6,rs13195636,24939493,30559493,4,27509493,3.33E-40,rs13195636,rs9461856,33395199,4.13E-14,rs8192589,32187637,2.27E-10,rs9368789,34001678,2.03E-07,rs10946808,26233387,1.11E-06,,,,,,,,
18,rs9636107,52666708,53680258,4,53200117,1.92E-13,"rs12969453,rs74914300,rs72934602,rs9636107",rs12969453,52751708,3.93E-11,rs78322266,53063676,6.89E-08,rs72934602,53568458,6.79E-07,rs9952704,53447690,2.60E-04,,,,,,,,
3,rs2710323,48132866,53603354,3,52815905,1.02E-21,"rs11715134,rs11917680,rs2710323,rs13080668",rs2236989,50505395,1.81E-10,rs7633840,48719638,4.41E-07,rs4687724,53408177,1.10E-04,,,,,,,,,,,
18,rs72980087,77324421,77981194,2,77632194,4.80E-16,"rs56040937,rs7238071,rs72980087,rs72970145",rs11664298,77578986,5.24E-16,rs12455965,77688830,0.0001109,,,,,,,,,,,,,,
3,rs13090130,161092491,161969035,2,161777035,

mikelove

library(plyranges)
x <- data.frame(seqnames=1,
                start=sample(1000,10,FALSE),
                width=1, id=1:10) %>%
  as_granges()
y <- data.frame(seqnames=1,
                start=sample(1000,10,FALSE),
                width=1, id=letters[1:10]) %>%
  as_granges()

mikelove

library(plyranges)

set.seed(1)
x <- data.frame(seqnames=1, start=0:9 * 100 + 1, 
                width=20, id=1:10) %>%
  as_granges()

y <- data.frame(seqnames=1, start=round(runif(4,100,900)), 
                width=10, id=letters[1:4]) %>%
  as_granges() %>%

mikelove

---
title: "Using RNA-seq DE methods to detect allele-specific expression"
author: "Michael Love"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

The question has come up a few times on the Bioconductor support site how to use
methods like DESeq2 or edgeR to detect allele-specific expression, 
and how to see if the allelic ratio differs across condition. Aaron Lun has written up 

mikelove

airway@assays@data[["counts"]][airway@rowRanges@elementMetadata@listData$gene_id == "ENSG00000000003", airway@colData@listData$dex == "trt"]
assays(airway)[["counts"]][names(rowRanges(airway)) == "ENSG00000000003", colData(airway)$dex == "trt"]
assay(airway, "counts")["ENSG00000000003", airway$dex == "trt"]
airway |> filter(symbol == "TSPAN6", dex == "trt") |> pull(counts)

mikelove

; mike love's .emacs

; general stuff
(menu-bar-mode -1)
(tool-bar-mode -1)
(global-set-key "\C-x\C-b" 'electric-buffer-list)
(global-unset-key (kbd "\C-x DEL") )
(global-unset-key (kbd "\C-t") )
(setq inhibit-startup-screen t)
(setq backup-directory-alist '(("" . "~/emacs-backup")))

mikelove

# scraped from https://doi.org/10.1038/nature11241
# note that intersections < 100 not included
banana_sets <- c(
	Musa = 759,
	Phoenix = 769,
	Sorghum = 827,
	Brachypodium = 387,
	Ozyza = 1246,
	Arabidopsis = 1187,
	`Musa&Phoenix` = 467,