Load library and data
Human-rhesus heart- DE genes
- Human versus rhesus heart- DE genes
- Human versus rhesus heart- non DE genes
Human-rhesus kidney- DE genes
- Human versus rhesus kidney- DE genes
- Human versus rhesus kidney- non DE genes
Liver
- Human versus rhesus liver- DE genes
- Human versus rhesus liver- non DE genes
Lung
- Human versus rhesus lung- DE genes
- Human versus rhesus lung- non DE genes
Make into figure

Load library 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")

FDR_level = 0.05

# 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/250_exp_avg_methyl_hcr_4155_genes.txt", sep="")


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

Human-rhesus heart- DE genes

# Check parameters

human_chimp_heart <- c(32, 36, 40, 44, 20, 24, 28)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]
rownames(hc_exprs_heart) <- exprs[,1]

methyl <- exprs_methyl[,human_chimp_heart]
rownames(methyl) <- exprs[,1]


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    3889     266

Human versus rhesus 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.3378995

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 versus rhesus 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.005228758

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-rhesus kidney- DE genes

# Check parameters

human_chimp_heart <- c(33, 37, 41, 45, 17, 21, 25, 29)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]
rownames(hc_exprs_heart) <- exprs[,1]

methyl <- exprs_methyl[,human_chimp_heart]
rownames(methyl) <- exprs[,1]


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    4073      82

Human versus rhesus 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.3333333

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

combine_dataa <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sig), 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_sigc")

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

Human versus rhesus 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.005989583

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

combine_datab <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigd), 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_sigd")

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

Liver

# Check parameters

human_chimp_heart <- c(34, 38, 42, 46, 18, 22, 26, 30)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]
rownames(hc_exprs_heart) <- exprs[,1]

methyl <- exprs_methyl[,human_chimp_heart]
rownames(methyl) <- exprs[,1]


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    4057      98

Human versus rhesus 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.4021352

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

combine_datac <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sig), 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_sigc")

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

Human versus rhesus 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.008776458

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

combine_datad <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigd), 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_sigd")

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

Lung

# Check parameters

human_chimp_heart <- c(35, 39, 43, 47, 19, 23, 27, 31)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]
rownames(hc_exprs_heart) <- exprs[,1]

methyl <- exprs_methyl[,human_chimp_heart]
rownames(methyl) <- exprs[,1]


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    4061      94

Human versus rhesus 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.2151394

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

combine_datae <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sig), 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_sigc")

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

Human versus rhesus 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.0002939447

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

combine_dataf <- as.data.frame(cbind(fit_ash$result$betahat, fit_ash$result$sebetahat, ind_sigd), 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_sigd")

hr_def <- ggplot(combine_dataf, aes(combine_dataf$Effect, combine_dataf$SE)) + geom_point(color = c(combine_dataf$ind_sigd)) + 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("11A. Heart DE Genes", "11B. Heart non-DE genes", "11C. Kidney DE Genes", "11D. Kidney non-DE genes", "11E. Liver DE genes", "11F. Liver non-DE genes", "11G. Lung DE genes", "11H. Lung non-DE genes", label_y = 0), ncol = 3, nrow = 5)

## 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/test11b.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
           )

Figure_Supplementary11

Lauren Blake

December 2, 2018

Load library and data

Human-rhesus heart- DE genes

Human versus rhesus heart- DE genes

Human versus rhesus heart- non DE genes

Human-rhesus kidney- DE genes

Human versus rhesus kidney- DE genes

Human versus rhesus kidney- non DE genes

Liver

Human versus rhesus liver- DE genes

Human versus rhesus liver- non DE genes

Lung

Human versus rhesus lung- DE genes

Human versus rhesus lung- non DE genes

Make into figure