DE_mix.Rmd

---
title: "Pseudo mixture analysis for differential transcript expression analysis"
author: "Xueyi Dong"
date: "28/09/2021"
output: 
  html_document: 
    toc: yes
    toc_float: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(Rsubread)
library(limma)
library(edgeR)
library(ggplot2)
library(DESeq2)
library(EBSeq)
library(NOISeq)
library(wesanderson)
library(cowplot)
```

```{r}
# Set output folder
OUT="./DTEmix"
```


```{r}
# Load Salmon quant TPM data
samples <-list.files("salmon_bs")
quant <- file.path("./salmon_bs", samples, "quant.sf")
library(tximport)
txi <- tximport(quant, type="salmon", txOut=TRUE, countsFromAbundance = "no")
counts <- as.data.frame(txi$counts, stringAsFactors = FALSE)
colnames(counts) <- samples

# deal with gene names
genes <- rownames(counts)
genes <- strsplit2(genes, "|", fixed=TRUE)
rownames(counts) <- genes[,1]

library(GenomicFeatures)
gtf <- "/wehisan/home/allstaff/d/dong.x/annotation/HumanSequins/gencode.v33.sequins.gtf"
txdb <- makeTxDbFromGFF(gtf)
# saveRDS(txdb, "txdb.RDS")
# txdb <- readRDS("txdb.RDS")
txdf <- select(txdb, keys(txdb, "GENEID"), "TXNAME", "GENEID") 
tab <- table(txdf$GENEID) 
txdf$ntx <- tab[match(txdf$GENEID, names(tab))]
```


## DE analysis using limma

```{r}
x <- DGEList(counts = counts)
#HCC827 proportion
group <- rep(c("000", "100", "075", "050", "025"), rep(3, 5))
x$samples$group <- group
x$samples$mixprop <- rep(c(0, 1, 0.75, 0.5, 0.25), rep(3, 5))
x$samples$pair <- rep(c(1, 2, 3), 5)
x$genes <- txdf[match(rownames(x), txdf$TXNAME),]
```

```{r}
# geneid <- strsplit2(rownames(x$counts), "\\.")[,1]
# genes <- select(Homo.sapiens, keys=geneid, columns=c("SYMBOL", "TXCHROM"), 
#                 keytype="ENSEMBL")
# x$genes <- genes[match(geneid, genes$ENSEMBL),]
```

```{r}
filter <- filterByExpr(x)
x <- x[filter, ,keep.lib.sizes=FALSE]
dim(x)

library(RColorBrewer)
col <- brewer.pal(5, "Set1")

x <- calcNormFactors(x)
cpm <- cpm(x, log=TRUE)
plotMDS(cpm, labels = x$samples$group)
x$samples

human.gene <- grep("^ENST", rownames(x))
sequin.gene <- grep("^R", rownames(x))
cpm.human <- cpm(x[human.gene, ], log=TRUE)
cpm.sequin <- cpm(x[sequin.gene, ], log=TRUE)
# par(mfrow=c(1,2))
pdf(file.path(OUT, "mds_human.pdf"), height = 4, width = 8)
par(mar=c(5.1, 5.1, 3.1, 6.1), xpd=TRUE)
plotMDS(cpm.human, pch = rep(c(1, 2), c(6, 9)),
        col = rep(col, rep(3, 5)), main="ONT human",
        cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)
legend("topright", inset=c(-0.16,0), c("000", "025", "050", "075", "100"), 
       pch = c(1, 2, 2, 2, 1),
       col = col[c(1, 5, 4, 3, 2)],
       title = "Group", title.col = "black",
       text.col = col[c(1, 5, 4, 3, 2)], cex = 1.25)
dev.off()
pdf(file.path(OUT, "mds_human_noLeg.pdf"), height = 4, width = 3.5)
plotMDS(cpm.human, pch = rep(c(1, 2), c(6, 9)),
        col = rep(col, rep(3, 5)), main="ONT human",
        cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)
