# 1 bp overlap
heart_liver_1 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 1E-9 -F 1E-9")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e32e22332f.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e36ed5d06b.bed -wao -f 1E-9 -F 1E-9

summary(as.integer(heart_liver_1[,9]) > 0)

   Mode   FALSE    TRUE 
logical   16965    5607 

# 10% overlap
heart_liver_10 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.10 -F 0.10")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3b81f65f.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3594ca6a6.bed -wao -f 0.10 -F 0.10

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17096    5472 

# 20% overlap
heart_liver_20 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.20 -F 0.20")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e356146bdf.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e32d3d1a5e.bed -wao -f 0.20 -F 0.20

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17374    5189 

# 25% overlap
heart_liver_25 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.25 -F 0.25")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e321e33fe7.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3504c6f64.bed -wao -f 0.25 -F 0.25

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17554    5007 

# 30% overlap
heart_liver_30 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.30 -F 0.30")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e348ca2905.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3707f3352.bed -wao -f 0.30 -F 0.30

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17737    4824 

#40% overlap
heart_liver_40 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.4 -F 0.4")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e316ec6cc5.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e379a1c588.bed -wao -f 0.4 -F 0.4

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   18161    4400 

# 50% overlap
heart_liver_50 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.50 -F 0.50")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35af7adf6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3646a08f0.bed -wao -f 0.50 -F 0.50

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   18703    3858 

#60% overlap
heart_liver_60 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.6 -F 0.6")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3b527156.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e34a0fee0e.bed -wao -f 0.6 -F 0.6

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   19286    3275 

# 70% overlap
heart_liver_70 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.70 -F 0.70")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e31f6a6b32.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e319b55a59.bed -wao -f 0.70 -F 0.70

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   19925    2636 

#75% overlap
heart_liver_75 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.75 -F 0.75")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e32c794932.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e364556bbc.bed -wao -f 0.75 -F 0.75

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   20256    2305 

# 80% overlap
heart_liver_80 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.80 -F 0.80")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e33ffeacc6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e34cda51f0.bed -wao -f 0.80 -F 0.80

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   20597    1964 

#90% overlap
heart_liver_90 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.9 -F 0.9")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e320beb9cd.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e330ca4304.bed -wao -f 0.9 -F 0.9

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   21223    1338 

#100% overlap
heart_liver_100 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 1 -F 1")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e34d3333b5.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e367dc9ea8.bed -wao -f 1 -F 1

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   21587     974 

Plot the overlap

# Find the number of tDMRs depending on the bp overlap
summary(as.integer(heart_liver_1[,9]) > 0)

   Mode   FALSE    TRUE 
logical   16965    5607 

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17096    5472 

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17374    5189 

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17554    5007 

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   17737    4824 

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   18161    4400 

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   18703    3858 

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   19286    3275 

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   19925    2636 

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   20256    2305 

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   20597    1964 

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   21223    1338 

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   21587     974 

# Plot the number of tDMRs
conserved_tdmr_heart_liver_heart_lung <- c(5607, 5472, 5189, 5007, 4824, 4400, 3858, 3275, 2636, 2305,1964,1338, 974)
perc_overlap <- c(0.000000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung, xlab = "Percentage of overlap (bp)", ylab = "Number of human tDMRs (HvLi and HvLu)",  pch = 19)

# Investigate the size of the tDMRs depending on the percentage overlap

check_1 <- heart_liver_10[which((as.integer(heart_liver_1[,9]) > 0) == TRUE),]
check_10 <- heart_liver_10[which((as.integer(heart_liver_10[,9]) > 0) == TRUE),]
check_20 <- heart_liver_20[which((as.integer(heart_liver_20[,9]) > 0) == TRUE),]
check_25 <- heart_liver_25[which((as.integer(heart_liver_25[,9]) > 0) == TRUE),]
check_30 <- heart_liver_30[which((as.integer(heart_liver_30[,9]) > 0) == TRUE),]
check_40 <- heart_liver_40[which((as.integer(heart_liver_40[,9]) > 0) == TRUE),]
check_50 <- heart_liver_50[which((as.integer(heart_liver_50[,9]) > 0) == TRUE),]
check_60 <- heart_liver_60[which((as.integer(heart_liver_60[,9]) > 0) == TRUE),]
check_70 <- heart_liver_70[which((as.integer(heart_liver_70[,9]) > 0) == TRUE),]
check_75 <- heart_liver_75[which((as.integer(heart_liver_75[,9]) > 0) == TRUE),]
check_80 <- heart_liver_80[which((as.integer(heart_liver_80[,9]) > 0) == TRUE),]
check_90 <- heart_liver_90[which((as.integer(heart_liver_90[,9]) > 0) == TRUE),]
check_100 <- heart_liver_100[which((as.integer(heart_liver_100[,9]) > 0) == TRUE),]

