Evalute the blunted immune response in MIA in the frontal cortex compared to the striatum
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
Differential gene expression for design: Group + Trt + Group:Trt
dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Group + Trt + Group:Trt)
samp.Ctx <- colData(dds)$Region == "Ctx"
samp.Str <- colData(dds)$Region == "Str"
Ctx <- dds[,samp.Ctx]
Str <- dds[,samp.Str]
# filter low expressed genes
keep.ctx <- rowSums(counts(Ctx)) >= 50
table(keep.ctx)
## keep.ctx
## FALSE TRUE
## 39214 16273
Ctx <- Ctx[keep.ctx,]
keep.str <- rowSums(counts(Str)) >= 50
table(keep.str)
## keep.str
## FALSE TRUE
## 38878 16609
Str <- Str[keep.str,]
Ctx <- DESeq(Ctx)
resultsNames(Ctx)
## [1] "Intercept" "Group_MIA_vs_CON" "Trt_LPS_vs_SAL" "GroupMIA.TrtLPS"
Str <- DESeq(Str)
resultsNames(Str)
## [1] "Intercept" "Group_MIA_vs_CON" "Trt_LPS_vs_SAL" "GroupMIA.TrtLPS"
Ctx LPS Response in CON vs MIA
# add in gene symbols
library(AnnotationHub)
ah <- AnnotationHub()
edb <- ah[["AH53222"]]
genes <- gsub("\\..*","",rownames(Ctx))
symbol <- mapIds(edb, keys=genes, column="SYMBOL", keytype="GENEID", multiVals="first")
LvSinCON <- results(Ctx, contrast=c("Trt","LPS","SAL"), alpha=0.05)
LvSinMIA <- results(Ctx, list( c("Trt_LPS_vs_SAL","GroupMIA.TrtLPS")), alpha=0.05)
LvSinCON$gene <- symbol
LvSinMIA$gene <- symbol
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
## 1426 530 1900 12417
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(color = "grey50") +
stat_smooth(method = "lm", col = "blue", size=2) +
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 )
))
}
ggplotRegression(lm(LPS.Response.final$log2FC.MIA ~ LPS.Response.final$log2FC.CON, data=LPS.Response.final)) + geom_abline(slope=1) + labs(x="CON log2FC (LPS/SAL)", y="MIA log2FC (LPS/SAL)") + xlim(-10, 11) + ylim(-10, 11) + theme_cowplot()
Concordance Correlation Coefficient (CCC) of MIA and CON LPS Response in Ctx
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.836 0.827 0.845
STR LPS Response in CON vs MIA
# add in gene symbols
ah <- AnnotationHub()
edb <- ah[["AH53222"]]
genes <- gsub("\\..*","",rownames(Str))
symbol <- mapIds(edb, keys=genes, column="SYMBOL", keytype="GENEID", multiVals="first")
LvSinCON <- results(Str, contrast=c("Trt","LPS","SAL"), alpha=0.05)
LvSinMIA <- results(Str, list( c("Trt_LPS_vs_SAL","GroupMIA.TrtLPS")), alpha=0.05)
LvSinCON$gene <- symbol
LvSinMIA$gene <- symbol
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
## 1443 785 1578 12803
LPS.Response.final<-subset(LPS.Response, LPS.Response$Group != 4)
ggplotRegression(lm(LPS.Response.final$log2FC.MIA ~ LPS.Response.final$log2FC.CON, data=LPS.Response.final)) + geom_abline(slope=1) + labs(x="CON log2FC (LPS/SAL)", y="MIA log2FC (LPS/SAL)") + xlim(-10, 11) + ylim(-10, 11) + theme_cowplot()
Concordance Correlation Coefficient (CCC) of MIA and CON LPS Response in Str
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.702 0.686 0.717
Differential gene expression for design: Region + Trt + Region:Trt
dds <- DESeqDataSetFromTximport(txi.rsem, colData = metadata, design = ~ Region + Trt + Region:Trt)
# filter low expressed genes
keep <- rowSums(counts(dds)) >= 50
table(keep)
## keep
## FALSE TRUE
## 37355 18132
dds <- dds[keep,]
dds <- DESeq(dds)
resultsNames(dds)
## [1] "Intercept" "Region_Str_vs_Ctx" "Trt_LPS_vs_SAL"
## [4] "RegionStr.TrtLPS"
# add in gene symbols
ah <- AnnotationHub()
## snapshotDate(): 2020-10-27
edb <- ah[["AH53222"]]
## loading from cache
genes <- gsub("\\..*","",rownames(dds))
symbol <- mapIds(edb, keys=genes, column="SYMBOL", keytype="GENEID", multiVals="first")
## Warning: Unable to map 539 of 18132 requested IDs.
LvSinCTX <- results(dds, contrast=c("Trt","LPS","SAL"), alpha=0.05)
LvSinSTR <- results(dds, list( c("Trt_LPS_vs_SAL","RegionStr.TrtLPS")), alpha=0.05)
LvSinCTX$gene <- symbol
LvSinSTR$gene <- symbol
LPS.Response <- merge(as.data.frame(LvSinCTX[,c(7,1,2,3,4,5,6)]),as.data.frame(LvSinSTR[,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.CTX","lfcSE.CTX","stat.CTX","pval.CTX","padj.CTX","log2FC.STR","lfcSE.STR","stat.STR","pval.STR","padj.STR")
#Assign the Group variable for coloring 4-way plot
for (i in c(1:nrow(LPS.Response))){
if (LPS.Response$padj.CTX[i]<0.05 & (!is.na(LPS.Response$padj.CTX[i]))){
LPS.Response$CTXSig[i]='YES'
}
else {
LPS.Response$CTXSig[i]='NO'
}
}
for (i in c(1:nrow(LPS.Response))){
if (LPS.Response$padj.STR[i]<0.05 & (!is.na(LPS.Response$padj.STR[i]))){
LPS.Response$STRSig[i]='YES'
}
else {
LPS.Response$STRSig[i]='NO'
}
}
for (i in c(1:nrow(LPS.Response))){
if (LPS.Response$CTXSig[i]=='YES' & LPS.Response$STRSig[i]=='NO'){
LPS.Response$Group[i]=1 }
else if (LPS.Response$CTXSig[i]=='NO' & LPS.Response$STRSig[i]=='YES'){
LPS.Response$Group[i]=2}
else if (LPS.Response$CTXSig[i]=='YES' & LPS.Response$STRSig[i]=='YES'){
LPS.Response$Group[i]=3}
else {
LPS.Response$Group[i]=4}
}
table(LPS.Response$Group)
##
## 1 2 3 4
## 928 1412 3382 12410
LPS.Response.final<-subset(LPS.Response, LPS.Response$Group != 4)
ggplot(data = LPS.Response.final, aes(x = log2FC.CTX, y = log2FC.STR, color=factor(Group))) + geom_point() + scale_color_manual(values = c("hotpink", "mediumpurple1","purple"), labels = c("CTX only","STR only","CTX & STR")) + geom_abline(slope=1, size=0.5) + labs(x="CTX logFC",y="STR logFC", color="") + xlim(-7,11) + ylim(-7,11) + theme_cowplot() + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX>0 & log2FC.STR >0 & log2FC.STR > log2FC.CTX)) # 1985
## [1] 1985 16
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX>0 & log2FC.STR >0 & log2FC.STR < log2FC.CTX)) # 1002
## [1] 1002 16
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX<0 & log2FC.STR <0 & log2FC.STR > log2FC.CTX)) # 1251
## [1] 1251 16
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX<0 & log2FC.STR <0 & log2FC.STR < log2FC.CTX)) # 1334
## [1] 1334 16
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX>0 & log2FC.STR <0)) # 68
## [1] 68 16
dim(LPS.Response.final %>% dplyr::filter(log2FC.CTX<0 & log2FC.STR >0)) # 82
## [1] 82 16
ggplotRegression(lm(LPS.Response.final$log2FC.STR ~ LPS.Response.final$log2FC.CTX, data=LPS.Response.final)) + geom_abline(slope=1) + labs(x="CTX log2FC (LPS/SAL)", y="STR log2FC (LPS/SAL)") + xlim(-7,11) + ylim(-7,11) + theme_cowplot()
Concordance Correlation Coefficient (CCC) of MIA and CON LPS Response
library(epiR)
# concordance correlation coefficient
epi.ccc(LPS.Response.final$log2FC.CTX, LPS.Response.final$log2FC.STR, ci = "z-transform", conf.level = 0.95)$rho.c
## est lower upper
## 1 0.858 0.851 0.864