dev.off()
pdf(file.path(OUT, "mds_sequin.pdf"), height = 5, width = 8)
plotMDS(cpm.sequin, pch = rep(c(1, 2), c(6, 9)),
        col = rep(col, rep(3, 5)), main = "ONT sequin",
        cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.5, cex = 1.5)
legend("topleft", c("000", "025", "050", "075", "100"), 
        text.col = col[c(1, 5, 4, 3, 2)], cex = 1.25, bty = "n")
dev.off()
```

### Make design and contrast

```{r}
comparison <- list(
  c100vs0=c("100", "000"),
  c75vs25=c("075", "025"),
  c50vs25=c("050", "025"),
  c75vs50=c("075", "050")
)


anno <- read.table("/wehisan/home/allstaff/d/dong.x/annotation/sequins/rnasequin_isoforms_2.4.tsv", header = TRUE, stringsAsFactors = FALSE)
anno$logFC <- log(anno$MIX_B / anno$MIX_A)
```

### fit a separate model for sequins genes

```{r}
x.sequin <- x[sequin.gene,]
x.sequin <- calcNormFactors(x.sequin)
tt.sequin <- lapply(comparison, function(x){
  x.tmp <- x.sequin[, x.sequin$samples$group %in% x]
  design <- model.matrix(~x.tmp$samples$group)
  v.sequin <- voom(x.tmp, design = design, plot=TRUE)
  fit.sequin <- lmFit(v.sequin)
  # fit.sequin <- contrasts.fit(fit.sequin, contrasts=contr)
  efit.sequin <- eBayes(fit.sequin)
  # dt.sequin <- decideTests(efit.sequin)
  topTable(efit.sequin, coef=ncol(design), n=Inf)
})
for(i in 1:4){
  write.table(tt.sequin[[i]], file = paste(OUT, "/topTableSequin", names(tt.sequin)[i], ".tsv", sep=""),
              sep = "\t")
}
```
### Fit a separate model for human genes

```{r}
x.human <- x[human.gene,]
x.human <- calcNormFactors(x.human)
tt.human <- lapply(comparison, function(x){
  x.tmp <- x.human[, x.human$samples$group %in% x]
  design <- model.matrix(~x.tmp$samples$group)
  v.human <- voom(x.tmp, design = design, plot=TRUE)
  fit.human <- lmFit(v.human)
  # fit.human <- contrasts.fit(fit.human, contrasts=contr)
  efit.human <- eBayes(fit.human)
  # dt.human <- decideTests(efit.human)
  topTable(efit.human, coef=ncol(design), n=Inf)
})
for(i in 1:4){
  write.table(tt.human[[i]], file = paste(OUT, "/topTableHuman", names(tt.human)[i], ".tsv", sep=""),
              sep = "\t")
}
```

Get DE genes

```{r}
DE.human.limma <- lapply(tt.human, function(x){
  rownames(x)[x$adj.P.Val < 0.05]
})
DE.sequin.limma <- lapply(tt.sequin, function(x){
  rownames(x)[x$adj.P.Val < 0.05]
})
tx.human.limma <- lapply(tt.human, function(x){
  return(x$TXNAME)
})
tx.sequin.limma <- lapply(tt.sequin, function(x){
  return(x$TXNAME)
})
```


```{r}
# mix.fit.human <- fitmixture(v.human$E, mixprop = x$samples$mixprop)
```



## DE analysis using edgeR


### Human genes

```{r}
ttag.human <- lapply(comparison, function(x){
  x.tmp <- x.human[, x.human$samples$group %in% x]
  design <- model.matrix(~x.tmp$samples$group)
  x.tmp <- estimateDisp(x.tmp, design)
  # plotBCV(x.human)
  qlfit.human <- glmQLFit(x.tmp, design)
  # plotQLDisp(qlfit.human)
  # res.human <- lapply(colnames(contr), function(x){ glmQLFTest(qlfit.human, contrast=contr[,x])})
  res.human <- glmQLFTest(qlfit.human)
  as.data.frame(topTags(res.human, n=Inf))
  # ttag.human <- ?lapply(res.human, function(x){topTags(x, n=Inf)})
})

