args=commandArgs(trailingOnly=TRUE) genome=args[1] genome_length=args[2] x=c(1:genome_length) all=x for (sample in c('embryon_8h','embryon_15h','embryon_36h','embryon_60h','femelle','male')) for (lib in 1:4) { input_file=paste(genome,'_mappers_in_',sample,'_',lib,'.dat',sep='') stat=read.table('/mnt/data/home/herve.seitz/Amphioxus/Small_RNA_mapping/Statistics.dat',header=T) depth=(stat$Genome.matching_reads-stat$ncRNA.matching_reads)[stat$Sample==sample & stat$Library==lib] data=read.table(input_file,header=T) length=c() y_sense=rep(0,times=length(x)) y_antisense=rep(0,times=length(x)) if (length(data$Sequence)>0) { for (i in 1:length(data$Sequence)) length=append(length,nchar(as.character(data$Sequence[i]))) starts=data$Position ends=starts+length-1 for (i in 1:length(data$Sequence)) for (bp in starts[i]:ends[i]) if (data$Orientation[i]==0) { y_sense[bp]=y_sense[bp]+1 } else y_antisense[bp]=y_antisense[bp]+1 y_sense=y_sense*1e6/depth y_antisense=-y_antisense*1e6/depth } all=cbind(all,y_sense[x],y_antisense[x]) } y_range=range(pretty(c(min(all[,2*c(1:24)+1]),max(all[,2*c(1:24)])))) colors=c('red','orange','green','cyan','blue','purple') pdf(paste(genome,'_coverage.pdf',sep=''),width=6,height=6) plot(x,y_sense,ty='n',ylim=y_range,xlab=paste('Coordinate in ',genome,' genome (bp)',sep=''),ylab='Read coverage (ppm)') for (i in 1:24) { j=(i+3)%/%4 k=(i+3)%%4+1 par(new=T) plot(x,all[,2*i],ty='l',lwd=2,col=colors[j],ylim=y_range,xlab='',ylab='',axes=F,lty=k) par(new=T) plot(x,all[,2*i+1],ty='l',lwd=2,col=colors[j],ylim=y_range,xlab='',ylab='',axes=F,lty=k) } #legend('topleft',c('8 hpf embryo','15hpf embryo','36 hpf embryo','60 hpf embryo','adult female','adult male'),pch='-',col=colors) #legend('topright',paste('lib.',1:4),lty=c(1:4),col='black') dev.off()