summary(as.integer(check_1[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0     0.0   120.0   280.7   468.8  3167.0       1 

summary(as.integer(check_10[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   228.0   419.0   501.7   697.0  3167.0 

summary(as.integer(check_20[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      6     251     440     522     717    3167 

summary(as.integer(check_25[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   265.0   454.0   534.6   727.0  3167.0 

summary(as.integer(check_30[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   275.8   470.0   546.2   738.0  3167.0 

summary(as.integer(check_40[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   300.0   500.0   570.6   764.0  3167.0 

summary(as.integer(check_50[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   323.0   533.0   598.4   806.0  3167.0 

summary(as.integer(check_60[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   344.0   574.0   627.7   845.0  3167.0 

summary(as.integer(check_70[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   350.0   601.0   649.6   885.5  3167.0 

summary(as.integer(check_75[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   346.0   604.0   655.4   904.0  3167.0 

summary(as.integer(check_80[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   332.8   599.0   650.3   907.8  3167.0 

summary(as.integer(check_90[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   276.5   522.0   601.2   839.0  3167.0 

summary(as.integer(check_100[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   223.2   410.5   502.5   707.8  2658.0 

# Plot the median tDMR length
conserved_tdmr_heart_liver_heart_lung_bp <- c(120, 419, 440, 454, 470, 500, 533, 574, 601, 604,599,522, 410.5)
perc_overlap <- c(0.000000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung_bp, xlab = "Percentage of overlap (bp)", ylab = "Median length of overlap (bp)", pch = 19, main = "Median length of overlap between human tDMRs (HvLi and HvLu)")

Chimp heart versus liver and chimp heart versus lung

# columns to pull
pull_col <- c(1:3, 16)

# make bed style format: change name from human to rhesus
humans_heart_liver_DMRs_bed <- chimps_heart_liver_DMRs[,pull_col]

# Get DMRs on the autosomal chromosomes only
humans_heart_liver_DMRs_bed_noX <- humans_heart_liver_DMRs_bed[which(humans_heart_liver_DMRs_bed[,1] != "chrX"),]

humans_heart_liver_DMRs_bed_noX[,1] <- as.character(humans_heart_liver_DMRs_bed_noX[,1])
humans_heart_liver_DMRs_bed_noX[,2] <- as.integer(humans_heart_liver_DMRs_bed_noX[,2])
humans_heart_liver_DMRs_bed_noX[,3] <- as.integer(humans_heart_liver_DMRs_bed_noX[,3])

summary(humans_heart_liver_DMRs_bed_noX[,3] - humans_heart_liver_DMRs_bed_noX[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0   266.0   498.0   617.6   849.0  4950.0 

sort_humans_heart_liver_DMRs_bed <- bedr.sort.region(humans_heart_liver_DMRs_bed_noX, check.chr = FALSE)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.

# make bed style format
humans_heart_lung_DMRs_bed <- chimps_heart_lung_DMRs[,pull_col]
humans_heart_lung_DMRs_bed_noX <- humans_heart_lung_DMRs_bed[which(humans_heart_lung_DMRs_bed[,1] != "chrX"),]

humans_heart_lung_DMRs_bed_noX[,1] <- as.character(humans_heart_lung_DMRs_bed_noX[,1])
humans_heart_lung_DMRs_bed_noX[,2] <- as.integer(humans_heart_lung_DMRs_bed_noX[,2])
humans_heart_lung_DMRs_bed_noX[,3] <- as.integer(humans_heart_lung_DMRs_bed_noX[,3])
sort_humans_heart_lung_DMRs_bed <- bedr.sort.region(humans_heart_lung_DMRs_bed_noX, check.chr = FALSE)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.

summary(humans_heart_lung_DMRs_bed_noX[,3] - humans_heart_lung_DMRs_bed_noX[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0   225.0   413.5   512.3   708.0  4186.0 

heart_lung <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.5 -F 0.5", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e322cb955.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e341d5deed.bed -wao -f 0.5 -F 0.5

summary(heart_lung[,3] - heart_lung[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0   266.0   498.0   617.6   849.0  4950.0 

summary(as.integer(heart_lung[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00    0.00   56.97    0.00 4186.00 

summary(as.integer(heart_lung[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26162    2607 

lung_heart <- bedr(input = list(a = sort_humans_heart_lung_DMRs_bed, b = sort_humans_heart_liver_DMRs_bed), method = "intersect",  params = "-wao -f 0.5 -F 0.5", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36d3b12c2.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e342e44a79.bed -wao -f 0.5 -F 0.5

summary(as.integer(lung_heart[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0     0.0   192.2   278.8  4186.0 

summary(as.integer(lung_heart[,9]) > 0)

   Mode   FALSE    TRUE 
logical    5919    2607 

# Make a Venn diagram based on the results

summary(as.integer(heart_lung[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26162    2607 

summary(as.integer(lung_heart[,9]) > 0)

   Mode   FALSE    TRUE 
logical    5919    2607 

grid.newpage()
draw.pairwise.venn(area1 = (26162+2607), area2 = (5919+2607), cross.area = 2607, category = c("Heart v. Liver", "Heart v. Lung"))

(polygon[GRID.polygon.20], polygon[GRID.polygon.21], polygon[GRID.polygon.22], polygon[GRID.polygon.23], text[GRID.text.24], text[GRID.text.25], text[GRID.text.26], text[GRID.text.27], text[GRID.text.28]) 

# 1 bp overlap
heart_liver_1 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.00000001 -F 0.00000001", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35d2cf2c5.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e35f708f26.bed -wao -f 0.00000001 -F 0.00000001

# 10% overlap
heart_liver_10 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.10 -F 0.10", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e362de3e60.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e319fed68.bed -wao -f 0.10 -F 0.10

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   24943    3837 

# 20% overlap
heart_liver_20 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.20 -F 0.20", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3aa4640d.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e365e71b74.bed -wao -f 0.20 -F 0.20

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25160    3614 

# 25% overlap
heart_liver_25 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.25 -F 0.25", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e31443c3f0.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e354f63da3.bed -wao -f 0.25 -F 0.25

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25301    3472 

# 30% overlap
heart_liver_30 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.30 -F 0.30", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e365db5c2c.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e369495f96.bed -wao -f 0.30 -F 0.30

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25439    3332 

#40% overlap
heart_liver_40 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.4 -F 0.4", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36f32eaea.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3794a197e.bed -wao -f 0.4 -F 0.4

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25795    2974 

# 50% overlap
heart_liver_50 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.50 -F 0.50", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36da7e2a6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3571beacc.bed -wao -f 0.50 -F 0.50

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26162    2607 

#60% overlap
heart_liver_60 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.6 -F 0.6", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e332110690.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3379becef.bed -wao -f 0.6 -F 0.6

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26589    2180 

# 70% overlap
heart_liver_70 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.70 -F 0.70", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3f3a284f.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3775bb6e7.bed -wao -f 0.70 -F 0.70

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27029    1740 

#75% overlap
heart_liver_75 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.75 -F 0.75", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e34f24d4a.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3f40eca.bed -wao -f 0.75 -F 0.75

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27267    1502 

# 80% overlap
heart_liver_80 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.80 -F 0.80", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e31beeaa88.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3435021e.bed -wao -f 0.80 -F 0.80

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27515    1254 

#90% overlap
heart_liver_90 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.9 -F 0.9", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e374fabe7b.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e342d70593.bed -wao -f 0.9 -F 0.9

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27915     854 

#100% overlap
heart_liver_100 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 1 -F 1", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36e98ab75.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e316d02663.bed -wao -f 1 -F 1

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   28148     621 

# Find the number of tDMRs depending on the bp overlap
summary(as.integer(heart_liver_1[,9]) > 0)

   Mode   FALSE    TRUE 
logical   24854    3940 

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   24943    3837 

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25160    3614 

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25301    3472 

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25439    3332 

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   25795    2974 

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26162    2607 

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   26589    2180 

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27029    1740 

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27267    1502 

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27515    1254 

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   27915     854 

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   28148     621 

# Plot the number of tDMRs
conserved_tdmr_heart_liver_heart_lung <- c(621, 3837, 3614, 3472, 3332, 2974, 2607, 2180, 1740, 1502,1254,854, 621)
perc_overlap <- c(0.00000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung, xlab = "Percentage of overlap (bp)", ylab = "Number of overlapping chimp tDMRs",  pch = 19, main = "Overlapping Chimp tDMRs (HvLi and HvLu)")

# Investigate the size of the tDMRs depending on the percentage overlap

check_1 <- heart_liver_1[which((as.integer(heart_liver_1[,9]) > 0) == TRUE),]
check_10 <- heart_liver_10[which((as.integer(heart_liver_10[,9]) > 0) == TRUE),]
check_20 <- heart_liver_20[which((as.integer(heart_liver_20[,9]) > 0) == TRUE),]
check_25 <- heart_liver_25[which((as.integer(heart_liver_25[,9]) > 0) == TRUE),]
check_30 <- heart_liver_30[which((as.integer(heart_liver_30[,9]) > 0) == TRUE),]
check_40 <- heart_liver_40[which((as.integer(heart_liver_40[,9]) > 0) == TRUE),]
check_50 <- heart_liver_50[which((as.integer(heart_liver_50[,9]) > 0) == TRUE),]
check_60 <- heart_liver_60[which((as.integer(heart_liver_60[,9]) > 0) == TRUE),]
check_70 <- heart_liver_70[which((as.integer(heart_liver_70[,9]) > 0) == TRUE),]
check_75 <- heart_liver_75[which((as.integer(heart_liver_75[,9]) > 0) == TRUE),]
check_80 <- heart_liver_80[which((as.integer(heart_liver_80[,9]) > 0) == TRUE),]
check_90 <- heart_liver_90[which((as.integer(heart_liver_90[,9]) > 0) == TRUE),]
check_100 <- heart_liver_100[which((as.integer(heart_liver_100[,9]) > 0) == TRUE),]

summary(as.integer(check_1[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    5.0   228.0   423.5   509.1   703.0  4186.0 

summary(as.integer(check_10[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    9.0   242.0   433.0   520.8   709.0  4186.0 

summary(as.integer(check_20[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   270.0   459.5   544.2   733.0  4186.0 

summary(as.integer(check_25[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   280.0   477.5   557.2   745.0  4186.0 

summary(as.integer(check_30[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   291.0   491.0   570.1   761.2  4186.0 

summary(as.integer(check_40[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   316.0   533.5   599.8   793.0  4186.0 

summary(as.integer(check_50[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   343.0   560.0   628.7   834.0  4186.0 

summary(as.integer(check_60[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   363.8   597.0   660.2   873.5  4186.0 

summary(as.integer(check_70[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   361.0   608.0   678.2   911.0  4186.0 

summary(as.integer(check_75[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   351.2   612.0   679.6   915.8  4186.0 

summary(as.integer(check_80[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   327.8   589.0   677.2   921.0  4186.0 

summary(as.integer(check_90[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   277.2   537.5   626.8   863.8  4186.0 

summary(as.integer(check_100[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   15.0   230.0   428.0   528.9   740.0  2637.0 

# Plot the median tDMR length
conserved_tdmr_heart_liver_heart_lung_bp <- c(301, 433, 459.5, 477.5, 491, 533.5, 560, 597, 608, 612,589,537.5, 428)
perc_overlap <- c(0.000000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung_bp, xlab = "Percentage of overlap (bp)", ylab = "Median length of overlap (bp)", pch = 19, main = "Median length of overlap between chimp tDMRs (HvLi and HvLu)")

Rhesus heart versus liver and rhesus heart versus lung

# Heart versus liver

# columns to pull
pull_col <- c(1:3, 16)

# make bed style format: change name from human to rhesus
humans_heart_liver_DMRs_bed <- rhesus_heart_liver_DMRs[,pull_col]

# Remove chr17_gl000204_random
humans_heart_liver_DMRs_bed_noX <- humans_heart_liver_DMRs_bed[which(humans_heart_liver_DMRs_bed[,1] != "chr17_gl000204_random"),]

humans_heart_liver_DMRs_bed_noX[,1] <- as.character(humans_heart_liver_DMRs_bed_noX[,1])
humans_heart_liver_DMRs_bed_noX[,2] <- as.integer(humans_heart_liver_DMRs_bed_noX[,2])
humans_heart_liver_DMRs_bed_noX[,3] <- as.integer(humans_heart_liver_DMRs_bed_noX[,3])

# Remove chr11_gl000202_random

humans_heart_liver_DMRs_bed_noX <- humans_heart_liver_DMRs_bed_noX[which(humans_heart_liver_DMRs_bed_noX[,1] != "chr11_gl000202_random"),]

humans_heart_liver_DMRs_bed_noX[,1] <- as.character(humans_heart_liver_DMRs_bed_noX[,1])
humans_heart_liver_DMRs_bed_noX[,2] <- as.integer(humans_heart_liver_DMRs_bed_noX[,2])
humans_heart_liver_DMRs_bed_noX[,3] <- as.integer(humans_heart_liver_DMRs_bed_noX[,3])


# Get DMRs on the autosomal chromosomes only
humans_heart_liver_DMRs_bed_noX <- humans_heart_liver_DMRs_bed_noX[which(humans_heart_liver_DMRs_bed_noX[,1] != "chrX"),]

humans_heart_liver_DMRs_bed_noX[,1] <- as.character(humans_heart_liver_DMRs_bed_noX[,1])
humans_heart_liver_DMRs_bed_noX[,2] <- as.integer(humans_heart_liver_DMRs_bed_noX[,2])
humans_heart_liver_DMRs_bed_noX[,3] <- as.integer(humans_heart_liver_DMRs_bed_noX[,3])

summary(humans_heart_liver_DMRs_bed_noX[,3] - humans_heart_liver_DMRs_bed_noX[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    4.0   277.0   524.0   658.6   906.0  6805.0 

sort_humans_heart_liver_DMRs_bed <- bedr.sort.region(humans_heart_liver_DMRs_bed_noX, check.chr = FALSE)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.

# Heart versus lung
humans_heart_lung_DMRs_bed <- rhesus_heart_lung_DMRs[,pull_col]

humans_heart_lung_DMRs_bed_noX <- humans_heart_lung_DMRs_bed[which(humans_heart_lung_DMRs_bed[,1] != "chr17_gl000204_random"),]

humans_heart_lung_DMRs_bed_noX[,1] <- as.character(humans_heart_lung_DMRs_bed_noX[,1])
humans_heart_lung_DMRs_bed_noX[,2] <- as.integer(humans_heart_lung_DMRs_bed_noX[,2])
humans_heart_lung_DMRs_bed_noX[,3] <- as.integer(humans_heart_lung_DMRs_bed_noX[,3])

humans_heart_lung_DMRs_bed_noX <- humans_heart_lung_DMRs_bed_noX[which(humans_heart_lung_DMRs_bed_noX[,1] != "chr11_gl000202_random"),]

humans_heart_lung_DMRs_bed_noX[,1] <- as.character(humans_heart_lung_DMRs_bed_noX[,1])
humans_heart_lung_DMRs_bed_noX[,2] <- as.integer(humans_heart_lung_DMRs_bed_noX[,2])
humans_heart_lung_DMRs_bed_noX[,3] <- as.integer(humans_heart_lung_DMRs_bed_noX[,3])



humans_heart_lung_DMRs_bed_noX <- humans_heart_lung_DMRs_bed_noX[which(humans_heart_lung_DMRs_bed_noX[,1] != "chrX"),]

humans_heart_lung_DMRs_bed_noX[,1] <- as.character(humans_heart_lung_DMRs_bed_noX[,1])
humans_heart_lung_DMRs_bed_noX[,2] <- as.integer(humans_heart_lung_DMRs_bed_noX[,2])
humans_heart_lung_DMRs_bed_noX[,3] <- as.integer(humans_heart_lung_DMRs_bed_noX[,3])
sort_humans_heart_lung_DMRs_bed <- bedr.sort.region(humans_heart_lung_DMRs_bed_noX, check.chr = FALSE)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.

summary(humans_heart_lung_DMRs_bed_noX[,3] - humans_heart_lung_DMRs_bed_noX[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    6.0   226.0   419.0   521.9   723.8  4002.0 

heart_lung <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.5 -F 0.5", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3529f0c7d.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e370fe6a48.bed -wao -f 0.5 -F 0.5

summary(heart_lung[,3] - heart_lung[,2])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    4.0   277.0   524.0   658.6   906.0  6805.0 

summary(as.integer(heart_lung[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00    0.00   35.26    0.00 3679.00 

summary(as.integer(heart_lung[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39051    2229 

lung_heart <- bedr(input = list(a = sort_humans_heart_lung_DMRs_bed, b = sort_humans_heart_liver_DMRs_bed), method = "intersect",  params = "-wao -f 0.5 -F 0.5", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35e96b10.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e363b20c8e.bed -wao -f 0.5 -F 0.5

summary(as.integer(lung_heart[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0     0.0   207.1   302.8  3679.0 

summary(as.integer(lung_heart[,9]) > 0)

   Mode   FALSE    TRUE 
logical    4797    2229 

# Make a Venn diagram based on the results

summary(as.integer(heart_lung[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39051    2229 

summary(as.integer(lung_heart[,9]) > 0)

   Mode   FALSE    TRUE 
logical    4797    2229 

grid.newpage()
draw.pairwise.venn(area1 = (39057+2229), area2 = (4798+2229), cross.area = 2229, category = c("Heart v. Liver", "Heart v. Lung"))

(polygon[GRID.polygon.29], polygon[GRID.polygon.30], polygon[GRID.polygon.31], polygon[GRID.polygon.32], text[GRID.text.33], text[GRID.text.34], text[GRID.text.35], text[GRID.text.36], text[GRID.text.37]) 

# 1 bp overlap
heart_liver_1 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.00000001 -F 0.00000001")

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36dec683b.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e378a8089c.bed -wao -f 0.00000001 -F 0.00000001

summary(as.integer(heart_liver_1[,9]) > 0)

   Mode   FALSE    TRUE 
logical   37645    3654 

# 10% overlap
heart_liver_10 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.10 -F 0.10", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e32ee745dc.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e338aec9dd.bed -wao -f 0.10 -F 0.10

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   37776    3520 

# 20% overlap
heart_liver_20 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.20 -F 0.20", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35399062a.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3140b2be4.bed -wao -f 0.20 -F 0.20

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38046    3245 

# 25% overlap
heart_liver_25 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.25 -F 0.25", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e36e168838.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e327d5f682.bed -wao -f 0.25 -F 0.25

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38205    3084 

# 30% overlap
heart_liver_30 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.30 -F 0.30", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e31b04eec5.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e35ad8399a.bed -wao -f 0.30 -F 0.30

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38353    2931 

#40% overlap
heart_liver_40 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.4 -F 0.4", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e379b6d28c.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e35f7a4e52.bed -wao -f 0.4 -F 0.4

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38703    2577 

# 50% overlap
heart_liver_50 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.50 -F 0.50", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3200e380a.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e32112e21f.bed -wao -f 0.50 -F 0.50

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39051    2229 

#60% overlap
heart_liver_60 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.6 -F 0.6", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e316151f87.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e32346f519.bed -wao -f 0.6 -F 0.6

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39421    1859 

# 70% overlap
heart_liver_70 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.70 -F 0.70", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e343596519.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e37e762ec4.bed -wao -f 0.70 -F 0.70

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39802    1478 

#75% overlap
heart_liver_75 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.75 -F 0.75", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e367436ef3.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e31470110c.bed -wao -f 0.75 -F 0.75

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39992    1288 

# 80% overlap
heart_liver_80 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.80 -F 0.80", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35d77e3c6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e320d3254a.bed -wao -f 0.80 -F 0.80

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40179    1101 

#90% overlap
heart_liver_90 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.9 -F 0.9", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e330bbc9e4.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e371a217f8.bed -wao -f 0.9 -F 0.9

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40542     738 

#100% overlap
heart_liver_100 <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 1 -F 1", check.chr = FALSE)

 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e332a59263.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e317bd74d1.bed -wao -f 1 -F 1

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40753     527 

# Find the number of tDMRs depending on the bp overlap
summary(as.integer(heart_liver_1[,9]) > 0)

   Mode   FALSE    TRUE 
logical   37645    3654 

summary(as.integer(heart_liver_10[,9]) > 0)

   Mode   FALSE    TRUE 
logical   37776    3520 

summary(as.integer(heart_liver_20[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38046    3245 

summary(as.integer(heart_liver_25[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38205    3084 

summary(as.integer(heart_liver_30[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38353    2931 

summary(as.integer(heart_liver_40[,9]) > 0)

   Mode   FALSE    TRUE 
logical   38703    2577 

summary(as.integer(heart_liver_50[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39051    2229 

summary(as.integer(heart_liver_60[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39421    1859 

summary(as.integer(heart_liver_70[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39802    1478 

summary(as.integer(heart_liver_75[,9]) > 0)

   Mode   FALSE    TRUE 
logical   39992    1288 

summary(as.integer(heart_liver_80[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40179    1101 

summary(as.integer(heart_liver_90[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40542     738 

summary(as.integer(heart_liver_100[,9]) > 0)

   Mode   FALSE    TRUE 
logical   40753     527 

# Plot the number of tDMRs
conserved_tdmr_heart_liver_heart_lung <- c(3654, 3520, 3245, 3472, 3332, 2974, 2607, 2180, 1740, 1502,1254,854, 621)
perc_overlap <- c(0.000000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung, xlab = "Percentage of overlap (bp)", ylab = "Number of overlapping rhesus tDMRs (HvLi and HvLu)",  pch = 19)

# Investigate the size of the tDMRs depending on the percentage overlap

check_1 <- heart_liver_1[which((as.integer(heart_liver_1[,9]) > 0) == TRUE),]
check_10 <- heart_liver_10[which((as.integer(heart_liver_10[,9]) > 0) == TRUE),]
check_20 <- heart_liver_20[which((as.integer(heart_liver_20[,9]) > 0) == TRUE),]
check_25 <- heart_liver_25[which((as.integer(heart_liver_25[,9]) > 0) == TRUE),]
check_30 <- heart_liver_30[which((as.integer(heart_liver_30[,9]) > 0) == TRUE),]
check_40 <- heart_liver_40[which((as.integer(heart_liver_40[,9]) > 0) == TRUE),]
check_50 <- heart_liver_50[which((as.integer(heart_liver_50[,9]) > 0) == TRUE),]
check_60 <- heart_liver_60[which((as.integer(heart_liver_60[,9]) > 0) == TRUE),]
check_70 <- heart_liver_70[which((as.integer(heart_liver_70[,9]) > 0) == TRUE),]
check_75 <- heart_liver_75[which((as.integer(heart_liver_75[,9]) > 0) == TRUE),]
check_80 <- heart_liver_80[which((as.integer(heart_liver_80[,9]) > 0) == TRUE),]
check_90 <- heart_liver_90[which((as.integer(heart_liver_90[,9]) > 0) == TRUE),]
check_100 <- heart_liver_100[which((as.integer(heart_liver_100[,9]) > 0) == TRUE),]

summary(as.integer(check_1[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.0   209.0   404.0   509.5   705.8  3679.0 

summary(as.integer(check_10[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   231.0   418.0   526.2   723.2  3679.0 

summary(as.integer(check_20[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   264.0   451.0   556.7   761.0  3679.0 

summary(as.integer(check_25[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   281.0   473.0   574.0   781.8  3679.0 

summary(as.integer(check_30[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   295.5   497.0   590.4   802.0  3679.0 

summary(as.integer(check_40[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   321.0   533.0   622.4   835.0  3679.0 

summary(as.integer(check_50[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   342.0   564.0   652.9   874.0  3679.0 

summary(as.integer(check_60[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   356.0   590.0   680.4   920.5  3679.0 

summary(as.integer(check_70[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   350.0   619.5   702.7   958.0  3679.0 

summary(as.integer(check_75[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   342.5   621.0   709.8   979.0  3679.0 

summary(as.integer(check_80[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     12     320     602     701     978    3590 

summary(as.integer(check_90[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   242.2   501.0   637.5   922.0  3590.0 

summary(as.integer(check_100[,9]))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   12.0   182.0   378.0   489.2   710.5  2266.0 

# Plot the median tDMR length
conserved_tdmr_heart_liver_heart_lung_bp <- c(404, 418, 450.5, 473, 497, 533, 564, 590, 619.5, 621,602,501, 378)
perc_overlap <- c(0.000000001, 0.1, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.80, 0.90, 1)

plot(perc_overlap, conserved_tdmr_heart_liver_heart_lung_bp, xlab = "Percentage of overlap (bp)", ylab = "Median length of overlap (bp)", pch = 19, main = "Median length of overlap between rhesus tDMRs (HvLi and HvLu)")

An overlap of 0.70-0.80 maximizes the median length of overlap between tDMRs. This result appears robust to the tissues and species.

within_species_tDMR_overlap <- function(filename1, filename2){
  
# make bed style format: change name from human to rhesus
humans_heart_liver_DMRs_bed <- filename1[,pull_col]

# Get DMRs on the autosomal chromosomes only
humans_heart_liver_DMRs_bed_noX <- humans_heart_liver_DMRs_bed[which(humans_heart_liver_DMRs_bed[,1] != "chrX"),]

humans_heart_liver_DMRs_bed_noX[,1] <- as.character(humans_heart_liver_DMRs_bed_noX[,1])
humans_heart_liver_DMRs_bed_noX[,2] <- as.integer(humans_heart_liver_DMRs_bed_noX[,2])
humans_heart_liver_DMRs_bed_noX[,3] <- as.integer(humans_heart_liver_DMRs_bed_noX[,3])

summary(humans_heart_liver_DMRs_bed_noX[,3] - humans_heart_liver_DMRs_bed_noX[,2])

sort_humans_heart_liver_DMRs_bed <- bedr.sort.region(humans_heart_liver_DMRs_bed_noX, check.chr = FALSE)

# make bed style format
humans_heart_lung_DMRs_bed <- filename2[,pull_col]
humans_heart_lung_DMRs_bed_noX <- humans_heart_lung_DMRs_bed[which(humans_heart_lung_DMRs_bed[,1] != "chrX"),]

humans_heart_lung_DMRs_bed_noX[,1] <- as.character(humans_heart_lung_DMRs_bed_noX[,1])
humans_heart_lung_DMRs_bed_noX[,2] <- as.integer(humans_heart_lung_DMRs_bed_noX[,2])
humans_heart_lung_DMRs_bed_noX[,3] <- as.integer(humans_heart_lung_DMRs_bed_noX[,3])
sort_humans_heart_lung_DMRs_bed <- bedr.sort.region(humans_heart_lung_DMRs_bed_noX, check.chr = FALSE)

summary(humans_heart_lung_DMRs_bed_noX[,3] - humans_heart_lung_DMRs_bed_noX[,2])

heart_lung <- bedr(input = list(a = sort_humans_heart_liver_DMRs_bed, b = sort_humans_heart_lung_DMRs_bed), method = "intersect",  params = "-wao -f 0.75 -F 0.75", check.chr = FALSE)

summary(heart_lung[,3] - heart_lung[,2])
summary(as.integer(heart_lung[,9]))
summary(as.integer(heart_lung[,9]) > 0)

lung_heart <- bedr(input = list(a = sort_humans_heart_lung_DMRs_bed, b = sort_humans_heart_liver_DMRs_bed), method = "intersect",  params = "-wao -f 0.75 -F 0.75", check.chr = FALSE)

summary(as.integer(lung_heart[,9]))
summary(as.integer(lung_heart[,9]) > 0)

# Make a Venn diagram based on the results

stats_filename1 <- as.array(summary(as.integer(heart_lung[,9]) > 0))
stats_filename2 <- as.array(summary(as.integer(lung_heart[,9]) > 0))

stats_filename12 <- as.data.frame(cbind(stats_filename1[2], stats_filename1[3], stats_filename1[4], stats_filename2[2], stats_filename2[3], stats_filename2[4]))

#stats_filename12[,1] <- as.integer(stats_filename12[,1])
#stats_filename12[,2] <- as.integer(stats_filename12[,2])
#stats_filename12[,3] <- as.integer(stats_filename12[,3])
#stats_filename12[,4] <- as.integer(stats_filename12[,4])
#stats_filename12[,5] <- as.integer(stats_filename12[,5])
#stats_filename12[,6] <- as.integer(stats_filename12[,6])


return(stats_filename12)

}

Heart versus liver and heart versus lung

within_species_tDMR_overlap(humans_heart_liver_DMRs, humans_heart_lung_DMRs)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e34dd967b6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e32ce3f501.bed -wao -f 0.75 -F 0.75
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e35b79bc20.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3149487ca.bed -wao -f 0.75 -F 0.75

         V1   V2   V3    V4   V5   V6
FALSE 20256 2305 <NA> 11903 2305 <NA>

grid.newpage()
draw.pairwise.venn(area1 = (20256+2305), area2 = (11903+2305), cross.area = 2305, category = c("Human Heart v. Liver", "Heart v. Lung"))

(polygon[GRID.polygon.38], polygon[GRID.polygon.39], polygon[GRID.polygon.40], polygon[GRID.polygon.41], text[GRID.text.42], text[GRID.text.43], text[GRID.text.44], text[GRID.text.45], text[GRID.text.46]) 

within_species_tDMR_overlap(chimps_heart_liver_DMRs, chimps_heart_lung_DMRs)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e3511f4114.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e35cf72fe.bed -wao -f 0.75 -F 0.75
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e346110dfb.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e32d37d33f.bed -wao -f 0.75 -F 0.75

         V1   V2   V3   V4   V5   V6
FALSE 27267 1502 <NA> 7024 1502 <NA>

grid.newpage()
draw.pairwise.venn(area1 = (27267+1502), area2 = (7024+1502), cross.area = 1502, category = c("Chimp Heart v. Liver", "Heart v. Lung"))

(polygon[GRID.polygon.47], polygon[GRID.polygon.48], polygon[GRID.polygon.49], polygon[GRID.polygon.50], text[GRID.text.51], text[GRID.text.52], text[GRID.text.53], text[GRID.text.54], text[GRID.text.55]) 

within_species_tDMR_overlap(rhesus_heart_liver_DMRs, rhesus_heart_lung_DMRs)

SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
SORTING
VALIDATE REGIONS
 * Checking input type... PASS
   Input is in bed format
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Overlapping regions can cause unexpected results.
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e32181876a.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e313de3ad8.bed -wao -f 0.75 -F 0.75
 * Processing input (1): a
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
 * Processing input (2): b
CONVERT TO BED
 * Checking input type... PASS
   Input is in bed format
VALIDATE REGIONS
 * Check if index is a string... PASS
 * Check index pattern... PASS
 * Check for missing values... PASS
 * Check for larger start position... PASS.
 * Check if zero based... PASS
 * Checking sort order... PASS
 * Checking for overlapping 'contiguous' regions... PASS
   bedtools intersect -a /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/a_c3e346542fc6.bed -b /var/folders/rf/qrcw6ncj05z1pc_pq9xzw3540000gn/T//Rtmph6wW6Y/b_c3e3492395ff.bed -wao -f 0.75 -F 0.75

         V1   V2   V3   V4   V5   V6
FALSE 39998 1288 <NA> 5739 1288 <NA>

grid.newpage()
draw.pairwise.venn(area1 = (39998+1288), area2 = (5739+1288), cross.area = 1288, category = c("Chimp Heart v. Liver", "Heart v. Lung"))

(polygon[GRID.polygon.56], polygon[GRID.polygon.57], polygon[GRID.polygon.58], polygon[GRID.polygon.59], text[GRID.text.60], text[GRID.text.61], text[GRID.text.62], lines[GRID.lines.63], text[GRID.text.64], text[GRID.text.65])