# write.table(ttag.human[,c("GeneID","logFC", "FDR")], file.path(OUT, "resEdgerHuman.tsv"), sep="\t", row.names = TRUE)

DE.human.edgeR <- lapply(ttag.human, function(x){
  x$TXNAME[x$FDR < 0.05]
})
tx.human.edgeR <- lapply(ttag.human, function(x){
  return(x$TXNAME)
})
```


### Sequin genes

```{r}
ttag.sequin <- lapply(comparison, function(x){
  x.tmp <- x.sequin[, x.sequin$samples$group %in% x]
  design <- model.matrix(~x.tmp$samples$group)
  x.tmp <- estimateDisp(x.tmp, design)
  # plotBCV(x.sequin)
  qlfit.sequin <- glmQLFit(x.tmp, design)
  # plotQLDisp(qlfit.sequin)
  # res.sequin <- lapply(colnames(contr), function(x){ glmQLFTest(qlfit.sequin, contrast=contr[,x])})
  res.sequin <- glmQLFTest(qlfit.sequin)
  as.data.frame(topTags(res.sequin, n=Inf))
  # ttag.sequin <- ?lapply(res.sequin, function(x){topTags(x, n=Inf)})
})

# write.table(ttag.sequin[,c("GeneID","logFC", "FDR")], file.path(OUT, "resEdgerHuman.tsv"), sep="\t", row.names = TRUE)

DE.sequin.edgeR <- lapply(ttag.sequin, function(x){
  x$TXNAME[x$FDR < 0.05]
})
tx.sequin.edgeR <- lapply(ttag.sequin, function(x){
  return(x$TXNAME)
})
```

## DE analysis using DESeq2
```{r}
human.gene <- grep("^ENST", rownames(counts))
sequin.gene <- grep("^R", rownames(counts))

deres.human <- lapply(comparison, function(x){
  sel = group %in% x
  groups <- group[sel]
  dds <- DESeqDataSetFromMatrix(countData=  round(counts[human.gene, sel]), 
                              colData = data.frame(
                                group = groups
                              ),
                              design = ~group)
  keep <- rowSums(counts(dds)) > 1 & rowSums(counts(dds) >= 10) >= 3
  dds <- dds[keep,]
  dds <- DESeq(dds)
  results(dds)
  })

deres.sequin <- lapply(comparison, function(x){
  sel = group %in% x
  groups <- group[sel]
  dds <- DESeqDataSetFromMatrix(countData=  round(counts[sequin.gene, sel]), 
                              colData = data.frame(
                                group = groups
                              ),
                              design = ~group)
  keep <- rowSums(counts(dds)) > 1 & rowSums(counts(dds) >= 10) >= 3
  dds <- dds[keep,]
  dds <- DESeq(dds)
  results(dds)
  })


DE.human.DESeq2 <- lapply(deres.human, function(x){
  x <- na.omit(x)
  rownames(x)[x$padj < 0.05]
})
DE.sequin.DESeq2 <- lapply(deres.sequin, function(x){
  x <- na.omit(x)
  rownames(x)[x$padj < 0.05]
})
tx.human.DESeq2 <- lapply(deres.human, function(x){
  return(rownames(x))
})
tx.sequin.DESeq2 <- lapply(deres.sequin, function(x){
  return(rownames(x))
})
```

## DE analysis using EBSeq

```{r}
human.gene <- grep("^ENST", rownames(counts))
sequin.gene <- grep("^R", rownames(counts))
counts.human <- as.matrix(counts[human.gene, ])
counts.sequin <- as.matrix(counts[sequin.gene, ])
Sizes.human = MedianNorm(counts.human)
Sizes.sequin = MedianNorm(counts.sequin)
Conditions = as.factor(rep(c("000", "100", "075", "050", "025"), rep(3, 5)))
comparisons = list(c("000", "100"), 
                          c("025", "075"),
                          c("025", "050"),
                          c("050", "075"))
