getwd()
## [1] "/Users/laurenblake/Desktop/Regulatory_Evol/ashlar-trial/analysis"
# 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")


# Load libraries
library("edgeR")
## Loading required package: limma
library("limma")
library("plyr")
library("ashr")
library("cowplot")
## Warning: package 'cowplot' was built under R version 3.4.4
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.4.4
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
library("vashr")
## Loading required package: SQUAREM
## Loading required package: qvalue
library("devtools")
## Warning: package 'devtools' was built under R version 3.4.4
#devtools::install_github("jhsiao999/mediation/pkg")
library("medinome")

Set FDR level/svalue

FDR_level <- 0.1
#FDR_add <- 0.10

Rhesus heart-kidney

# Check parameters

human_chimp_heart <- c(32, 33, 36, 37, 40, 41, 44, 45)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3665     490

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.2034062

Human heart versus kidney- conserved DE genes

# Species no DE
species_no_DE <- c(17, 20, 21, 24, 25, 28, 29)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     586     391
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.4002047
# Species no DE
species_no_DE <- c(17, 20, 21, 24, 25, 28, 29)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    1879     256
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1199063

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0

Chimp heart-liver

# Check parameters

human_chimp_heart <- c(32, 34, 36, 38, 40, 42, 44, 46)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3529     626

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.2702249

Human heart versus kidney- conserved DE genes

# Species no DE
species_no_DE <- c(18, 20, 22, 24, 26, 28, 30)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     482     482
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.5
# Species no DE
species_no_DE <- c(18, 20, 22, 24, 26, 28, 30)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    1692     323
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1602978

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0

Chimp heart-lung

# Check parameters

human_chimp_heart <- c(32, 35, 36, 39, 40, 43, 44, 47)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3623     532

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.2430625

Human heart versus lung- conserved DE genes

# Species no DE
species_no_DE <- c(19, 20, 23, 24, 27, 28, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     537     390
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.420712
# Species no DE
species_no_DE <-  c(19, 20, 23, 24, 27, 28, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    1723     341
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1652132

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0

Chimp kidney-liver

# Check parameters

human_chimp_heart <- c(33, 34, 37, 38, 41, 42, 45, 46)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3786     369

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1690838

Human heart versus lung- conserved DE genes

# Species no DE
species_no_DE <- c(17, 18, 21, 22, 25, 26, 29, 30)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     660     314
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.3223819
# Species no DE
species_no_DE <-  c(17, 18, 21, 22, 25, 26, 29, 30)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    1677     176
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.09498111

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.000753012

Chimp kidney-lung

# Check parameters

human_chimp_heart <- c(33, 35, 37, 39, 41, 43, 45, 47)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3686     469

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.2034121

Human heart versus lung- conserved DE genes

# Species no DE
species_no_DE <- c(17, 19, 21, 23, 25, 27, 29, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     340     282
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.4533762
# Species no DE
species_no_DE <-  c(17, 19, 21, 23, 25, 27, 29, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    2090     336
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1384996

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0

Chimp liver-lung

# Check parameters

human_chimp_heart <- c(34, 35, 38, 39, 42, 43, 46, 47)

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

# HC heart only


hc_exprs_heart <- hc_exprs[,human_chimp_heart]


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


Y <- hc_exprs_heart
two_species <- samples$Tissue[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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash_normal$result$svalue < 0.05)
##    Mode   FALSE    TRUE 
## logical    3564     591

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.2383387

Human heart versus lung- conserved DE genes

# Species no DE
species_no_DE <- c(18, 19, 22, 23, 26, 27, 30, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] <  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical     534     465
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.4654655
# Species no DE
species_no_DE <-  c(18, 19, 22, 23, 26, 27, 30, 31)

species_no_DE_exprs <- hc_exprs[,species_no_DE]

species_no_DE_methyl <- exprs_methyl[,species_no_DE]

two_species <- samples$Tissue[species_no_DE]
no_DE_tissue <- droplevels.factor(two_species)

no_DE_RIN <- RIN[species_no_DE]

design <- model.matrix(~ as.factor(no_DE_tissue) + no_DE_RIN)
  fit_species_no_DE <- lmFit(cpm.voom.cyclic[,species_no_DE], design)
  fit_species_no_DE <- eBayes(fit_species_no_DE)

  HvC_Heart_fit_species_no_DE = topTable(fit_species_no_DE, coef=2, adjust="BH", number=Inf,
sort.by="none")

# Species DE 
  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_species_no_DE = cbind(HvC_Heart_fit_species_no_DE, HvC_Heart_fit_all)
  
 
      HvC_Heart_fit_all_5perc <- HvC_Heart_fit_species_no_DE[which(HvC_Heart_fit_species_no_DE[,6] >  FDR_level & HvC_Heart_fit_species_no_DE[,13] < 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")
## Due to absence of package REBayes, switching to EM algorithm
summary(fit_ash$result$svalue < FDR_level)
##    Mode   FALSE    TRUE 
## logical    1841     290
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0.1360863

Human heart versus 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")
## Due to absence of package REBayes, switching to EM algorithm
calc <- as.array(summary(fit_ash$result$svalue < FDR_level))
calc2 <- as.numeric(calc[[2]])

(nrow(Y)-calc2)/nrow(Y)
## [1] 0