calculate_regional_disorder.R

## pre flight: break up reference genome into 200 bp windows using bedtools & use bedtools map to get average per read level PDN values for each 200 bp window 

##handling the per genomic window levels of PDN - each of the 200 bp windows of the genome has the average PDN value of the reads covering it (reads have to be 51% or more within the region) and
##the number of reads which cover that region

library(dplyr)
library(stringr)
library(readr)


##set unique file ID's
unique_fileID <- c("SRR13012807",
"SRR13012806",
"SRR13012822",
"SRR13012821",
"SRR13012820",
"SRR13012819",
"SRR13012818",
"SRR13012817",
"SRR13012816",
"SRR13012815",
"SRR13012814",
"SRR13012813",
"SRR13012812",
"SRR13012811",
"SRR13012810",
"SRR13012809",
"SRR13012808",
"SRR13012838",
"SRR13012837",
"SRR13012836",
"SRR13012835",
"SRR13012834",
"SRR13012833",
"SRR13012832",
"SRR13012831",
"SRR13012830",
"SRR13012829",
"SRR13012828",
"SRR13012827",
"SRR13012826",
"SRR13012825",
"SRR13012824",
"SRR13012823",
"SRR13012855",
"SRR13012854",
"SRR13012853",
"SRR13012852",
"SRR13012851",
"SRR13012850",
"SRR13012849",
"SRR13012848",
"SRR13012847",
"SRR13012846",
"SRR13012845",
"SRR13012844",
"SRR13012843",
"SRR13012842",
"SRR13012841",
"SRR13012840",
"SRR13012839",
"SRR13012866",
"SRR13012865",
"SRR13012864",
"SRR13012863",
"SRR13012862",
"SRR13012861",
"SRR13012860",
"SRR13012859",
"SRR13012858",
"SRR13012857",
"SRR13012856",
"SRR13012876",
"SRR13012875",
"SRR13012874",
"SRR13012873",
"SRR13012872",
"SRR13012871",
"SRR13012870",
"SRR13012869",
"SRR13012868",
"SRR13012867",
"SRR13012891",
"SRR13012890",
"SRR13012889",
"SRR13012888",
"SRR13012887",
"SRR13012886",
"SRR13012885",
"SRR13012884",
"SRR13012883",
"SRR13012882",
"SRR13012881",
"SRR13012880",
"SRR13012879",
"SRR13012878",
"SRR13012877",
"SRR13012909",
"SRR13012908",
"SRR13012907",
"SRR13012906",
"SRR13012905",
"SRR13012904",
"SRR13012903",
"SRR13012902",
"SRR13012901",
"SRR13012900",
"SRR13012899",
"SRR13012898",
"SRR13012897",
"SRR13012896",
"SRR13012895",
"SRR13012894",
"SRR13012893",
"SRR13012892",
"SRR13012924",
"SRR13012923",
"SRR13012922",
"SRR13012921",
"SRR13012920",
"SRR13012919",
"SRR13012918",
"SRR13012917",
"SRR13012916",
"SRR13012915",
"SRR13012914",
"SRR13012913",
"SRR13012912",
"SRR13012911",
"SRR13012910",
"SRR13012938",
"SRR13012937",
"SRR13012936",
"SRR13012935",
"SRR13012934",
"SRR13012933",
"SRR13012932",
"SRR13012931",
"SRR13012930",
"SRR13012929",
"SRR13012928",
"SRR13012927",
"SRR13012926",
"SRR13012925",
"SRR13012805")

##set file paths
filePaths <- paste0("/scratch/emb19132/GENOME_PDR/RAT/PDN_output/", unique_fileID, "_PDN_overlap.txt")

##loop to calculate PDN
for(i in 1:length(unique_fileID)){
  ##read in overlap file written as a plain text file with headers removed
  overlap <- read.table(filePaths[i], header = FALSE, fill = TRUE)
  ##replace any regions with less than 5 reads covering it (number of reads is V5)
  overlap2 <- overlap %>%  mutate_at(vars(V4, V5), ~replace(., V5 <= 4, NA)) ##less than or equal to 4
  ##rename the avg PDN with sample name
  assign(paste(unique_fileID[i], "PDN", sep = "_"), overlap2$V4)
  assign(paste(unique_fileID[i], "reads", sep = "_"), overlap2$V5)
  rm(overlap) ##remove large files from memory to reduce load
  rm(overlap2)
  }

##now bind all sample PDN values together with original coordinate values
coord <- read.table("/scratch/emb19132/GENOME_PDR/RAT/PDN_output/SRR13012807_PDN_overlap.txt",  header = FALSE)
##give each set of coordinates a name (chrom - V1, start - V2, end - V3)
coord$region.name <- paste(coord$V1,coord$V2,coord$V3, sep=":")
assign(paste("genomic_region"), coord$region.name)

