Bedtools experimentation

library(readr)
## Reads in the Macau data.
result_assoc <- read_delim("~/Downloads/result.assoc.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)

Parsed with column specification:
cols(
  id = col_character(),
  n = col_integer(),
  acpt_rate = col_double(),
  beta = col_double(),
  se_beta = col_double(),
  pvalue = col_double(),
  h = col_double(),
  se_h = col_double(),
  sigma2 = col_double(),
  se_sigma2 = col_double(),
  alpha0 = col_double(),
  se_alpha0 = col_double()
)

## Removes all regions with associated p-values greater than 0.05
sig.results <- result_assoc[result_assoc$pvalue <= 0.05,]
## Splits the Scaffold # and location # apart using the "-" character between them ad then saves it to colums 13 and 14
for(i in 1:nrow(sig.results))   {
  temp <- strsplit(sig.results$id[i], "-") 
  sig.results[i,13] <- temp[[1]][1]
  sig.results[i,14] <- temp[[1]][2]
  
}
temp <- rep(10, length.out = nrow(sig.results))
## Names the two columns
colnames(sig.results)[c(13,14)] <- c("Scaffold", "Loc")
## Pulls out the scaffold name, location, and h values associated with all significant DMRs
bed.test <- as.data.frame(cbind(sig.results[,13], (sig.results[,14]), (sig.results[,14]), sig.results$h))
colnames(bed.test)[4] <- "h"
write.table(bed.test, file = "bed-test.txt", sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)
## insert some excel wizardry here to remove the row label column and column labels so genomeCoverageBed wouldn't get mad
## FAI indexes the draft genome
system("samtools faidx Panopea_generosa-Scaff-10k.fa.txt")
## Trims out the two columns we care about from the FAI indexed file
system("cut -f 1,2 Panopea_generosa-Scaff-10k.fa.txt.fai > chrom.sizes")
## Turns the Macau DMR results in to a bedgraph
system("genomeCoverageBed -bg -trackline -i bed-test.txt -g chrom.sizes > test.bedgraph")

So the above works… mostly. For some reason the genomeCoverageBed command does not transver over the data column (of h values) to the bedgraph. To fix this, I cheated and just pulled the header from the generated bedgraph file and appended it to the original bed-text.txt file, and renamed it to bed-text.bedgraph. IGV seemed happy, though I’m still working through in my head what that actually means. It’s progress though, I think?