Count the conserved bases and conserved region within a genomic region from phastCons44way UCSC track

conservedbases.R

cmd_args <- commandArgs(TRUE)

chr <- cmd_args[1]
cbeg <- cmd_args[2]
cend <- cmd_args[3]

tmp <- read.table("tmp",header=F)
nbases <- length(which(tmp$V1 != 0))
vec <- tmp$V1
i <- 1
region <- 0
while(i <= length(vec)-1){
if(vec[i]==0 & vec[i+1]!=0){
 region <- region +1}
i<-i+1
}

line <- paste(chr,cbeg,cend,nbases,region,sep="\t")
write(line,"coverage.txt",append=T)