compiled <- cbind(genomic_region, SRR13012807_PDN,
SRR13012806_PDN,
SRR13012822_PDN,
SRR13012821_PDN,
SRR13012820_PDN,
SRR13012819_PDN,
SRR13012818_PDN,
SRR13012817_PDN,
SRR13012816_PDN,
SRR13012815_PDN,
SRR13012814_PDN,
SRR13012813_PDN,
SRR13012812_PDN,
SRR13012811_PDN,
SRR13012810_PDN,
SRR13012809_PDN,
SRR13012808_PDN,
SRR13012838_PDN,
SRR13012837_PDN,
SRR13012836_PDN,
SRR13012835_PDN,
SRR13012834_PDN,
SRR13012833_PDN,
SRR13012832_PDN,
SRR13012831_PDN,
SRR13012830_PDN,
SRR13012829_PDN,
SRR13012828_PDN,
SRR13012827_PDN,
SRR13012826_PDN,
SRR13012825_PDN,
SRR13012824_PDN,
SRR13012823_PDN,
SRR13012855_PDN,
SRR13012854_PDN,
SRR13012853_PDN,
SRR13012852_PDN,
SRR13012851_PDN,
SRR13012850_PDN,
SRR13012849_PDN,
SRR13012848_PDN,
SRR13012847_PDN,
SRR13012846_PDN,
SRR13012845_PDN,
SRR13012844_PDN,
SRR13012843_PDN,
SRR13012842_PDN,
SRR13012841_PDN,
SRR13012840_PDN,
SRR13012839_PDN,
SRR13012866_PDN,
SRR13012865_PDN,
SRR13012864_PDN,
SRR13012863_PDN,
SRR13012862_PDN,
SRR13012861_PDN,
SRR13012860_PDN,
SRR13012859_PDN,
SRR13012858_PDN,
SRR13012857_PDN,
SRR13012856_PDN,
SRR13012876_PDN,
SRR13012875_PDN,
SRR13012874_PDN,
SRR13012873_PDN,
SRR13012872_PDN,
SRR13012871_PDN,
SRR13012870_PDN,
SRR13012869_PDN,
SRR13012868_PDN,
SRR13012867_PDN,
SRR13012891_PDN,
SRR13012890_PDN,
SRR13012889_PDN,
SRR13012888_PDN,
SRR13012887_PDN,
SRR13012886_PDN,
SRR13012885_PDN,
SRR13012884_PDN,
SRR13012883_PDN,
SRR13012882_PDN,
SRR13012881_PDN,
SRR13012880_PDN,
SRR13012879_PDN,
SRR13012878_PDN,
SRR13012877_PDN,
SRR13012909_PDN,
SRR13012908_PDN,
SRR13012907_PDN,
SRR13012906_PDN,
SRR13012905_PDN,
SRR13012904_PDN,
SRR13012903_PDN,
SRR13012902_PDN,
SRR13012901_PDN,
SRR13012900_PDN,
SRR13012899_PDN,
SRR13012898_PDN,
SRR13012897_PDN,
SRR13012896_PDN,
SRR13012895_PDN,
SRR13012894_PDN,
SRR13012893_PDN,
SRR13012892_PDN,
SRR13012924_PDN,
SRR13012923_PDN,
SRR13012922_PDN,
SRR13012921_PDN,
SRR13012920_PDN,
SRR13012919_PDN,
SRR13012918_PDN,
SRR13012917_PDN,
SRR13012916_PDN,
SRR13012915_PDN,
SRR13012914_PDN,
SRR13012913_PDN,
SRR13012912_PDN,
SRR13012911_PDN,
SRR13012910_PDN,
SRR13012938_PDN,
SRR13012937_PDN,
SRR13012936_PDN,
SRR13012935_PDN,
SRR13012934_PDN,
SRR13012933_PDN,
SRR13012932_PDN,
SRR13012931_PDN,
SRR13012930_PDN,
SRR13012929_PDN,
SRR13012928_PDN,
SRR13012927_PDN,
SRR13012926_PDN,
SRR13012925_PDN,
SRR13012805_PDN)

##remove rows with more than 20% NAs
compiled <- as.data.frame(compiled)
row.names(compiled) <- compiled$genomic_region
compiled <- compiled[,c(-1)]
compiled <- compiled[!rowSums(is.na(compiled)) > ncol(compiled)*.2,] ##remove regions not covered by at least 20% of samples

write.csv(compiled, file="/scratch/emb19132/GENOME_PDR/RAT/PDN_output/rat_merged_region_PDN.csv")