system("scp *.vcf ~/Documents/Genewiz_hdd/Day10/SNPvcf")
system("scp *.output ~/Documents/Genewiz_hdd/Day10/CGoutput")
setwd("~/Documents/Genewiz_hdd/Day145/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 > day145-mergedSNV.vcf")
system("vcftools --relatedness --vcf day145-mergedSNV.vcf --out day145-relatedness")
system("head day145-relatedness.relatedness")
system("head ~/Documents/Genewiz_hdd/Day10/macau/day10-relate-finished.txt")
0.678322 -0.0766125 -0.0673996 -0.0233873 -0.052746 -0.0379132 -0.0721357 -0.131874 -0.0886885 -0.00624352 -0.0839034 -0.0456071
-0.0766125 0.567962 -0.0273555 -0.114651 -0.101214 -0.06135 -0.0732031 -0.0259374 -0.0922275 -0.0651632 0.0114194 -0.0946184
-0.0673996 -0.0273555 0.611564 -0.0870945 -0.0581762 -0.0896531 -0.0723983 -0.0318017 -0.0183919 -0.0388794 -0.0922383 -0.0553631
-0.0233873 -0.114651 -0.0870945 0.672552 -0.0591877 -0.0923132 -0.00924771 -0.0494841 -0.0875719 -0.0303991 -0.062276 -0.118044
-0.052746 -0.101214 -0.0581762 -0.0591877 0.662125 -0.0737998 -0.0907036 -0.0181994 -0.00401629 -0.132427 -0.0360628 -0.0532848
-0.0379132 -0.06135 -0.0896531 -0.0923132 -0.0737998 0.66284 -0.047892 -0.0806077 -0.0450114 -0.134666 -0.0893109 -0.0408424
-0.0721357 -0.0732031 -0.0723983 -0.00924771 -0.0907036 -0.047892 0.71912 -0.118116 -0.0457757 -0.0899252 -0.01862 -0.0610602
-0.131874 -0.0259374 -0.0318017 -0.0494841 -0.0181994 -0.0806077 -0.118116 0.64739 -0.093684 -0.0756179 -0.0403376 -0.0241352
-0.0886885 -0.0922275 -0.0183919 -0.0875719 -0.00401629 -0.0450114 -0.0457757 -0.093684 0.694795 -0.034276 -0.00944736 -0.0608305
-0.00624352 -0.0651632 -0.0388794 -0.0303991 -0.132427 -0.134666 -0.0899252 -0.0756179 -0.034276 0.633542 -0.0834649 -0.0292018
-0.0839034 0.0114194 -0.0922383 -0.062276 -0.0360628 -0.0893109 -0.01862 -0.0403376 -0.00944736 -0.0834649 0.754264 0.0244691
-0.0456071 -0.0946184 -0.0553631 -0.118044 -0.0532848 -0.0408424 -0.0610602 -0.0241352 -0.0608305 -0.0292018 0.0244691 0.675882
library(readr)
library(data.table)
setwd("~/Documents/Genewiz_hdd/Day145/CGoutput")
data <- list.files(path = ".", pattern = "*.output")
print(data)
db <- read_delim("~/Documents/Genewiz_hdd/Day145/CGoutput/EPI-175_CG.output", "\t", escape_double = FALSE, trim_ws = TRUE, na = ".")
db <- as.data.table(db)
db$EPI175Meth <- rowSums(cbind(db[,4], db[,7]), na.rm = TRUE)
db$EPI175TotCov <- 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/Day145/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/Day145/macau")The working directory was changed to /home/sean/Documents/Genewiz_hdd/Day145/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.meth.db2 <- meth.db[complete.cases(meth.db[,c(4:15)]),]
totcov.db2 <- totcov.db[complete.cases(totcov.db[,c(4:15)]),]
#meth.db3 <- meth.db2[apply(meth.db2[, c(4:15)], MARGIN = 1, function(x) all(x > 10)), ]
totcov.db3 <- totcov.db2[apply(totcov.db2[, c(4:15)], MARGIN = 1, function(x) all(x > 10)), ]
meth.db4 <- as.data.table(cbind(paste0(meth.db2$`#CHROM`, "-", meth.db2$POS), meth.db2[,4:ncol(meth.db2)]))
colnames(meth.db4)[1] <- "Site"
totcov.db4 <- as.data.table(cbind(paste0(totcov.db3$`#CHROM`, "-", totcov.db3$POS), totcov.db3[,4:ncol(totcov.db3)]))
colnames(totcov.db4)[1] <- "Site"
meth.db5 <- meth.db4[meth.db4$Site %in% totcov.db4$Site,]
nrow(meth.db5)[1] 4582nrow(totcov.db4)[1] 4582length(totcov.db4$Site %in% meth.db5$Site)[1] 4582length(meth.db5$Site %in% totcov.db4$Site)[1] 4582#write.table(meth.db5, file = "macau.meth.txt", sep = " ", quote = FALSE, row.names = FALSE)
#write.table(totcov.db5, file = "macau.totalcov.txt", sep = " ", quote = FALSE, row.names = FALSE)setwd("~/Documents/Genewiz_hdd/Day145/macau/")The working directory was changed to /home/sean/Documents/Genewiz_hdd/Day145/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.result_assoc <- read_delim("~/Documents/Genewiz_hdd/Day145/macau/output/Init-Pred.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(),
alpha1 = col_double(),
se_alpha1 = col_double()
)sig.results <- result_assoc[result_assoc$pvalue <= 0.05,]
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]
}
colnames(sig.results)[c(13,14)] <- c("Scaffold", "Loc")
bed.test <- as.data.frame(cbind(sig.results[,13], (sig.results[,14]), (sig.results[,14]), sig.results$beta))
colnames(bed.test)[4] <- "beta"
write.table(bed.test, file = "Init-Pred-DiffMethRegions.bed", sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)
nrow(sig.results)[1] 183setwd("~/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.result_assoc <- read_delim("~/Documents/Genewiz_hdd/Day145/macau/output/Sec-Pred.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(),
alpha1 = col_double(),
se_alpha1 = col_double()
)sig.results <- result_assoc[result_assoc$pvalue <= 0.05,]
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]
}
colnames(sig.results)[c(13,14)] <- c("Scaffold", "Loc")
bed.test <- as.data.frame(cbind(sig.results[,13], (sig.results[,14]), (sig.results[,14]), sig.results$beta))
colnames(bed.test)[4] <- "beta"
write.table(bed.test, file = "Sec-Pred-DiffMethRegions.bed", sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)
nrow(sig.results)[1] 203Add 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).