HMM call NAD
setwd('/lustre/user/liclab/pengt/projects/nucleus/bamfiles/dnabam/bedgraph/')
###########WGS
wgs1 <- fread("hela_wgs1.q30.10k.bedgraph.norm.BedGraph")
wgs2 <- fread("hela_wgs2.q30.10k.bedgraph.norm.BedGraph")
n1 <- fread("hela_n1_wgs.q30.10k.bedgraph.norm.BedGraph")
n2 <- fread("hela_n2_wgs.q30.10k.bedgraph.norm.BedGraph")
wgs <- cbind(wgs1[,c(1,2,3,6)], wgs2[,6])
wgs$mean <- rowMeans(wgs[,4:5])
nu = cbind(n1[,c(1,2,3,6)], n2[,6])
nu$mean <- rowMeans(nu[,4:5])
wgs = data.frame(wgs)
row.names(wgs) = paste0(wgs$V1,':',wgs$V2,'-',wgs$V3)
nu = data.frame(nu)
row.names(nu) = paste0(nu$V1,':',nu$V2,'-',nu$V3)
mat <- cbind(wgs[,c(1:3,6)], nu[,6])
names(mat) <- c("chromosome", "start",'end', "wgs","nucleous")
mat$ratio <- log2(mat$nucleous/mat$wgs)
mat = na.omit(mat)
mat = data.frame(mat)
mat = mat[is.finite(mat$ratio)==T,]
write.table(mat[,c(1,2,3,6)],'/lustre/user/liclab/pengt/projects/nucleus/new/NAD/wgs_all.bedgraph', row.names = F,quote = F,col.names = F,sep='\t')
###########hic
insitu2 = fread("../hicbam/bedgraph/Hela-hic-2_DO15080209_HF3TYCCXX_L4_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
insitu3 = fread("../hicbam/bedgraph/Hela-hic-3_DO15081636_HCKK5CCXX_L2_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
insitu4 = fread("../hicbam/bedgraph/Hela-hic-4_DO15081637_HCKK5CCXX_L2_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
nhic1 = fread("../hicbam/bedgraph/Hela-n1-hic_DO15081210_HF3V3CCXX_L3_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
nhic2 = fread("../hicbam/bedgraph/Hela-n2-hic_DO15081211_HF3V3CCXX_L5_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
nhic3 = fread("../hicbam/bedgraph/Hela-n3-hic_DO15081635_HCKK5CCXX_L2_clean_hg19_rDNA.q30.10k.bedgraph.norm.BedGraph")
insitu = cbind(insitu4$V6,insitu2$V6,insitu3$V6)
insitumean = rowMeans(insitu)
nhic <- cbind(nhic1$V6, nhic2$V6, nhic3$V6)
nmean <- rowMeans(nhic)
nhic <- cbind(nhic1[,1:3], nmean)
insitu = cbind(insitu2[,1:3],insitumean)
mat <- cbind(insitu[,c(1:2,3,4)], nhic[,4])
names(mat) <- c("chromosome", "start",'end', "hic","nhic")
mat$ratio <- log2(mat$nhic/mat$hic)
mat = na.omit(mat)
mat = data.frame(mat)
mat = mat[is.finite(mat$ratio)==T,]
write.table(mat[,c(1,2,3,6)],'/lustre/user/liclab/pengt/projects/nucleus/new/NAD/nhic_all.bedgraph', row.names = F,quote = F,col.names = F,sep='\t')
#########call NAD
mat <- data.frame(mat)
wgmat <- split(mat, f = mat$chromosome)
wgregions <- data.frame(chr = as.character(), start = as.numeric(), end = as.numeric())
for(i in 1:23){
if(i == 23){
chr = "chrX"
}else{
chr = paste0("chr", i)
}
chrmat <- wgmat[[chr]]
wgregions <- rbind(wgregions, getNAD(chrmat))
print(i)
}