Day 10 Samples

setwd("~/Documents/Genewiz_hdd/Day145")

The working directory was changed to /home/sean/Documents/Genewiz_hdd/Day145 inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.

system("mkdir SNPvcf")
system("mkdir CGoutput")
samples.vec <- c("EPI-175", "EPI-176", "EPI-181", "EPI-182", "EPI-184", "EPI-185", "EPI-187", "EPI-188", "EPI-193", "EPI-195", "EPI-199", "EPI-200", "EPI-205", "EPI-206", "EPI-208", "EPI-209", "EPI-214", "EPI-215", "EPI-220", "EPI-221", "EPI-226", "EPI-227", "EPI-229", "EPI-230")

system("scp *.vcf ~/Documents/Genewiz_hdd/Day10/SNPvcf")
system("scp *.output ~/Documents/Genewiz_hdd/Day10/CGoutput")

setwd("~/Documents/Genewiz_hdd/Day10/SNPvcf")

for(i in 1:length(samples.vec))   {
  
  system(paste0("bgzip ", samples.vec[i], "_SNVs.vcf"))
  system(paste0("tabix ", samples.vec[i], "_SNVs.vcf.gz"))
  
}

system("vcf-merge *.gz > day10-mergedSNV.vcf")

system("vcftools --relatedness --vcf day10-mergedSNV.vcf --out day10-relatedness")

system("head day10-relatedness.relatedness")

system("head ~/Documents/Genewiz_hdd/Day10/macau/day10-relate-finished.txt")

library(readr)
setwd("~/Documents/Genewiz_hdd/Day10/CGoutput")

data <- list.files(path = ".", pattern = "*.output")
print(data)

db <- read_delim("~/Documents/Genewiz_hdd/Day10/CGoutput/EPI-103_CG.output", "\t", escape_double = FALSE, trim_ws = TRUE, na = ".")

db <- as.data.table(db)

db$EPI103Meth <- rowSums(cbind(db[,4], db[,7]), na.rm = TRUE)
db$EPI103TotCov <- rowSums(cbind(db[,5], db[,8]), na.rm = TRUE)
meth.db <- db[,c(1,2,3,10)]
totcov.db <- db[,c(1,2,3,11)]

for(i in 2:length(data))   {
  
  temp <- read_delim(paste0("~/Documents/Genewiz_hdd/Day10/CGoutput/", data[i]), "\t", escape_double = FALSE, trim_ws = TRUE, na = ".")
  temp <- as.data.table(temp)
  temp$Meth <- rowSums(cbind(temp[,4], temp[,7]), na.rm = TRUE)
  temp$TotCov <- rowSums(cbind(temp[,5], temp[,8]), na.rm = TRUE)

  meth.db <- merge(meth.db, temp[,c(1,2,3,10)], by = names(meth.db[,c(1,2,3)]), all = TRUE)
  colnames(meth.db)[ncol(meth.db)] <- paste0(substr(data[i], 1, 7), "Meth")
  meth.db <- as.data.table(meth.db)
  
  totcov.db <- merge(totcov.db, temp[,c(1,2,3,11)], by = names(totcov.db[,c(1,2,3)]), all = TRUE)
  colnames(totcov.db)[ncol(totcov.db)] <- paste0(substr(data[i], 1, 7), "TotCov")
  totcov.db <- as.data.table(totcov.db)
}

setwd("~/Documents/Genewiz_hdd/Day10/macau")

The working directory was changed to /home/sean/Documents/Genewiz_hdd/Day10/macau inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.

write.table(meth.db5, file = "macau.meth2.txt", sep = " ", quote = FALSE, row.names = FALSE)
write.table(totcov.db4, file = "macau.totalcov2.txt", sep = " ", quote = FALSE, row.names = FALSE)

setwd("~/Documents/Genewiz_hdd/Day10/macau")

The working directory was changed to /home/sean/Documents/Genewiz_hdd/Day10/macau inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the the working directory for notebook chunks.

system("macau -g macau.meth2.txt -t macau.totalcov2.txt -p day10-pred.csv -k day10-relate-finished.csv -o output3 -s 10000 -seed 123")

Reading Files ... 
## number of total individuals = 12
## number of analyzed individuals = 12
## number of covariates = 1
## number of total genes/sites = 15558
Performing Analysis                                                   0.00%
Performing Analysis ===                                               6.43%
Performing Analysis ======                                            12.86%
Performing Analysis =========                                         19.28%
Performing Analysis ============                                      25.71%
Performing Analysis ================                                  32.14%
Performing Analysis ===================                               38.57%
Performing Analysis ======================                            45.00%
Performing Analysis =========================                         51.42%
Performing Analysis ============================                      57.85%
Performing Analysis ================================                  64.28%
Performing Analysis ===================================               70.71%
Performing Analysis ======================================            77.14%
Performing Analysis =========================================         83.56%
Performing Analysis ============================================      89.99%
Performing Analysis ================================================  96.42%
Performing Analysis ==================================================100.00%

nrow(sig.results)

[1] 267

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).