data=read.csv('Spike-in_counts.csv')
extract=subset(data,data$Strategy=='Grep_full-length')
wt_range=range(pretty(c(0,max(extract$wt_spike_count))))
m2_range=range(pretty(c(0,max(extract$m2_spike_count))))
m6_range=range(pretty(c(0,max(extract$m6_spike_count))))

pdf('Spike-in_reproducibility.pdf',width=11,height=6)
par(mfrow=c(1,2))
cor_coef=signif(as.numeric(cor.test(extract$wt_spike_count,extract$m2_spike_count,method='pearson')$estimate),digits=3)
pval=signif(cor.test(extract$wt_spike_count,extract$m2_spike_count,method='pearson')$p.value,digits=3)
plot(extract$wt_spike_count,extract$m2_spike_count,xlim=wt_range,ylim=m2_range,xlab='Number of wt spike-in reads',ylab='Number of m2 spike-in reads',main=paste("Pearson's correlation coefficient: ",cor_coef,"\n(p-value=",pval,")",sep=''),axes=F)
axis(1);axis(2)
cor_coef=signif(as.numeric(cor.test(extract$wt_spike_count,extract$m6_spike_count,method='pearson')$estimate),digits=3)
pval=signif(cor.test(extract$wt_spike_count,extract$m6_spike_count,method='pearson')$p.value,digits=3)
plot(extract$wt_spike_count,extract$m6_spike_count,xlim=wt_range,ylim=m6_range,xlab='Number of wt spike-in reads',ylab='Number of m6 spike-in reads',main=paste("Pearson's correlation coefficient: ",cor_coef,"\n(p-value=",pval,")",sep=''),axes=F)
axis(1);axis(2)
dev.off()