EBOut.human = lapply(comparisons,
                     function(x){
                       sel = which(as.character(Conditions) %in% x)
                       out = EBTest(Data=counts.human[,sel], 
                                    Conditions = factor(Conditions[sel]),
                                    sizeFactors = Sizes.human[sel],
                                    maxround = 5)
                       out})
EBDERes.human = lapply(EBOut.human, function(x){GetDEResults(x, FDR = 0.05)})
EBOut.sequin = lapply(comparisons,
                     function(x){
                       sel = which(as.character(Conditions) %in% x)
                       out = EBTest(Data=counts.sequin[,sel], 
                                    Conditions = factor(Conditions[sel]),
                                    sizeFactors = Sizes.sequin[sel],
                                    maxround = 5)
                       out})
EBDERes.sequin = lapply(EBOut.sequin, function(x){GetDEResults(x, FDR = 0.05)})

GeneFC.human = lapply(EBOut.human, PostFC)
GeneFC.sequin = lapply(EBOut.sequin, PostFC)

DE.human.EBSeq <- lapply(EBDERes.human, function(x){x$DEfound})
DE.sequin.EBSeq <- lapply(EBDERes.sequin, function(x){x$DEfound})
tx.human.EBSeq <- lapply(EBDERes.human, function(x){
  names(x$Status)[x$Status %in% c("DE", "EE")]
})
tx.sequin.EBSeq <- lapply(EBDERes.sequin, function(x){
  names(x$Status)[x$Status %in% c("DE", "EE")]
})
```


## DE analysis using NOISeq

```{r}
myfilt.human = filtered.data(counts.human, factor = x$samples$group, norm=FALSE, method = 1)
myfilt.sequin = filtered.data(counts.sequin, factor = x$samples$group, norm=FALSE, method = 1)
mydata.human <- readData(
  data = myfilt.human,
  factors = data.frame(
    group = rep(c("000", "100", "075", "050", "025"), rep(3, 5))
  )
)
mydata.sequin <- readData(
  data = myfilt.sequin,
  # length = mylength[sequin.gene,, drop=FALSE],
  factors = data.frame(
    group = rep(c("000", "100", "075", "050", "025"), rep(3, 5))
  )
)

noires.human <- lapply(comparison, function(x){
  mydata.human = mydata.human[, mydata.human$group %in% x]
  mynoiseqbio.human = noiseqbio(mydata.human,  norm="tmm", factor="group", )
  noires.human <- (mynoiseqbio.human@results[[1]])
  noires.human$FDR=1-noires.human$prob
  noires.human
})

noires.sequin <- lapply(comparison, function(x){
  mydata.sequin = mydata.sequin[, mydata.sequin$group %in% x]
  mynoiseqbio.sequin = noiseqbio(mydata.sequin,  norm="tmm", factor="group", )
  noires.sequin <- (mynoiseqbio.sequin@results[[1]])
  noires.sequin$FDR=1-noires.sequin$prob
  noires.sequin
})

