Install libraries and data
Set FDR level/svalue
Human-chimp heart
Human-chimp heart- DE genes
Human-chimp heart- non DE genes
Human-chimp kidney
Human-chimp kidney- DE genes
Human-chimp kidney- non DE genes
Human-chimp liver
Human-chimp liver- DE genes
Human-chimp liver- non DE genes
Human-chimp lung
Human-chimp lung- DE genes
Human-chimp lung- non DE genes
Make into figure

Install libraries and data

# Load libraries
library("edgeR")

## Loading required package: limma

library("limma")
library("plyr")
library("ashr")

## Warning: package 'ashr' was built under R version 3.4.4

library("cowplot")

## Warning: package 'cowplot' was built under R version 3.4.4

## Loading required package: ggplot2

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggplot2':
## 
##     ggsave

library("vashr")

## Loading required package: SQUAREM

## Loading required package: qvalue

library("devtools")
#devtools::install_github("jhsiao999/mediation/pkg")
library("medinome")

# import sample labels
samples <- read.delim("../../../Reg_Evo_Primates/data/Sample_info_RNAseq_RIN.txt")
RIN <- samples$RIN
# expression
exprs <- read.table("../../../Reg_Evo_Primates/data/human_chimp_orth_exp_methyl_7725_hum.txt", sep="")
# methylation data
methyl <- read.csv("../../../Reg_Evo_Primates/data/chimp_human_orth_7725_avg_methyl_per_ts_gene.txt", sep="")

# Normalized gene expression data
cpm.voom.cyclic <- readRDS("../../../Reg_Evo_Primates/data/human_chimp_orth_cpm_voom_cyclic.rds")

Set FDR level/svalue

FDR_level <- 0.05

Human-chimp heart

# Check parameters

human_chimp_heart <- c(1, 5, 9, 13, 20, 24, 28)

# Methylation
hc_exprs <- exprs[,2:48]
exprs_methyl <- exprs[,49:79]

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


methyl <- exprs_methyl[,human_chimp_heart]



Y <- hc_exprs_heart
two_species <- samples$Species[human_chimp_heart]
X <- droplevels.factor(two_species)
M <- methyl
RIN_subset <- RIN[human_chimp_heart]



check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

## Warning: package 'assertthat' was built under R version 3.4.4

fit <- vash(check_values$d_se,df=5, singlecomp = T, unimodal = "auto") 
fit_ash_normal <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")
summary(fit_ash_normal$result$svalue < 0.05)

##    Mode   FALSE    TRUE 
## logical    7688      37

Human-chimp heart- DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val < FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:7]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:7]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=5, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0.1141498

ind_sigg <- (fit_ash$result$svalue < 0.05)
ind_sigg[ind_sigg == TRUE] <- "red"
ind_sigg[ind_sigg == FALSE] <- "black"

combine_datag <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigg), stringsAsFactors = FALSE)

combine_datag[,1] <- as.numeric(combine_datag[,1])
combine_datag[,2] <- as.numeric(combine_datag[,2])
colnames(combine_datag) <- c("Effect", "SE", "ind_sigg")

hr_deg <- ggplot(combine_datag, aes(combine_datag$Effect, combine_datag$SE)) + geom_point(color = c(combine_datag$ind_sigg)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp heart- non DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val > FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:7]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:7]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=5, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0

ind_sigh <- (fit_ash$result$svalue < 0.05)
ind_sigh[ind_sigh == TRUE] <- "red"
ind_sigh[ind_sigh == FALSE] <- "black"

combine_datah <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigh), stringsAsFactors = FALSE)

combine_datah[,1] <- as.numeric(combine_datah[,1])
combine_datah[,2] <- as.numeric(combine_datah[,2])
colnames(combine_datah) <- c("Effect", "SE", "ind_sigh")

hr_deh <- ggplot(combine_datah, aes(combine_datah$Effect, combine_datah$SE)) + geom_point(color = c(combine_datah$ind_sigh)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp kidney

# Check parameters

human_chimp_heart <- c(2, 6, 10, 14, 17, 21, 25, 29)

# Methylation
hc_exprs <- exprs[,2:48]
exprs_methyl <- exprs[,49:79]

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


methyl <- exprs_methyl[,human_chimp_heart]



Y <- hc_exprs_heart
two_species <- samples$Species[human_chimp_heart]
X <- droplevels.factor(two_species)
M <- methyl
RIN_subset <- RIN[human_chimp_heart]



check_values <- mediate.test.regressing(Y, X, M, RIN_subset)
fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash_normal <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")
summary(fit_ash_normal$result$svalue < 0.05)

##    Mode   FALSE    TRUE 
## logical    7711      14

Human-chimp kidney- DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val < FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0.1222571

ind_siga <- (fit_ash$result$svalue < 0.05)
ind_siga[ind_siga == TRUE] <- "red"
ind_siga[ind_siga == FALSE] <- "black"

combine_dataa <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_siga), stringsAsFactors = FALSE)

