DESeq2_MACS_Adult_Bulk

load("../1.align/txi.rsem.RData")
metadata <- read.csv("metadata.csv", stringsAsFactors = TRUE)
metadata$Trt <- factor(metadata$Trt, levels=c("SAL","LPS"))
metadata$condition1 <- as.factor(paste(metadata$Group, metadata$Region, metadata$Trt, sep = "."))
metadata$condition2 <- factor(paste(metadata$Group, metadata$Trt, sep = "."), levels=c("CON.SAL", "CON.LPS", "MIA.SAL", "MIA.LPS"))
rownames(metadata) <- metadata$ID
all.equal(colnames(txi.rsem$counts), as.character(metadata$ID))

## [1] TRUE

referece: https://rpubs.com/ge600/deseq2

Differential gene expression for design: Group + Trt + Region + Group:Trt

dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Group + Trt + Region + Group:Trt)
dds <- dds[,-21] # exclude MSL3 outlier by cluster

# filter low expressed genes
keep <- rowSums(counts(dds)) >= 50
table(keep)

## keep
## FALSE  TRUE 
## 37491 17996

dds <- dds[keep,]

# add in gene symbols
genes <- gsub("\\..*","",rownames(dds))
# symbols <- mapIds(org.Mm.eg.db, keys=genes, column="SYMBOL", keytype="ENSEMBL", multiVals="first")
library(AnnotationHub)
ah <- AnnotationHub()
edb <- ah[["AH53222"]]
symbol <- mapIds(edb, keys=genes, column="SYMBOL", keytype="GENEID", multiVals="first")
mcols(dds) <- cbind(mcols(dds), symbol)

dds <- DESeq(dds)
resultsNames(dds)

## [1] "Intercept"         "Group_MIA_vs_CON"  "Trt_LPS_vs_SAL"   
## [4] "Region_Str_vs_Ctx" "GroupMIA.TrtLPS"

# The effect of treatment in CON (the main effect)
summary(results(dds, contrast=c("Trt","LPS","SAL"), alpha=0.05))

## 
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2970, 17%
## LFC < 0 (down)     : 2654, 15%
## outliers [1]       : 51, 0.28%
## low counts [2]     : 339, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# up: 2970 down: 2654

# The effect of treatment in MIA (the main effect)
summary(results(dds, list( c("Trt_LPS_vs_SAL","GroupMIA.TrtLPS")), alpha=0.05))

## 
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 2288, 13%
## LFC < 0 (down)     : 1996, 11%
## outliers [1]       : 51, 0.28%
## low counts [2]     : 684, 3.8%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# up: 2288 down: 1996

# difference between MIA and CON in SAL
summary(results(dds, name="Group_MIA_vs_CON", alpha=0.05))

## 
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 5, 0.028%
## LFC < 0 (down)     : 2, 0.011%
## outliers [1]       : 51, 0.28%
## low counts [2]     : 1030, 5.7%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# up: 5 down: 2

# difference between MIA and CON in LPS
summary(results(dds, list( c("Group_MIA_vs_CON","GroupMIA.TrtLPS")), alpha=0.05))

## 
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 96, 0.53%
## LFC < 0 (down)     : 305, 1.7%
## outliers [1]       : 51, 0.28%
## low counts [2]     : 4502, 25%
## (mean count < 17)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# up: 96 down: 305

# different Trt response between 2 groups
summary(results(dds, name="GroupMIA.TrtLPS", alpha=0.05))

## 
## out of 17996 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 21, 0.12%
## LFC < 0 (down)     : 76, 0.42%
## outliers [1]       : 51, 0.28%
## low counts [2]     : 2768, 15%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# up: 21 down: 76

LvSinCON <- results(dds, contrast=c("Trt","LPS","SAL"), alpha=0.05)
LvSinMIA <- results(dds, list( c("Trt_LPS_vs_SAL","GroupMIA.TrtLPS")), alpha=0.05)
MvCinSAL <- results(dds, name="Group_MIA_vs_CON", alpha=0.05)
MvCinLPS <- results(dds, list( c("Group_MIA_vs_CON","GroupMIA.TrtLPS")), alpha=0.05)
MLPSvCLPS <- results(dds, name="GroupMIA.TrtLPS", alpha=0.05)

LvSinCON$gene <- symbol
LvSinMIA$gene <- symbol
MvCinSAL$gene <- symbol
MvCinLPS$gene <- symbol
MLPSvCLPS$gene <- symbol

Heatmaps (ED2)

vsd <- vst(dds, blind=FALSE)
rownames(vsd) <- symbol
m<-plotPCA(vsd, intgroup=c("Group", "Trt"), returnData = TRUE)
ggplot(m, aes(PC1, PC2, color = group)) + geom_point(size = 3) + scale_color_manual(values = c("black", "green3", "red", "magenta")) + theme_cowplot() + xlab("PC1: 69% variance") + ylab("PC2: 8% variance")

library("pheatmap")
df <- as.data.frame(colData(dds)[,c("Group","Region","Trt")])
ann_colors = list(
    Group = c(CON = "black",MIA = "red"),
    Trt = c(SAL = "blue", LPS = "green"),
    Region = c(Ctx = "darkviolet", Str = "deeppink"))

DEG <- order(MvCinSAL$pvalue)[1:20]
pheatmap(assay(vsd)[DEG,], cluster_rows=TRUE, cluster_cols=TRUE, annotation_col=df, scale="row", border_color = NA, col = colorRampPalette(c("navy", "white", "firebrick3"))(50), annotation_colors = ann_colors, main="MIA vs CON in SAL")