# # write.table(noires.human[,c("gene", "log2FC", "FDR")], file.path(OUT, "resNoiseqHuman.tsv"), sep="\t", row.names = TRUE)
# write.table(noires.sequin[,c("gene", "log2FC", "FDR")], file.path(OUT, "resNoiseqSequin.tsv"), sep="\t", row.names = TRUE)
DE.human.NOISeq <- lapply(noires.human, function(x){
  rownames(x)[x$FDR < 0.05]
})
DE.sequin.NOISeq <- lapply(noires.sequin, function(x){
  rownames(x)[x$FDR < 0.05]
})
tx.human.NOISeq <- lapply(noires.human, function(x){
  return(rownames(x))
})
tx.sequin.NOISeq <- lapply(noires.sequin, function(x){
  return(rownames(x))
})
```

## Results

```{r}
DE.human <- list(DE.human.limma, DE.human.edgeR, DE.human.DESeq2, DE.human.EBSeq, DE.human.NOISeq)
DE.sequin <- list(DE.sequin.limma, DE.sequin.edgeR, DE.sequin.DESeq2, DE.sequin.EBSeq, DE.sequin.NOISeq)
saveRDS(DE.human, file.path(OUT, "DE.human.RDS"))
saveRDS(DE.sequin, file.path(OUT, "DE.sequin.RDS"))
# see which tx were used in analysis
tx.human <- list(tx.human.limma, tx.human.edgeR, tx.human.DESeq2, tx.human.EBSeq, tx.human.NOISeq)
tx.sequin <- list(tx.sequin.limma, tx.sequin.edgeR, tx.sequin.DESeq2, tx.sequin.EBSeq, tx.sequin.NOISeq)
saveRDS(tx.human, file.path(OUT, "tx.DE.human.RDS"))
saveRDS(tx.sequin, file.path(OUT, "tx.DE.sequin.RDS"))
```


### Calculate recovery rate

```{r}
calcRecovery <- function(DE){
  tmp <- sapply(2:length(DE), function(x){
    sum(DE[[1]] %in% DE[[x]])
  }, simplify=TRUE)
  c(length(DE[[1]]), tmp)/length(DE[[1]])
}

rec.human <- as.data.frame(sapply(DE.human,
                    calcRecovery,
                    simplify = TRUE))