combine_dataa[,1] <- as.numeric(combine_dataa[,1])
combine_dataa[,2] <- as.numeric(combine_dataa[,2])
colnames(combine_dataa) <- c("Effect", "SE", "ind_siga")

hr_dea <- ggplot(combine_dataa, aes(combine_dataa$Effect, combine_dataa$SE)) + geom_point(color = c(combine_dataa$ind_siga)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp kidney- non DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val > FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0

ind_sigb <- (fit_ash$result$svalue < 0.05)
ind_sigb[ind_sigb == TRUE] <- "red"
ind_sigb[ind_sigb == FALSE] <- "black"

combine_datab <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigb), stringsAsFactors = FALSE)

combine_datab[,1] <- as.numeric(combine_datab[,1])
combine_datab[,2] <- as.numeric(combine_datab[,2])
colnames(combine_datab) <- c("Effect", "SE", "ind_sigb")

hr_deb <- ggplot(combine_datab, aes(combine_datab$Effect, combine_datab$SE)) + geom_point(color = c(combine_datab$ind_sigb)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp liver

# Check parameters

human_chimp_heart <- c(3, 7, 11, 15, 18, 22, 26, 30)

# Methylation
hc_exprs <- exprs[,2:48]
exprs_methyl <- exprs[,49:79]

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


methyl <- exprs_methyl[,human_chimp_heart]



Y <- hc_exprs_heart
two_species <- samples$Species[human_chimp_heart]
X <- droplevels.factor(two_species)
M <- methyl
RIN_subset <- RIN[human_chimp_heart]



check_values <- mediate.test.regressing(Y, X, M, RIN_subset)
fit <- vash(check_values$d_se,df=5, singlecomp = T, unimodal = "auto") 
fit_ash_normal <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")
summary(fit_ash_normal$result$svalue < 0.05)

##    Mode   FALSE    TRUE 
## logical    7702      23

Human-chimp liver- DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val < FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0.2562814

ind_siga <- (fit_ash$result$svalue < 0.05)
ind_siga[ind_siga == TRUE] <- "red"
ind_siga[ind_siga == FALSE] <- "black"

combine_datac <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_siga), stringsAsFactors = FALSE)

combine_datac[,1] <- as.numeric(combine_datac[,1])
combine_datac[,2] <- as.numeric(combine_datac[,2])
colnames(combine_datac) <- c("Effect", "SE", "ind_siga")

hr_dec <- ggplot(combine_datac, aes(combine_datac$Effect, combine_datac$SE)) + geom_point(color = c(combine_datac$ind_siga)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp liver- non DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val > FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0

ind_sigb <- (fit_ash$result$svalue < 0.05)
ind_sigb[ind_sigb == TRUE] <- "red"
ind_sigb[ind_sigb == FALSE] <- "black"

combine_datad <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigb), stringsAsFactors = FALSE)

combine_datad[,1] <- as.numeric(combine_datad[,1])
combine_datad[,2] <- as.numeric(combine_datad[,2])
colnames(combine_datad) <- c("Effect", "SE", "ind_sigb")

hr_ded <- ggplot(combine_datad, aes(combine_datad$Effect, combine_datad$SE)) + geom_point(color = c(combine_datad$ind_sigb)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp lung

# Check parameters

human_chimp_heart <- c(4, 8, 12, 16, 19, 23, 27, 31)

# Methylation
hc_exprs <- exprs[,2:48]
exprs_methyl <- exprs[,49:79]

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


methyl <- exprs_methyl[,human_chimp_heart]



Y <- hc_exprs_heart
two_species <- samples$Species[human_chimp_heart]
X <- droplevels.factor(two_species)
M <- methyl
RIN_subset <- RIN[human_chimp_heart]



check_values <- mediate.test.regressing(Y, X, M, RIN_subset)
fit <- vash(check_values$d_se,df=5, singlecomp = T, unimodal = "auto") 
fit_ash_normal <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")
summary(fit_ash_normal$result$svalue < 0.05)

##    Mode   FALSE    TRUE 
## logical    7712      13

Human-chimp lung- DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val < FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0.1155378

ind_siga <- (fit_ash$result$svalue < 0.05)
ind_siga[ind_siga == TRUE] <- "red"
ind_siga[ind_siga == FALSE] <- "black"

combine_datae <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_siga), stringsAsFactors = FALSE)