DEG <- order(MvCinLPS$pvalue)[1:400]
pheatmap(assay(vsd)[DEG,], cluster_rows=TRUE, show_rownames=FALSE, cluster_cols=TRUE, annotation_col=df, scale="row", col = colorRampPalette(c("navy", "white", "firebrick3"))(50), annotation_colors = ann_colors, main="MIA vs CON in LPS")

DEG <- order(MLPSvCLPS$pvalue)[1:100]
pheatmap(assay(vsd)[DEG,], cluster_rows=TRUE, show_rownames=TRUE, cluster_cols=TRUE, annotation_col=df, border_color = NA, scale="row", col = colorRampPalette(c("navy", "white", "firebrick3"))(50), annotation_colors = ann_colors, main="LPS Response Differnce between MIA & CON")

# scaled transformed counts

LPS Response in CON vs MIA

LPS.Response <- merge(as.data.frame(LvSinCON[,c(7,1,2,3,4,5,6)]),as.data.frame(LvSinMIA[,c(2,3,4,5,6)]), by="row.names", all.x=TRUE)
rownames(LPS.Response) <- LPS.Response$Row.names
colnames(LPS.Response) <- c("gene.id","gene.name","baseMean","log2FC.CON","lfcSE.CON","stat.CON","pval.CON","padj.CON","log2FC.MIA","lfcSE.MIA","stat.MIA","pval.MIA","padj.MIA")

#Assign the Group variable for coloring 4-way plot 
for (i in c(1:nrow(LPS.Response))){
  if (LPS.Response$padj.CON[i]<0.05 & (!is.na(LPS.Response$padj.CON[i]))){
    LPS.Response$CONSig[i]='YES'
  }
  else {
    LPS.Response$CONSig[i]='NO'
  }
}

for (i in c(1:nrow(LPS.Response))){
  if (LPS.Response$padj.MIA[i]<0.05 & (!is.na(LPS.Response$padj.MIA[i]))){
    LPS.Response$MIASig[i]='YES'
  }
  else {
    LPS.Response$MIASig[i]='NO'
  }
}

for (i in c(1:nrow(LPS.Response))){
  if (LPS.Response$CONSig[i]=='YES' & LPS.Response$MIASig[i]=='NO'){
    LPS.Response$Group[i]=1 }
  else if (LPS.Response$CONSig[i]=='NO' & LPS.Response$MIASig[i]=='YES'){
    LPS.Response$Group[i]=2}
  else if (LPS.Response$CONSig[i]=='YES' & LPS.Response$MIASig[i]=='YES'){
    LPS.Response$Group[i]=3}
  else {
    LPS.Response$Group[i]=4}
  
}
table(LPS.Response$Group)

## 
##     1     2     3     4 
##  2267   927  3357 11445

LPS.Response.final<-subset(LPS.Response, LPS.Response$Group != 4)

ggplotRegression <- function (fit) {
require(ggplot2)
ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
  geom_point() +
  stat_smooth(method = "lm", col = "blue", size=2, se = TRUE) +
  labs(subtitle = paste(" Slope =",signif(fit$coef[[2]], 3),
                        " P =",signif(summary(fit)$coef[2,4], 3),
                        "Adj R2 = ",signif(summary(fit)$adj.r.squared, 3),
                        "Intercept =",signif(fit$coef[[1]],3 )
                     ))
}

ggplot(data = LPS.Response.final, aes(x = log2FC.CON, y = log2FC.MIA, color=factor(Group))) + geom_point() + scale_color_manual(values = c("black", "red","hotpink"), labels = c("CON only","MIA only","CON & MIA")) + geom_abline(slope=1, size=0.5) + labs(x="Con logFC",y="MIA logFC", color="") + xlim(-10,11) + ylim(-10,10) + theme_cowplot() + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON>0 & log2FC.MIA >0 & log2FC.MIA > log2FC.CON)) # 916

## [1] 916  16

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON>0 & log2FC.MIA >0 & log2FC.MIA < log2FC.CON)) # 2260

## [1] 2260   16

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON<0 & log2FC.MIA <0 & log2FC.MIA > log2FC.CON)) # 1946

## [1] 1946   16

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON<0 & log2FC.MIA <0 & log2FC.MIA < log2FC.CON)) # 1050

## [1] 1050   16

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON>0 & log2FC.MIA <0)) # 219

## [1] 219  16

dim(LPS.Response.final %>% dplyr::filter(log2FC.CON<0 & log2FC.MIA >0)) # 160

## [1] 160  16

ggplotRegression(lm(LPS.Response.final$log2FC.MIA ~ LPS.Response.final$log2FC.CON, data=LPS.Response.final)) + geom_abline(slope=1) + labs(x="CON t-value (LPS/SAL)", y="MIA t-value (LPS/SAL)") + xlim(-10, 11) + ylim(-10, 11) + theme_cowplot()

## `geom_smooth()` using formula 'y ~ x'

Concordance Correlation Coefficient (CCC) of MIA and CON LPS Response

library(epiR)
# concordance correlation coefficient
epi.ccc(LPS.Response.final$log2FC.CON, LPS.Response.final$log2FC.MIA, ci = "z-transform", conf.level = 0.95)$rho.c

##     est lower upper
## 1 0.811 0.803 0.819