rownames(rec.human) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(rec.human) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
rec.human <- melt(rec.human, variable.name = "Method", value.name = "Recovery rate")
rec.human$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "humanRecovery.pdf"), height = 5, width = 8)
ggplot(rec.human, aes(x=Comparison, y=`Recovery rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "longdash", "dashed", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
dev.off()

rec.sequin <- as.data.frame(sapply(DE.sequin,
                    calcRecovery,
                    simplify = TRUE))
rownames(rec.sequin) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(rec.sequin) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
rec.sequin <- melt(rec.sequin, variable.name = "Method", value.name = "Recovery rate")
rec.sequin$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "sequinRecovery.pdf"), height = 5, width = 8)
ggplot(rec.sequin, aes(x=Comparison, y=`Recovery rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "solid", "dashed", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw()  +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
dev.off()
```

### Calculate inconsistancy rate

```{r}
calcInconsistancy <- function(DE){
  tmp <- sapply(1:length(DE), function(x){
    sum(!(DE[[x]] %in% DE[[1]]))
  }, simplify=TRUE)
  tmp/length(DE[[1]])
}

inc.human <- as.data.frame(sapply(DE.human,
                    calcInconsistancy,
                    simplify = TRUE))
rownames(inc.human) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(inc.human) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
inc.human <- melt(inc.human, variable.name = "Method", value.name = "Inconsistancy rate")
inc.human$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "humanInconsistancy.pdf"), height = 5, width = 8)
humanInc <- 
  ggplot(inc.human, aes(x=Comparison, y=`Inconsistancy rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "longdash", "dashed", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
print(humanInc)
dev.off()

inc.sequin <- as.data.frame(sapply(DE.sequin,
                    calcInconsistancy,
                    simplify = TRUE))
rownames(inc.sequin) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(inc.sequin) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
inc.sequin <- melt(inc.sequin, variable.name = "Method", value.name = "Inconsistancy rate")
inc.sequin$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "sequinInconsistancy.pdf"), height = 5, width = 8)
ggplot(inc.sequin, aes(x=Comparison, y=`Inconsistancy rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "solid", "longdash", "dashed", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
dev.off()
```

### Power (number of DE genes)

```{r}
power.human <- as.data.frame(sapply(DE.human,
                    function(x){sapply(x, length, simplify = TRUE)},
                    simplify = TRUE))
rownames(power.human) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(power.human) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
power.human <- melt(power.human, variable.name = "Method", value.name = "Number of DE transcripts")
power.human$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "humanPower.pdf"), height = 5, width = 8)
humanPow <- 
  ggplot(power.human, aes(x=Comparison, y=`Number of DE transcripts`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "longdash", "dashed", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
print(humanPow)
dev.off()
```

```{r}
power.sequin <- as.data.frame(sapply(DE.sequin,
                    function(x){sapply(x, length, simplify = TRUE)},
                    simplify = TRUE))
rownames(power.sequin) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(power.sequin) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
library(data.table)
power.sequin <- melt(power.sequin, variable.name = "Method", value.name = "Number of DE transcripts")
power.sequin$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "sequinPower.pdf"), height = 5, width = 8)
ggplot(power.sequin, aes(x=Comparison, y=`Number of DE transcripts`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "longdash", "dashed", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
dev.off()
```

### Sequin FDR

```{r}
trueDE <- anno$NAME[anno$logFC != 0]
calcFDR <- function(DE){
  sapply(DE, function(x){
    FD = sum(!(x %in% trueDE))
    FD / length(x)
  }, simplify=TRUE)
}
                    
FDR.sequin <- as.data.frame(sapply(DE.sequin,
                     calcFDR, 
                     simplify = TRUE))

rownames(FDR.sequin) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(FDR.sequin) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
FDR.sequin <- melt(FDR.sequin, variable.name = "Method", value.name = "False discovery rate")
FDR.sequin$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "sequinFDR.pdf"), height = 5, width = 8)
sequinFDR <- 
  ggplot(FDR.sequin, aes(x=Comparison, y=`False discovery rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "dashed", "solid", "longdash", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
print(sequinFDR)
dev.off()
```

### Sequin TPR

```{r}
calcTPR <- function(DE){
  sapply(DE, function(x){
    TP = sum(x %in% trueDE)
    TP / length(trueDE)
  }, simplify=TRUE)
}
                    
TPR.sequin <- as.data.frame(sapply(DE.sequin,
                     calcTPR, 
                     simplify = TRUE))

rownames(TPR.sequin) <- c("100vs000", "075vs025", "050vs025", "075vs050")
colnames(TPR.sequin) <- c("limma", "edgeR", "DESeq2", "EBSeq", "NOISeq")
TPR.sequin <- melt(TPR.sequin, variable.name = "Method", value.name = "True positive rate")
TPR.sequin$Comparison <- factor(rep(c("100vs000", "075vs025", "050vs025", "075vs050"), 5), levels = c("100vs000", "075vs025", "050vs025", "075vs050"))

pdf(file.path(OUT, "sequinTPR.pdf"), height = 5, width = 8)
sequinTPR <- 
  ggplot(TPR.sequin, aes(x=Comparison, y=`True positive rate`, colour=Method, group=Method)) +
  geom_line(linetype=rep(c("solid", "dashed", "longdash", "solid", "solid"), rep(4, 5))) +
  geom_point(size=3, shape=rep(c(1, 2, 3, 4, 5), rep(4, 5))) +
  theme_bw() +
  scale_colour_manual(values = c("#D96A70", "#476937",  "#D5A2CB", "#708FA6", "#9FC675"))+
  theme(text = element_text(size = 20))
print(sequinTPR)
dev.off()
```
### Compose figure to be used in main figure

```{r}
library(cowplot)
leg <- get_legend(humanPow + theme(legend.position = "bottom", text = element_text(size = 20)))
pdf(file.path(OUT, "analysisMix.pdf"), height = 8, width = 16)
plot_grid(
  humanPow + theme(legend.position = "none"),
  humanInc + theme(legend.position = "none"),
  sequinTPR + theme(legend.position = "none"),
  sequinFDR + theme(legend.position = "none"),
  leg,
  ncol = 2,
  rel_heights = c(1, 1, 0.1),
  labels = c("d", "e", "f", "g", ""),
  label_size = 24
)
dev.off()
```


### Sequin logFC expected vs estimated

Only calculate for pure samples comparisons

```{r}
GeneFC.sequin = PostFC(EBOut.sequin[[1]])
anno$logFC_limma <- tt.sequin[[1]]$logFC[match(anno$NAME, tt.sequin[[1]]$TXNAME)]
anno$logFC_edgeR <- ttag.sequin[[1]]$logFC[match(anno$NAME, ttag.sequin[[1]]$TXNAME)]
anno$logFC_DESeq2 <- deres.sequin[[1]]$log2FoldChange[match(anno$NAME, rownames(deres.sequin[[1]]))]
anno$logFC_EBSeq <- -log2(GeneFC.sequin$PostFC)[match(anno$NAME, names(GeneFC.sequin$PostFC))]
anno$logFC_NOISeq <- -noires.sequin[[1]]$log2FC[match(anno$NAME, rownames(noires.sequin[[1]]))]
anno$DE_limma <- anno$NAME %in% DE.sequin.limma[[1]]
anno$DE_edgeR <- anno$NAME %in% DE.sequin.edgeR[[1]]
anno$DE_DESeq2 <- anno$NAME %in% DE.sequin.DESeq2[[1]]
anno$DE_EBSeq <- anno$NAME %in% DE.sequin.EBSeq[[1]]
anno$DE_NOISeq <- anno$NAME %in% DE.sequin.NOISeq[[1]]
```

```{r}
lfclm <- lm(anno$logFC~anno$logFC_limma)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2,
                 list(r2 = format(summary(lfclm)$r.squared, digits = 3)))))
# pdf(file.path(OUT, "sequinlogfclimma.pdf"), height = 4)
p1 <- ggplot(anno, aes(x=logFC, y=logFC_limma))+
  geom_point(aes(colour=DE_limma)) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  annotate(geom="text", x = -4, y = 8, label = eq, parse = TRUE, size = 6) +
  # geom_abline(intercept = 0, slope = 1) +
  theme_bw()+
  theme(text = element_text(size = 20), legend.position = "none") +
  ggtitle("limma")
  
# dev.off()

lfclm <- lm(anno$logFC~anno$logFC_edgeR)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2,
                 list(r2 = format(summary(lfclm)$r.squared, digits = 3)))))
# pdf(file.path(OUT, "sequinlogfcedgeR.pdf"), height = 4)
p2 <- ggplot(anno, aes(x=logFC, y=logFC_edgeR))+
  geom_point(aes(colour = DE_edgeR)) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  annotate(geom="text", x = -4, y = 8, label = eq, parse = TRUE, size = 6) +
  # geom_abline(intercept = 0, slope = 1) +
  theme_bw()+
  theme(text = element_text(size = 20), legend.position = "none") +
  ggtitle("edgeR")
# dev.off()

lfclm <- lm(anno$logFC~anno$logFC_DESeq2)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2,
                 list(r2 = format(summary(lfclm)$r.squared, digits = 3)))))
# pdf(file.path(OUT, "sequinlogfcDESeq2.pdf"), height = 4)
p3 <- ggplot(anno, aes(x=logFC, y=logFC_DESeq2))+
  geom_point(aes(colour=DE_DESeq2)) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  annotate(geom="text", x = -4, y = 8, label = eq, parse = TRUE, size = 6) +
  # geom_abline(intercept = 0, slope = 1) +
  theme_bw()+
  theme(text = element_text(size = 20), legend.position = "none") +
  ggtitle("DESeq2")
# dev.off()

lfclm <- lm(anno$logFC~anno$logFC_EBSeq)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2,
                 list(r2 = format(summary(lfclm)$r.squared, digits = 3)))))
# pdf(file.path(OUT, "sequinlogfcEBSeq.pdf"), height = 4)
p4 <- ggplot(anno, aes(x=logFC, y=logFC_EBSeq))+
  geom_point(aes(colour=DE_EBSeq)) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x)+
  annotate(geom="text", x = -4, y = 8, label = eq, parse = TRUE, size = 6) +
  # geom_abline(intercept = 0, slope = 1) +
  theme_bw()+
  theme(text = element_text(size = 20), legend.position = "none") +
  ggtitle("EBSeq")
# dev.off()

lfclm <- lm(anno$logFC~anno$logFC_NOISeq)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2,
                 list(r2 = format(summary(lfclm)$r.squared, digits = 3)))))
# pdf(file.path(OUT, "sequinlogfcNOISeq.pdf"), height = 4)
p5 <- ggplot(anno, aes(x=logFC, y=logFC_NOISeq))+
  geom_point(aes(colour=DE_NOISeq)) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  annotate(geom="text", x = -4, y = 8, label = eq, parse = TRUE, size = 6) +
  # geom_abline(intercept = 0, slope = 1) +
  theme_bw()+
  theme(text = element_text(size = 20)) +
  ggtitle("NOISeq") +
  labs(colour="is DTE")
# dev.off()

leg <- get_legend(p5)
p5 <- p5 +
  theme(legend.position = "none")
pdf(file.path(OUT, "sequinlogfc.pdf"), height = 10, width = 15)
plot_grid(p1, p2, p3, p4, p5, leg, ncol=3)
dev.off()
```

### sequin count vs annotation
```{r}
m <- match(anno$NAME, rownames(cpm.sequin))
annoData <- cbind(anno[,1:4], cpm.sequin[m,1:6])
annoData <- data.table::melt(annoData, id.vars = 1:4)
annoData$exp <-annoData$MIX_A
annoData$exp[annoData$variable %in% c("barcode04", "barcode05", "barcode06")] <- annoData$MIX_B[annoData$variable %in% c("barcode04", "barcode05", "barcode06")]
colnames(annoData)[5] <- "sample"
cpmablm <- lm(log(annoData$exp)~annoData$value)
conf <- confintr::ci_rsquared(cpmablm)
eq <- as.character(as.expression(substitute(italic(r)^2~"="~r2~"("*int1*"-"*int2*")",
                 list(r2 = format(summary(cpmablm)$r.squared, digits = 3),
                      int1 = format(conf$interval[1], digits = 3),
                      int2 = format(conf$interval[2], digits = 3)))))
library(scales)
library(MetBrewer)
pdf(file.path(OUT, "sequinCPMvsAbundance.pdf"), height = 5, width = 8)
p <- ggplot(na.omit(annoData), aes(x=exp, y=value, colour=sample)) +
  geom_point() +
  theme_bw() +
  scale_x_continuous(trans = "log10", labels = comma) +
  # scale_y_continuous(trans = "log10", labels=comma, limits=c(0.5, 180000)) +
  labs(x = "Expected abundance", y = "log CPM") +
  scale_colour_manual(values = met.brewer("Troy", 6), 
                      name = "Sample", 
                      labels = paste(
                        rep(c("H1975", "HCC827"), c(3,3)),
                        c(1, 2, 3, 1, 2, 5),
                        sep = "-"), )+
  theme(text = element_text(size = 20)) +
  annotate(geom="text", x = 1, y = 11, label = eq, parse = TRUE, size = 6) 
  # annotate(geom="text", x=300, y=160000, label=paste0("Pearson's r=", round(cor(annoData$exp, annoData$value, use = "complete.obs"), 2)), size=8) 
plot(p)
dev.off()
saveRDS(p, "sequinCPMvsAbundancePlot.RDS")

```

Dig a little bit further into the false positives / not well estimated logFCs

```{r}
anno$AverageAbundance <- (anno$MIX_A + anno$MIX_B) / 2
ggplot(anno, aes(x=logFC, y=logFC_limma))+
  geom_point(aes(colour=AverageAbundance, shape=DE_limma), size=3) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  theme_bw()+
  theme(text = element_text(size = 20)) +
  scale_colour_viridis_c(trans="log")
ggplot(anno, aes(x=logFC, y=logFC_limma))+
  geom_point(aes(colour=LENGTH, shape=DE_limma), size=3) +
  geom_smooth(alpha=0.5, method='lm', formula= y~x) +
  theme_bw()+
  theme(text = element_text(size = 20)) +
  scale_colour_viridis_c(trans="log")

summary(abs(anno$logFC_limma-anno$logFC))
na.omit(anno[abs(anno$logFC_limma-anno$logFC) > 1,])
```

## Session information

```{r}
save.image("DE_mix.RData")
sessionInfo()
```