combine_datae[,1] <- as.numeric(combine_datae[,1])
combine_datae[,2] <- as.numeric(combine_datae[,2])
colnames(combine_datae) <- c("Effect", "SE", "ind_siga")

hr_dee <- ggplot(combine_datae, aes(combine_datae$Effect, combine_datae$SE)) + geom_point(color = c(combine_datae$ind_siga)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Human-chimp lung- non DE genes

# Run the linear model in limma
  design_1 <- model.matrix(~ X + RIN_subset)
  fit_all <- lmFit(cpm.voom.cyclic[,human_chimp_heart], design_1)
  fit_all <- eBayes(fit_all)

# Get results
  HvC_Heart_fit_all = topTable(fit_all, coef=2, adjust="BH", number=Inf, sort.by="none")
  HvC_Heart_fit_all_5perc <-
HvC_Heart_fit_all[which(HvC_Heart_fit_all$adj.P.Val > FDR_level), ]
  
 human_chimp_heart1 <-  rownames(methyl) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(methyl, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
M <- counts_human_chimp_heart_subset

human_chimp_heart1 <-  rownames(hc_exprs_heart) %in% rownames( HvC_Heart_fit_all_5perc)
human_chimp_heart1 <- as.data.frame(human_chimp_heart1)
counts_genes_in <- cbind(hc_exprs_heart, human_chimp_heart1)
counts_genes_in_cutoff <- subset(counts_genes_in, human_chimp_heart1 == "TRUE")
counts_human_chimp_heart_subset <- counts_genes_in_cutoff[,1:8]
Y <- counts_human_chimp_heart_subset

check_values <- mediate.test.regressing(Y, X, M, RIN_subset)

fit <- vash(check_values$d_se,df=6, singlecomp = T, unimodal = "auto") 
fit_ash <- ash(as.vector(check_values$d), fit$sd.post, mode = 0, mixcompdist = "normal")

calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)

## [1] 0

ind_sigb <- (fit_ash$result$svalue < 0.05)
ind_sigb[ind_sigb == TRUE] <- "red"
ind_sigb[ind_sigb == FALSE] <- "black"

combine_dataf <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigb), stringsAsFactors = FALSE)

combine_dataf[,1] <- as.numeric(combine_dataf[,1])
combine_dataf[,2] <- as.numeric(combine_dataf[,2])
colnames(combine_dataf) <- c("Effect", "SE", "ind_sigb")

hr_def <- ggplot(combine_dataf, aes(combine_dataf$Effect, combine_dataf$SE)) + geom_point(color = c(combine_dataf$ind_sigb)) + ylab("Standard error") + xlab("Effect size difference") + ylim(0, 6.5) + xlim(-8, 8)

Make into figure

eight_plots <- plot_grid(hr_deg, hr_deh, hr_dea, hr_deb, hr_dec, hr_ded, hr_dee, hr_def, labels = c("10A. Heart DE Genes", "10. Heart non-DE genes", "10C. Kidney DE Genes", "10D. Kidney non-DE genes", "10E. Liver DE genes", "10F. Liver non-DE genes", "10G. Lung DE genes", "10H. Lung non-DE genes", label_y = 0), ncol = 3, nrow = 5)

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

plot_fig4 <- plot_grid(eight_plots, ncol = 1, rel_heights=c(0.1, 1))

save_plot("/Users/laurenblake/Desktop/tissue_paper/test11.pdf", plot_fig4,
           ncol = 3, # we're saving a grid plot of 2 columns
           nrow = 5, # and 2 rows
           # each individual subplot should have an aspect ratio of 1.3
           base_aspect_ratio = 1
           )

Final_effect_size_figure_S10

Lauren Blake

December 2, 2018

Install libraries and data

Set FDR level/svalue

Human-chimp heart

Human-chimp heart- DE genes

Human-chimp heart- non DE genes

Human-chimp kidney

Human-chimp kidney- DE genes

Human-chimp kidney- non DE genes

Human-chimp liver

Human-chimp liver- DE genes

Human-chimp liver- non DE genes

Human-chimp lung

Human-chimp lung- DE genes

Human-chimp lung- non DE genes

Make into figure