R/replicationdomains.R

In dvera/arthur: What the package does (short line)

replicationdomains <-
function(bgfile ){


	ttrname<-paste0(basename(removeext(bgfile)),"_ttr.bed")
	rdname<-paste0(basename(removeext(bgfile)),"_rds.bed")



	cat(basename(bgfile),": loading bedGraph\n")
	r<-tsvRead(bgfile)
	rl<-nrow(r)

	cat(basename(bgfile),": finding TTRs\n")

	#find sign-changes
	r$prod=c(r[,4][1:(rl-1)]*r[,4][2:rl],1)

	#calculate slope at each point
	slopes=diff(r[,4])/diff(r[,2])
	r$slope=c(slopes,slopes[rl-1])
	rm(slopes)

	#find slop inversions
	r$slopeprod=c(r$slope[1],r$slope[1:(rl-1)]*r$slope[2:rl])
	slopeinvs<-which(r$slopeprod<0)
	signchanges<-which(r$prod<0)

	cat(basename(bgfile),": finding left TTR boundaries\n")
	leftttrb<-unlist(lapply(1:length(signchanges), function(x) signchanges[x]-which(r$slopeprod[signchanges[x]:1]<0)[1]+1 ))
	cat(basename(bgfile),": finding right TTR boundaries\n")
	ritettrb<-unlist(lapply(1:length(signchanges), function(x) signchanges[x]+which(r$slopeprod[(signchanges[x]):(signchanges[length(signchanges)])]<0)[1]-1 ))

	#t<-data.frame(signchanges,c(NA,leftttrb),c(ritettrb,NA))
	t<-data.frame(signchanges,leftttrb,ritettrb)




	cat(basename(bgfile),": removing incomplete TTRs at chromosome ends\n")
	cbreaks1 <- which(  r[t[,2],2] > r[t[,1],2]  )
	if(length(cbreaks1)>0){ t[cbreaks1,2] <- NA }

	cbreaks2 <- which(  r[t[,3],2] < r[t[,1],2]  )
	if(length(cbreaks2)>0){ t[cbreaks2,3] <- NA }

	# t=t[complete.cases(t),]
	# t=t[which(t[,2]!=t[,3]),]

	cat(basename(bgfile),": finding TTR directionality\n")
	t$strand <- NA
	t[which(r[signchanges,4]<0),4]<-"-"
	t[which(r[signchanges,4]>0),4]<-"+"

	cat(basename(bgfile),": saving TTRs to ",ttrname,"\n")
	ttr<-data.frame("V1"=r[t[,1],1],"V2"=r[t[,2],2],"V3"=r[t[,3],2],"V4"=1:nrow(t),"V5"=1,"V6"=t[,4],"V7"=r[t[,1],2],"V8"=r[t[,1],2]+1)
	ttr<-ttr[which(complete.cases(ttr)),]
	write.tsv(ttr,file=ttrname)

	cat(basename(bgfile),": saving RDs to ",rdname,"\n")
	rd<-data.frame("V1"=r[t[1:(nrow(t)-1),1],1],"V2"=r[t[1:(nrow(t)-1),1],2],"V3"=r[t[2:(nrow(t)),1],2], "V4"=1:(nrow(t)-1) , "V5"=1 , "V6"=t[2:(nrow(t)),4] , "V7"=r[t[1:(nrow(t)-1),3],2],"V8"=r[t[2:(nrow(t)),2],3])
	rd<-rd[which(rd[,3]>rd[,2]),]
	rd<-rd[which(complete.cases(rd)),]
	write.tsv(rd,file=rdname)
	return(c(ttrname,rdname)
	#chromends<-prerd[which(prerd[,3]<prerd[,2]),]
	#chrombegs<-chromends
	#chrombegs[,2]<-1
	#chromends[,3]<-gf[match(chromends$V1,gf$V1),2]
		#otherchrom=r[t[2:(nrow(t)),1],1]
	#rd<-rd[which(rd$V1==otherchrom),]
	#rm(otherchrom)

}