https://github.com/HarisIqbal88/PlotNeuralNet
Showing posts with label plot. Show all posts
Showing posts with label plot. Show all posts
Wednesday, January 12, 2022
Saturday, January 19, 2019
customerized hclustfun or distfun in heatmap
for( mymethod in mymethods ) {
hd = hclust( dist(ctb2), mymethod);
# plot( hd, main="hamming distance, ward linkage" )
coat.cat = cutree(hd, numclus ) ###<=== change is here
col.palette = c("red","brown","blue","green");
coat.color = col.palette[coat.cat]
library(RColorBrewer);
#hmcol = colorRampPalette(brewer.pal(10,"RdBu"))(256);
hmcol = colorRampPalette(brewer.pal(5,"RdBu"))(16);
#heatmap( ctb2, col=hmcol, scale="none", margins = c(5,10) );
heatmap( ctb2, col=hmcol, scale="none", margins = c(5,10),
RowSideColors=coat.color, ColSideColors = spec.colors,
hclustfun = function(c) hclust( c, method=mymethod),
distfun = function(c) as.dist(hamming.distance(c)) #Hamming is less pleasant than Euclidean
main = mymethod
);
}
Tuesday, January 26, 2016
H2O2-LOH error bars, _1a.batchPlotH2O2-LOH.2016Jan25.R
# 20160125, try to add error bars to plots
rm=(list=ls())
setwd("~/github/LOH_H2O2_2016/analysis")
debug = 0;
#all data are in 'data.H2O2-LOH'
FileList = list.files( path="../data.H2O2-LOH/"); FileList;
if( debug > 5) {FileList = FileList[1:2]}
# infile = FileList[2]
for( infile in FileList) {
# infile='M1-2,20111207.modified.csv' #debug for NA, and low lead concentration
# infile = 'M2-8,08172011.H2O2.LOH.modified.csv'; #viability plot error
# infile = 'M2-8,06062011.H2O2onLOH.csv'; #b plot error
fullFileName = paste('../data.H2O2-LOH/',infile, sep='');
mylabel = infile
tb = read.csv(fullFileName, colClasses=c("character",NA, NA, "character", rep("numeric",8 ), NA));
names(tb) = c("Strain", "OD600", "Dilution","Date","H2O2stock", "White", "Black", "halfBlack", "quarterBlack", "ThreeQBlack", "QQBlack", "Other", "Notes")
######## set zeros
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other");
for( i in 1:length(tb[,1])) {
for ( j in mycolumns) {
if( is.na(tb[i,j]) ) { tb[i,j]= 0 }
}
}
#H2O2 working stock is 2X
tb$H2O2 = tb$H2O2stock/2
tb$tot = tb$White + tb$Black + tb$halfBlack + tb$quarterBlack + tb$ThreeQBlack + tb$QQBlack + tb$Other
tb.ori = tb;
tb = tb[ ! is.na(tb$White), ]
tb$Dilution = tb$Dilution / tb$Dilution[1]
######## normalize all data
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other","tot");
for ( j in mycolumns) {
tb[,j] = tb[,j] * tb$Dilution
}
####### find out means
H2O2 = sort(unique( tb$H2O2))
#s = H2O2
tbm = data.frame(cbind(H2O2))
for ( i in 1:length(H2O2)) {
c = H2O2[i]
tmp = tb[ tb$H2O2==c, ]
tbm$tot[i] = mean(tmp$tot, na.rm=T)
tbm$tot.sd[i] = sd(tmp$tot, na.rm=T)
tbm$White[i] = mean(tmp$White, na.rm=T)
tbm$White.sd[i] = sd(tmp$White, na.rm=T)
tbm$Black[i] = mean(tmp$Black, na.rm=T)
tbm$Black.sd[i] = sd(tmp$Black, na.rm=T)
tbm$halfBlack[i] = mean(tmp$halfBlack, na.rm=T)
tbm$halfBlack.sd[i] = sd(tmp$halfBlack, na.rm=T)
tbm$quarterBlack[i] = mean(tmp$quarterBlack, na.rm=T)
tbm$ThreeQBlack[i] = mean(tmp$ThreeQBlack, na.rm=T)
tbm$QQBlack[i] = mean(tmp$QQBlack, na.rm=T);
}
tbm = tbm[tbm$tot>1, ] #remove plates with zero colonies
###### some manual curations here
#tbm$halfBlack[tbm$halfBlack==0 & tbm$H2O2>0] = NA;
###### calculate fractions
tbf = tbm;
tbf$s = tbf$tot / max(tbf$tot)
tbf$s.sd = tbm$tot.sd / max(tbm$tot)
for ( j in 3:11) {
tbf[, j] = tbf[,j] / tbf$tot
}
tbf$Black[tbf$Black<0]=NA; #remove weird experimental data, such as low-lead concentration effect
pdf(paste("sandbox/",infile, ".batch.sd.pdf", sep=''), width=5, height=5)
### full black plot
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='', ylim=c(-0.05, 1.2), type='p');
lines( tbf$s ~ tbf$H2O2, col="blue",lty=2)
arrows( tbf$H2O2, (tbf$s - tbf$s.sd), tbf$H2O2, (tbf$s + tbf$s.sd), length=0.1, angle=90,code=3, lty=1, col="blue" );
axis( 4, at=pretty(c(0, 1.2)), col='black')
legend ( max(H2O2)*0.7, 0.9, c("viability","black"), col=c("blue","black"), lty=c(2,1), pch=c(1,16) )
par(new=T)
plot( tbf$Black ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="black", ylim=c(-0.01, max(tbf$Black, na.rm = T)*3))
lines(tbf$Black ~ tbf$H2O2)
bottoms = tbf$Black -tbf$Black.sd
bottoms[bottoms<0]= 0.001
arrows( tbf$H2O2, bottoms, tbf$H2O2, (tbf$Black + tbf$Black.sd), length=0.1, angle=90,code=3, lty=1, col="gray" );
title(mylabel)
## full black -log plot
#tbf$H2O2[tbf$H2O2==0] = min(H2O2[-1])/10;
#with( tbf, plot( tot ~ log10(H2O2), col="blue"));
#par(new=T)
#with( tbf, plot( Black ~ log10(H2O2), pch=16))
#tbf$H2O2[tbf$H2O2==0] = min(H2O2[-1])/10;
#with( tbf, plot( s ~ log10(H2O2), col="blue", axes=F, xlab="H2O2", ylab='viability'), ylim=c(-0.5, 1.1) );
#xlabels = log10(c(0.001,0.001, 0.005, 0.01, 0.025, 0.05, 0.1, 0.15))
#axis(1, at = xlabels, labels= 10^xlabels)
#ylabels = c(0, 0.25, 0.5, 0.75, 1.0)
#axis(2, at=ylabels, labels=ylabels);
#par(new=T)
#with( tbf, plot( Black ~ log10(H2O2), pch=16, axes=F, xlab='', ylab=''))
#axis(4, pretty(range(tbf$Black)))
#title(mylabel)
### half black plot
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='', ylim=c(0, 1.2), type='p', main=mylabel);
lines( tbf$s ~ tbf$H2O2, col="blue",lty=2)
axis( 4, at=pretty(c(0, 1.2)))
legend ( max(H2O2)*0.7, 0.5, c("viability","half-black"), col=c("blue","red"), lty=c(2,1), pch=c(1,16) )
par(new=T)
plot( tbf$halfBlack ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="half-black", col='red', ylim=c(-0.005, max(tbf$halfBlack, na.rm=T)*2 ))
lines(tbf$halfBlack ~ tbf$H2O2, col='red')
bottoms = tbf$halfBlack - tbf$halfBlack.sd
bottoms[bottoms<0]= 0
arrows( tbf$H2O2, bottoms, tbf$H2O2, (tbf$halfBlack + tbf$halfBlack.sd), length=0.1, angle=90,code=3, lty=1, col="gray" );
title(mylabel)
### half/full plot
tbf$half.vs.full = tbf$halfBlack / tbf$Black
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='', ylim=c(0, 1.2), type='p', main=mylabel);
lines( tbf$s ~ tbf$H2O2, col="blue",lty=2)
axis( 4, at=pretty(c(0, 1.2)))
legend ( max(H2O2)*0.7, 0.5, c("viability","half/full"), col=c("blue","green"), lty=c(2,1), pch=c(1,16) )
par(new=T)
plot( tbf$half.vs.full ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="half/full", col='green')
lines(tbf$half.vs.full ~ tbf$H2O2, col='green')
bottoms = tbf$Black -tbf$Black.sd
bottoms[bottoms<0]= 0.001
title(mylabel)
dev.off() #end pdf
} #infile loop
#quit("yes")
Friday, November 20, 2015
Wednesday, February 4, 2015
R plot of population heterozygoy
https://github.com/cooplab/popgen-notes/blob/master/Rcode/Loss_of_heterozyg_varying_pop.R
Tuesday, December 10, 2013
Plot H2O2-LOH figure with error bars, YPS 163, 20110413 expt
# 2013 Dec 10, use layout and mar() to do layered plots (not overlays)
#YPS163 2011April 13, plot
rm=(list=ls())
setwd("~/github/LOH_H2O2_2012-master/figures")
debug = 0;
FileList = list.files( path="../data.H2O2-LOH/", pattern="YPS163"); FileList;
filename = 'YPS163,Apr13,2011.H2O2.LOH.csv';
fullFileName = paste('../data.H2O2-LOH/',filename, sep='');
mylabel = 'YPS163 20110413'
tb = read.csv(fullFileName, colClasses=c("character",NA, NA, "character", rep("numeric",8 ), NA));
names(tb) = c("Strain", "OD600", "Dilution","Date","H2O2stock", "White", "Black", "halfBlack", "quarterBlack", "ThreeQBlack", "QQBlack", "Other", "Notes")
######## set zeros
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other");
for( i in 1:length(tb[,1])) {
for ( j in mycolumns) {
if( is.na(tb[i,j]) ) { tb[i,j]= 0 }
}
}
#adjust for low-counts
tb$Black[tb$Black<=0] = 0.5
tb$halfBlack[tb$halfBlack<=0] = 0.5
tb$Black[tb$Black<0]=NA; #remove weird experimental data, such as low-lead concentration effect
tb$H2O2 = tb$H2O2stock/2
tb$tot = tb$White + tb$Black + tb$halfBlack + tb$quarterBlack + tb$ThreeQBlack + tb$QQBlack + tb$Other
tb.ori = tb;
tb = tb[ ! is.na(tb$White), ]
tb$Dilution = tb$Dilution / tb$Dilution[1]
######## normalize all data
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other","tot");
for ( j in mycolumns) {
tb[,j] = tb[,j] * tb$Dilution
}
####### find out means, sd
H2O2 = sort(unique( tb$H2O2))
#s = H2O2
tbm = data.frame(cbind(H2O2))
for ( i in 1:length(H2O2)) {
c = H2O2[i]
tmp = tb[ tb$H2O2==c, ]
tbm$tot[i] = mean(tmp$tot, na.rm=T)
tbm$tot.sd[i] = sd(tmp$tot, na.rm=T)
tbm$White[i] = mean(tmp$White, na.rm=T)
tbm$Black[i] = mean(tmp$Black, na.rm=T)
tbm$Black.sd[i] = sd(tmp$Black, na.rm=T)
tbm$halfBlack[i] = mean(tmp$halfBlack, na.rm=T)
tbm$halfBlack.sd[i] = sd(tmp$halfBlack, na.rm=T)
tbm$quarterBlack[i] = mean(tmp$quarterBlack, na.rm=T)
tbm$ThreeQBlack[i] = mean(tmp$ThreeQBlack, na.rm=T)
tbm$QQBlack[i] = mean(tmp$QQBlack, na.rm=T);
}
tbm = tbm[tbm$tot>2, ] #remove plates with <2 colonies. At least 1 white and 1 black is needed.
###### some manual curations here
# ...
###### calculate fractions
tbf = tbm;
tbf$s = tbf$tot / max(tbf$tot)
for ( j in c("Black","halfBlack", "quarterBlack","ThreeQBlack","QQBlack")) {
tbf[, j] = tbf[,j] / tbf$tot
}
for ( j in c("tot.sd", "Black.sd","halfBlack.sd")) {
tbf[, j] = tbf[,j] / max(tbf$tot); #bug fixed, 2013 Dec 10, 18:50
tbf[, j] = ifelse(tbf[, j]<=1/700, 1/700, tbf[, j] ) #at least 0.2% error
}
head(tbf)
plot(tbf$tot.sd ~ tbf$H2O2)
#It is not clear whether color figures would be charged more or not by Aging Cell.
#### define functions
logistical.viability <- function( T, w, t ) { ret <- 1 /( 1 + ( t / T )^ w ); }
logistical.black <- function(b.max, b.min, T, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
# genome.integrity <- function(b.max, b.min, T, w, t) { 2 * (b.max - b.min) / (1 + (t/T)^w) + (1 - 2 * b.max); }
######## calclulate Cv, Cb
require(nlme)
t= tbf$H2O2; s= tbf$s; b = tbf$Black
model1 = function(v,w) { 1/( 1 + ( t / v )^ w ) }
fm.s <- gnls( s ~ model1(v,w) , start = list( v = 0.02, w=2) );
t2 = seq(0,1,by=0.001)
fit.s = logistical.viability ( fm.s$coefficients[1], fm.s$coefficients[2], t2 );
Cv = fm.s$coefficients[1]
#estimate starting values
b.tmp = sort( b[ b>0 ] );
b.min = sum(b.tmp[1])/1
b.tmp = rev(b.tmp);
b.max = sum(b.tmp[1:3])/3
modelBlack = function(v, w) { b.max - (b.max-b.min)/(1 + (t/v)^w) }
modelBlack1 = function(v) { b.max - (b.max-b.min)/(1 + (t/v)) }
#fm.b <- gnls( b ~ modelBlack1(v) , start = list( v = 0.05) );
fm.b <- gnls( b ~ modelBlack(v,w) , start = list( v = 0.005, w= 4 )); fm.b
fit.b = logistical.black(b.max, b.min, fm.b$coef[1], fm.b$coef[2], t2 )
Cb = fm.b$coef[1]
####################################################################
# do the plot
####################################################################
#### specificy the layout
#pdf(paste("_20131210", filename, "pdf", sep="."), width=6,height=9);
pdf(paste("_20131210", filename, "pdf", sep="."), width=5,height=8);
#tiff(paste("_20131210", filename, "tif", sep="."), width=240,height=500);
mat = matrix( seq(1,3), nrow=3, ncol= 1 );
layout(mat, heights= c( 1.15, 1, 1.2) );
### viability
par(mar=c(0,5,5,1))
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='Viability', ylim=c(-0.1, 1.2), type='p', pch=19);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
axis( 2, at=pretty(tbf$s), tcl=0.2, las=2 )
lines( fit.s ~ t2, col="blue",lty=1);
#legend ( max(H2O2)*0.5, 0.5, c("viability","black", "half black"), col=c("blue","black"), lty=c(2,1), pch=c(1,16) )
arrows( tbf$H2O2, (tbf$s - tbf$tot.sd), tbf$H2O2, (tbf$s + tbf$tot.sd), length=0.1, angle=90,code=3, lty=2 );
#legend ( max(H2O2)*0.5, 0.5, c("viability"), col=c("blue"), lty=c(1), pch=c(16) )
points( Cv, 0.5, pch=15, col="red", cex=1.2 );
arrows( Cv, 0.5, Cv, -1, lty=2, col="red");
mtext( "Cv",side=1,at=c(Cv*0.85), line=-1, cex=0.8 );
box()
title(mylabel)
############ full blacks
par(mar=c(0,5,0,1))
plot( tbf$Black ~ tbf$H2O2, pch=16, xlab='',ylab="Black", log='y', ylim=c(0.8E-3, 0.15), axes=F)
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
arrows( tbf$H2O2, (tbf$Black - tbf$Black.sd), tbf$H2O2, (tbf$Black + tbf$Black.sd), length=0.1, angle=90,code=3, lty=2 );
lines( fit.b ~ t2, lty=1, col='black');
points ( Cb, ( b.max/2 + b.min/2), pch=15, col="red", cex=1.2);
arrows( Cb, (b.max/2 + b.min/2), Cb, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb",side=1,at=c(Cb*1.2), line=-1, cex=0.8 );
box()
####################half black plots
#estimate starting values
b0.5 = tbf$halfBlack
# Do I need to adjust the dilution bias in [0.02,0.03]? No. Plot as is.
# b0.5[6:8] = b0.5[6:8] + 0.01
b0.5.tmp = sort( b0.5[ b0.5>0 ] );
b0.5.min = sum(b0.5.tmp[1])*1.3
b0.5.tmp = rev(b0.5.tmp);
b0.5.max = sum(b0.5.tmp[1])*1.05 #account for increasing trend
modelBlack0.5B = function(v) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^2) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5B(v) , start = list( v = 0.05) ); fm.b0.5
#fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], 2, t2 )
#logistical.black <- function(b.max, b.min, v, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
modelBlack0.5 = function(v, w) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^w) }
fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5(v,w) , start = list( v = 0.008, w=1.5) ); fm.b0.5
fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], t2 )
#logistical.blackC <- function(b.max, b.min, v, w, u, t ) { ret <- b.max - (b.max - b.min) /( u + ( t / v )^ w ); }
#modelBlack0.5C = function(v, w, u) { b0.5.max - (b0.5.max-b0.5.min)/(u + (t/v)^w) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5C(v,w,u) , start = list( v = 0.008, w=1.2, u=1.1) ); fm.b0.5
#fit.b0.5 = logistical.blackC(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], fm.b0.5$coef[3], t2 )
#plot(b0.5 ~ t )
#lines(fit.b0.5 ~ t2)
Cb0.5 = fm.b0.5$coef[1]
#par(new=T)
par(mar=c(5,5,0,1))
#plot( tbf$halfBlack ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="half-black", col='red')
plot( b0.5 ~ t, pch=16, xlab='H2O2',ylab="Half-black", col='black', log='y', ylim=c(0.8E-3,0.15), axes=F)
axis( 1, at = pretty(t), labels=T, tcl=0.3);
#axis( 2, at = c(0.001,0.01,0.05,0.1), labels=T, tcl=0.2, las=2);
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
box()
lines(fit.b0.5 ~ t2, col='black')
arrows( tbf$H2O2, (tbf$halfBlack - tbf$halfBlack.sd), tbf$H2O2, (tbf$halfBlack + tbf$halfBlack.sd),
length=0.1, angle=90,code=3, lty=2, lwd=1 );
points ( Cb0.5, ( b0.5.max/2 + b0.5.min/2), pch=15, col="red", cex=1.2);
arrows( Cb0.5, (b0.5.max/2 + b0.5.min/2), Cb0.5, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb0.5",side=1,at=c(Cb*1.2),line=-1, cex=0.8 );
dev.off()
#YPS163 2011April 13, plot
rm=(list=ls())
setwd("~/github/LOH_H2O2_2012-master/figures")
debug = 0;
FileList = list.files( path="../data.H2O2-LOH/", pattern="YPS163"); FileList;
filename = 'YPS163,Apr13,2011.H2O2.LOH.csv';
fullFileName = paste('../data.H2O2-LOH/',filename, sep='');
mylabel = 'YPS163 20110413'
tb = read.csv(fullFileName, colClasses=c("character",NA, NA, "character", rep("numeric",8 ), NA));
names(tb) = c("Strain", "OD600", "Dilution","Date","H2O2stock", "White", "Black", "halfBlack", "quarterBlack", "ThreeQBlack", "QQBlack", "Other", "Notes")
######## set zeros
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other");
for( i in 1:length(tb[,1])) {
for ( j in mycolumns) {
if( is.na(tb[i,j]) ) { tb[i,j]= 0 }
}
}
#adjust for low-counts
tb$Black[tb$Black<=0] = 0.5
tb$halfBlack[tb$halfBlack<=0] = 0.5
tb$Black[tb$Black<0]=NA; #remove weird experimental data, such as low-lead concentration effect
tb$H2O2 = tb$H2O2stock/2
tb$tot = tb$White + tb$Black + tb$halfBlack + tb$quarterBlack + tb$ThreeQBlack + tb$QQBlack + tb$Other
tb.ori = tb;
tb = tb[ ! is.na(tb$White), ]
tb$Dilution = tb$Dilution / tb$Dilution[1]
######## normalize all data
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other","tot");
for ( j in mycolumns) {
tb[,j] = tb[,j] * tb$Dilution
}
####### find out means, sd
H2O2 = sort(unique( tb$H2O2))
#s = H2O2
tbm = data.frame(cbind(H2O2))
for ( i in 1:length(H2O2)) {
c = H2O2[i]
tmp = tb[ tb$H2O2==c, ]
tbm$tot[i] = mean(tmp$tot, na.rm=T)
tbm$tot.sd[i] = sd(tmp$tot, na.rm=T)
tbm$White[i] = mean(tmp$White, na.rm=T)
tbm$Black[i] = mean(tmp$Black, na.rm=T)
tbm$Black.sd[i] = sd(tmp$Black, na.rm=T)
tbm$halfBlack[i] = mean(tmp$halfBlack, na.rm=T)
tbm$halfBlack.sd[i] = sd(tmp$halfBlack, na.rm=T)
tbm$quarterBlack[i] = mean(tmp$quarterBlack, na.rm=T)
tbm$ThreeQBlack[i] = mean(tmp$ThreeQBlack, na.rm=T)
tbm$QQBlack[i] = mean(tmp$QQBlack, na.rm=T);
}
tbm = tbm[tbm$tot>2, ] #remove plates with <2 colonies. At least 1 white and 1 black is needed.
###### some manual curations here
# ...
###### calculate fractions
tbf = tbm;
tbf$s = tbf$tot / max(tbf$tot)
for ( j in c("Black","halfBlack", "quarterBlack","ThreeQBlack","QQBlack")) {
tbf[, j] = tbf[,j] / tbf$tot
}
for ( j in c("tot.sd", "Black.sd","halfBlack.sd")) {
tbf[, j] = tbf[,j] / max(tbf$tot); #bug fixed, 2013 Dec 10, 18:50
tbf[, j] = ifelse(tbf[, j]<=1/700, 1/700, tbf[, j] ) #at least 0.2% error
}
head(tbf)
plot(tbf$tot.sd ~ tbf$H2O2)
#It is not clear whether color figures would be charged more or not by Aging Cell.
#### define functions
logistical.viability <- function( T, w, t ) { ret <- 1 /( 1 + ( t / T )^ w ); }
logistical.black <- function(b.max, b.min, T, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
# genome.integrity <- function(b.max, b.min, T, w, t) { 2 * (b.max - b.min) / (1 + (t/T)^w) + (1 - 2 * b.max); }
######## calclulate Cv, Cb
require(nlme)
t= tbf$H2O2; s= tbf$s; b = tbf$Black
model1 = function(v,w) { 1/( 1 + ( t / v )^ w ) }
fm.s <- gnls( s ~ model1(v,w) , start = list( v = 0.02, w=2) );
t2 = seq(0,1,by=0.001)
fit.s = logistical.viability ( fm.s$coefficients[1], fm.s$coefficients[2], t2 );
Cv = fm.s$coefficients[1]
#estimate starting values
b.tmp = sort( b[ b>0 ] );
b.min = sum(b.tmp[1])/1
b.tmp = rev(b.tmp);
b.max = sum(b.tmp[1:3])/3
modelBlack = function(v, w) { b.max - (b.max-b.min)/(1 + (t/v)^w) }
modelBlack1 = function(v) { b.max - (b.max-b.min)/(1 + (t/v)) }
#fm.b <- gnls( b ~ modelBlack1(v) , start = list( v = 0.05) );
fm.b <- gnls( b ~ modelBlack(v,w) , start = list( v = 0.005, w= 4 )); fm.b
fit.b = logistical.black(b.max, b.min, fm.b$coef[1], fm.b$coef[2], t2 )
Cb = fm.b$coef[1]
####################################################################
# do the plot
####################################################################
#### specificy the layout
#pdf(paste("_20131210", filename, "pdf", sep="."), width=6,height=9);
pdf(paste("_20131210", filename, "pdf", sep="."), width=5,height=8);
#tiff(paste("_20131210", filename, "tif", sep="."), width=240,height=500);
mat = matrix( seq(1,3), nrow=3, ncol= 1 );
layout(mat, heights= c( 1.15, 1, 1.2) );
### viability
par(mar=c(0,5,5,1))
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='Viability', ylim=c(-0.1, 1.2), type='p', pch=19);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
axis( 2, at=pretty(tbf$s), tcl=0.2, las=2 )
lines( fit.s ~ t2, col="blue",lty=1);
#legend ( max(H2O2)*0.5, 0.5, c("viability","black", "half black"), col=c("blue","black"), lty=c(2,1), pch=c(1,16) )
arrows( tbf$H2O2, (tbf$s - tbf$tot.sd), tbf$H2O2, (tbf$s + tbf$tot.sd), length=0.1, angle=90,code=3, lty=2 );
#legend ( max(H2O2)*0.5, 0.5, c("viability"), col=c("blue"), lty=c(1), pch=c(16) )
points( Cv, 0.5, pch=15, col="red", cex=1.2 );
arrows( Cv, 0.5, Cv, -1, lty=2, col="red");
mtext( "Cv",side=1,at=c(Cv*0.85), line=-1, cex=0.8 );
box()
title(mylabel)
############ full blacks
par(mar=c(0,5,0,1))
plot( tbf$Black ~ tbf$H2O2, pch=16, xlab='',ylab="Black", log='y', ylim=c(0.8E-3, 0.15), axes=F)
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
arrows( tbf$H2O2, (tbf$Black - tbf$Black.sd), tbf$H2O2, (tbf$Black + tbf$Black.sd), length=0.1, angle=90,code=3, lty=2 );
lines( fit.b ~ t2, lty=1, col='black');
points ( Cb, ( b.max/2 + b.min/2), pch=15, col="red", cex=1.2);
arrows( Cb, (b.max/2 + b.min/2), Cb, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb",side=1,at=c(Cb*1.2), line=-1, cex=0.8 );
box()
####################half black plots
#estimate starting values
b0.5 = tbf$halfBlack
# Do I need to adjust the dilution bias in [0.02,0.03]? No. Plot as is.
# b0.5[6:8] = b0.5[6:8] + 0.01
b0.5.tmp = sort( b0.5[ b0.5>0 ] );
b0.5.min = sum(b0.5.tmp[1])*1.3
b0.5.tmp = rev(b0.5.tmp);
b0.5.max = sum(b0.5.tmp[1])*1.05 #account for increasing trend
modelBlack0.5B = function(v) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^2) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5B(v) , start = list( v = 0.05) ); fm.b0.5
#fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], 2, t2 )
#logistical.black <- function(b.max, b.min, v, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
modelBlack0.5 = function(v, w) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^w) }
fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5(v,w) , start = list( v = 0.008, w=1.5) ); fm.b0.5
fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], t2 )
#logistical.blackC <- function(b.max, b.min, v, w, u, t ) { ret <- b.max - (b.max - b.min) /( u + ( t / v )^ w ); }
#modelBlack0.5C = function(v, w, u) { b0.5.max - (b0.5.max-b0.5.min)/(u + (t/v)^w) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5C(v,w,u) , start = list( v = 0.008, w=1.2, u=1.1) ); fm.b0.5
#fit.b0.5 = logistical.blackC(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], fm.b0.5$coef[3], t2 )
#plot(b0.5 ~ t )
#lines(fit.b0.5 ~ t2)
Cb0.5 = fm.b0.5$coef[1]
#par(new=T)
par(mar=c(5,5,0,1))
#plot( tbf$halfBlack ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="half-black", col='red')
plot( b0.5 ~ t, pch=16, xlab='H2O2',ylab="Half-black", col='black', log='y', ylim=c(0.8E-3,0.15), axes=F)
axis( 1, at = pretty(t), labels=T, tcl=0.3);
#axis( 2, at = c(0.001,0.01,0.05,0.1), labels=T, tcl=0.2, las=2);
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
box()
lines(fit.b0.5 ~ t2, col='black')
arrows( tbf$H2O2, (tbf$halfBlack - tbf$halfBlack.sd), tbf$H2O2, (tbf$halfBlack + tbf$halfBlack.sd),
length=0.1, angle=90,code=3, lty=2, lwd=1 );
points ( Cb0.5, ( b0.5.max/2 + b0.5.min/2), pch=15, col="red", cex=1.2);
arrows( Cb0.5, (b0.5.max/2 + b0.5.min/2), Cb0.5, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb0.5",side=1,at=c(Cb*1.2),line=-1, cex=0.8 );
dev.off()
Rotated axis labels in R
axis( 2, at=pretty(tbf$s), tcl=0.2, las=2 )
http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-create-rotated-axis-labels_003f
7.27 How can I create rotated axis labels?
To rotate axis labels (using base graphics), you need to use
text(), rather than mtext(), as the latter does not support par("srt").## Increase bottom margin to make room for rotated labels par(mar = c(7, 4, 4, 2) + 0.1) ## Create plot with no x axis and no x axis label plot(1 : 8, xaxt = "n", xlab = "") ## Set up x axis with tick marks alone axis(1, labels = FALSE) ## Create some text labels labels <- paste("Label", 1:8, sep = " ") ## Plot x axis labels at default tick marks text(1:8, par("usr")[3] - 0.25, srt = 45, adj = 1, labels = labels, xpd = TRUE) ## Plot x axis label at line 6 (of 7) mtext(1, text = "X Axis Label", line = 6)
When plotting the x axis labels, we use
srt = 45 for text rotation angle, adj = 1 to place the right end of text at the tick marks, and xpd = TRUE to allow for text outside the plot region. You can adjust the value of the 0.25 offset as required to move the axis labels up or down relative to the x axis. See ?par for more information.Monday, June 17, 2013
R plot in bold font
tiff("lnR-G-diploid-20130617.tif",width=480,heigh=480)
par(font=2)
plot( log(tb$a) ~ tb$b, pch=19, col='blue', xlim=c(0.05, 0.22),ylim=c(-8,-3.5),xlab="G",ylab="lnR")
text( tb$b*1.02, log(tb$a*1.1), as.character(tb$strain ) )
dev.off()
Friday, October 10, 2008
Aligned chromosome plots in R
### 101008 use my annotation, there are errors in Combimatrix's chr annotation files, such as orf19.5320
### 042308 remove hist part. remove x-axis. tick inside
### based 021908-chr-hist.R
### after background correction and non-chr5 nomalization
### histogram of expression ratio along chromosome, then histogram
### modified from 020508 version
rm(list=ls())
#left ajusted
MyBox <- function( x0,y0,h,w,col) {
xin = c( x0, x0, x0+w, x0+w );
yin = c( y0, y0+h, y0+h, y0 );
polygon(xin,yin, col=col, border=NA);
}
#centered
MyBoxC <- function( x0,y0,h,w,col) {
xin = c( x0-w/2, x0-w/2, x0+w/2, x0+w/2 );
yin = c( y0-h/2, y0+h/2, y0+h/2, y0-h/2 );
polygon(xin,yin, col=col, border=NA);
}
# Exactly form of mulinomial sampling
# vector of x and p
loglh.multinomial.sampling <- function( x, p ) {
total = sum( x );
if ( length(x) == length(p) ) {
y = p ^ x;
ret <- lfactorial(total) - sum( lfactorial(x) ) + sum(log(y) );
} else {
ret <- NA;
}
}
################################
#source ("source/021908-sor55-norm-merge.b.R");
# or
# load ("021908.bkgrdCorrected.nonChr5Normalized.RData");
load ("old/022008.bkgrdCorrected.nonChr5Normalized.RData");
out$chr.combimatrix = out$chr; #store the original chr annotation
################# the expression results
cen.orfs = c('CEN1','CENR','CEN2','CEN3','CEN4','CEN5','CEN6','CEN7');
chrs = c('Chr1','ChrR','Chr2','Chr3','Chr4','Chr5','Chr6','Chr7');
locs = c( '1', 'R','2','3','4','5','6','7');
dic= data.frame( cbind( chrs, locs) )
dic$chrs = as.character( dic$chrs );
dic$locs = as.character( dic$locs );
#tb = read.delim("cal21-annotation.csv");
tb = read.delim("cal21-annotation.csv",
colClasses = c("character","character","character","character","integer","integer","integer",
"character") );
#tb$chr = as.character( tb$chr );
cens = tb[ match(cen.orfs, tb$orf), ];
## now combine tables 101008 change
out$start = tb$start[ match( out$orf, tb$orf ) ]; # 100208
out$end = tb$end[ match( out$orf, tb$orf ) ]; # 100208
out$strand= tb$strand[ match( out$orf, tb$orf ) ]; # 100208
out$chr2 = tb$chr[ match( out$orf, tb$orf ) ]; # 100208
out$chr = dic$chrs[ match( out$chr2, dic$locs) ]
out$loc = out$start /2 + out$end/2
table(out$chr); table( out$chr2); #check, OK
save.image ("101008.bkgrdCorrected.nonChr5Normalized.myChr.RData");
xscale = 3.2E6;
xlim=c(0, xscale);
yscale = 4;
ylim = c(0.1, 6);
yy = c(0.1, 0.2, 0.5, 1, 2, 4);
xx = pretty( range( tb$end, na.rm=T ) )
cutoffs = quantile( out$fold.mut.by.wt, probs=c(0.05,0.95) );
window.bp = 0.5E5; window = 15; pcutoff = 1E-4; #
expecteds = c(0.05, 0.9, 0.05) * window;
pdf( "_101008-chrs-chisq.pdf",height=9,width=6);
#jpeg( "_42308-chrs-chisq.jpg",height=2000,width=1500, horizontal=T);
#jpeg( "_42308-chrs-chisq.jpg",height=900,width=600, horizontal=T);
mat = matrix( seq(1,length(chrs)), nrow=length( chrs), ncol= 2 );
layout(mat, heights= c( 1.15, rep(1, nrow(mat)-2), 1.2) );
#i=4; i=5
for ( i in 1:length(chrs)) {
mychr = chrs[ i ];
exp.chrtb = out[ grep( mychr, out$chr),] ; ###change 100208, not work for >10 chrs
summary(exp.chrtb) #25 NAs ??, maybe,
exp.chrtb[is.na(exp.chrtb$start),]
u = 1; #u = mean( exp.chrtb$fold.mut.by.wt); # [1] 0.7891938
if ( i == length(chrs) ) {
par(mar=c(2,2,0,1));
plot( exp.chrtb$fold.mut.by.wt ~ exp.chrtb$loc, col="white",xlab='',ylab='sor55/wt',log='y',
xlim=xlim,ylim=ylim, axes=F);
axis( 1, at = xx, tcl=0.2, line= 0);
} else {
if ( i == 1) {
par(mar=c(0,2,1,1) );
} else {
par(mar=c(0,2,0,1) );
}
plot( exp.chrtb$fold.mut.by.wt ~ exp.chrtb$loc, col="white",xlab='bp',ylab='sor55/wt',log='y',
xlim=xlim,ylim=ylim, axes=F);
axis( 1, at = xx, labels=F, tcl=0.2); #no labels, tick inside
}
box( ); #set the boxes
axis(2, at= yy, tcl=0.4, line= 0);
#axis(2, at=pretty(range(ylim)) );
# main= paste( "expression ratio sor55/wt along", mychr, sep=' ') );
for( j in 1:length(exp.chrtb$orf) ) {
yin = exp.chrtb$fold.mut.by.wt[j];
mycol = 'NA';
if (( yin > cutoffs[2] ) && ( exp.chrtb$WtMean>5)) { mycol = 'blue'; ####remove outliers??!!!!
} else if ( yin < cutoffs[1] ) { mycol = 'red';
} else { mycol= 'green'; }
exp.chrtb$type[j] = mycol;
y = ifelse ( yin > u, yin-u, yin-u );
MyBox( exp.chrtb$loc[j], u, y, xscale/1000, col= mycol);
}
MyBox( 0, u, 0.01, max(exp.chrtb$loc,na.rm=T), col="black");
points( cens$start[cens$orf==cen.orfs[i]], u+0.005 ,pch=19,cex=0.8, col="black");
#text( -5E4, 2, chrs[i] );
text( -5E3, 4, chrs[i] );
#sliding chi-sq test
sorted.locs = sort( exp.chrtb$loc );
labels = c('orf','loc', 'type', 'pchisq');
sliding.test = data.frame( matrix(nrow=length(sorted.locs), ncol=length(labels) ) );
names( sliding.test) = labels;
sliding.test$loc = sorted.locs;
sliding.test$orf = exp.chrtb$orf [match( sorted.locs, exp.chrtb$loc )];
sliding.test$type = exp.chrtb$type[match( sorted.locs, exp.chrtb$loc )];
sliding.test$clustertype = "green"; ### 101008 change
sliding.test[1:5,];
exp.chrtb[exp.chrtb$orf=='orf19.5692',] #check passed.
for( j in 1:(length(sliding.test$loc) - window + 1) ) {
xs = j:(j+window-1);
mylocs = sliding.test$loc[xs]; mylocs = mylocs[ ! is.na(mylocs) ];
if (( max(mylocs) - min(mylocs)) < window.bp ) {
types = sliding.test$type[xs];
downs = length( types[ types=="red"] );
means = length( types[ types=="green"] );
ups = length( types[ types=="blue"] );
obs = c( downs, means, ups);
obs[is.na(obs)] = 0;
#expecteds = c(0.5, 9, 0.5);
k2 = sum( (obs - expecteds)^2 / expecteds );
center = floor( j + window/2 );
sliding.test$pchisq[center] = pchisq( k2, length(obs)-1, lower.tail = F);
yplot= NA;
if ( sliding.test$pchisq[ center ] < pcutoff ) {
if ((obs[1]> expecteds[1]) && (obs[3]< expecteds[3]) ){
sliding.test$clustertype[center] = "red"; yplot = 0.15;
#rect( sliding.test$loc[j], yplot, sliding.test$loc[j+window-1], yplot+yscale/100,col="red",
# border=NA );
} else if ((obs[1] < expecteds[1]) && (obs[3] > expecteds[3]) ) {
sliding.test$clustertype[center] = "blue"; yplot = 0.15;
}
MyBoxC( sliding.test$loc[center],yplot,yscale/100,xscale/200,
col=sliding.test$clustertype[center] );
# rect( sliding.test$loc[j], yplot, sliding.test$loc[j+window-1], yplot+yscale/100,
# border=NA,sliding.test$clustertype[center] );
} else {
sliding.test$clustertype[center] = "green"; #greens are not shown
}
} ## window.bi if loop
}
#sliding.test[1:20,]
#par(new=T);
#plot ( sliding.test$pchisq ~ sliding.test$loc, col='purple', type='l',xlab='',ylab='',log='y', xlim=xlim, axes=F,lty=2);
}#chr loop
dev.off();
q("yes");
### 042308 remove hist part. remove x-axis. tick inside
### based 021908-chr-hist.R
### after background correction and non-chr5 nomalization
### histogram of expression ratio along chromosome, then histogram
### modified from 020508 version
rm(list=ls())
#left ajusted
MyBox <- function( x0,y0,h,w,col) {
xin = c( x0, x0, x0+w, x0+w );
yin = c( y0, y0+h, y0+h, y0 );
polygon(xin,yin, col=col, border=NA);
}
#centered
MyBoxC <- function( x0,y0,h,w,col) {
xin = c( x0-w/2, x0-w/2, x0+w/2, x0+w/2 );
yin = c( y0-h/2, y0+h/2, y0+h/2, y0-h/2 );
polygon(xin,yin, col=col, border=NA);
}
# Exactly form of mulinomial sampling
# vector of x and p
loglh.multinomial.sampling <- function( x, p ) {
total = sum( x );
if ( length(x) == length(p) ) {
y = p ^ x;
ret <- lfactorial(total) - sum( lfactorial(x) ) + sum(log(y) );
} else {
ret <- NA;
}
}
################################
#source ("source/021908-sor55-norm-merge.b.R");
# or
# load ("021908.bkgrdCorrected.nonChr5Normalized.RData");
load ("old/022008.bkgrdCorrected.nonChr5Normalized.RData");
out$chr.combimatrix = out$chr; #store the original chr annotation
################# the expression results
cen.orfs = c('CEN1','CENR','CEN2','CEN3','CEN4','CEN5','CEN6','CEN7');
chrs = c('Chr1','ChrR','Chr2','Chr3','Chr4','Chr5','Chr6','Chr7');
locs = c( '1', 'R','2','3','4','5','6','7');
dic= data.frame( cbind( chrs, locs) )
dic$chrs = as.character( dic$chrs );
dic$locs = as.character( dic$locs );
#tb = read.delim("cal21-annotation.csv");
tb = read.delim("cal21-annotation.csv",
colClasses = c("character","character","character","character","integer","integer","integer",
"character") );
#tb$chr = as.character( tb$chr );
cens = tb[ match(cen.orfs, tb$orf), ];
## now combine tables 101008 change
out$start = tb$start[ match( out$orf, tb$orf ) ]; # 100208
out$end = tb$end[ match( out$orf, tb$orf ) ]; # 100208
out$strand= tb$strand[ match( out$orf, tb$orf ) ]; # 100208
out$chr2 = tb$chr[ match( out$orf, tb$orf ) ]; # 100208
out$chr = dic$chrs[ match( out$chr2, dic$locs) ]
out$loc = out$start /2 + out$end/2
table(out$chr); table( out$chr2); #check, OK
save.image ("101008.bkgrdCorrected.nonChr5Normalized.myChr.RData");
xscale = 3.2E6;
xlim=c(0, xscale);
yscale = 4;
ylim = c(0.1, 6);
yy = c(0.1, 0.2, 0.5, 1, 2, 4);
xx = pretty( range( tb$end, na.rm=T ) )
cutoffs = quantile( out$fold.mut.by.wt, probs=c(0.05,0.95) );
window.bp = 0.5E5; window = 15; pcutoff = 1E-4; #
expecteds = c(0.05, 0.9, 0.05) * window;
pdf( "_101008-chrs-chisq.pdf",height=9,width=6);
#jpeg( "_42308-chrs-chisq.jpg",height=2000,width=1500, horizontal=T);
#jpeg( "_42308-chrs-chisq.jpg",height=900,width=600, horizontal=T);
mat = matrix( seq(1,length(chrs)), nrow=length( chrs), ncol= 2 );
layout(mat, heights= c( 1.15, rep(1, nrow(mat)-2), 1.2) );
#i=4; i=5
for ( i in 1:length(chrs)) {
mychr = chrs[ i ];
exp.chrtb = out[ grep( mychr, out$chr),] ; ###change 100208, not work for >10 chrs
summary(exp.chrtb) #25 NAs ??, maybe,
exp.chrtb[is.na(exp.chrtb$start),]
u = 1; #u = mean( exp.chrtb$fold.mut.by.wt); # [1] 0.7891938
if ( i == length(chrs) ) {
par(mar=c(2,2,0,1));
plot( exp.chrtb$fold.mut.by.wt ~ exp.chrtb$loc, col="white",xlab='',ylab='sor55/wt',log='y',
xlim=xlim,ylim=ylim, axes=F);
axis( 1, at = xx, tcl=0.2, line= 0);
} else {
if ( i == 1) {
par(mar=c(0,2,1,1) );
} else {
par(mar=c(0,2,0,1) );
}
plot( exp.chrtb$fold.mut.by.wt ~ exp.chrtb$loc, col="white",xlab='bp',ylab='sor55/wt',log='y',
xlim=xlim,ylim=ylim, axes=F);
axis( 1, at = xx, labels=F, tcl=0.2); #no labels, tick inside
}
box( ); #set the boxes
axis(2, at= yy, tcl=0.4, line= 0);
#axis(2, at=pretty(range(ylim)) );
# main= paste( "expression ratio sor55/wt along", mychr, sep=' ') );
for( j in 1:length(exp.chrtb$orf) ) {
yin = exp.chrtb$fold.mut.by.wt[j];
mycol = 'NA';
if (( yin > cutoffs[2] ) && ( exp.chrtb$WtMean>5)) { mycol = 'blue'; ####remove outliers??!!!!
} else if ( yin < cutoffs[1] ) { mycol = 'red';
} else { mycol= 'green'; }
exp.chrtb$type[j] = mycol;
y = ifelse ( yin > u, yin-u, yin-u );
MyBox( exp.chrtb$loc[j], u, y, xscale/1000, col= mycol);
}
MyBox( 0, u, 0.01, max(exp.chrtb$loc,na.rm=T), col="black");
points( cens$start[cens$orf==cen.orfs[i]], u+0.005 ,pch=19,cex=0.8, col="black");
#text( -5E4, 2, chrs[i] );
text( -5E3, 4, chrs[i] );
#sliding chi-sq test
sorted.locs = sort( exp.chrtb$loc );
labels = c('orf','loc', 'type', 'pchisq');
sliding.test = data.frame( matrix(nrow=length(sorted.locs), ncol=length(labels) ) );
names( sliding.test) = labels;
sliding.test$loc = sorted.locs;
sliding.test$orf = exp.chrtb$orf [match( sorted.locs, exp.chrtb$loc )];
sliding.test$type = exp.chrtb$type[match( sorted.locs, exp.chrtb$loc )];
sliding.test$clustertype = "green"; ### 101008 change
sliding.test[1:5,];
exp.chrtb[exp.chrtb$orf=='orf19.5692',] #check passed.
for( j in 1:(length(sliding.test$loc) - window + 1) ) {
xs = j:(j+window-1);
mylocs = sliding.test$loc[xs]; mylocs = mylocs[ ! is.na(mylocs) ];
if (( max(mylocs) - min(mylocs)) < window.bp ) {
types = sliding.test$type[xs];
downs = length( types[ types=="red"] );
means = length( types[ types=="green"] );
ups = length( types[ types=="blue"] );
obs = c( downs, means, ups);
obs[is.na(obs)] = 0;
#expecteds = c(0.5, 9, 0.5);
k2 = sum( (obs - expecteds)^2 / expecteds );
center = floor( j + window/2 );
sliding.test$pchisq[center] = pchisq( k2, length(obs)-1, lower.tail = F);
yplot= NA;
if ( sliding.test$pchisq[ center ] < pcutoff ) {
if ((obs[1]> expecteds[1]) && (obs[3]< expecteds[3]) ){
sliding.test$clustertype[center] = "red"; yplot = 0.15;
#rect( sliding.test$loc[j], yplot, sliding.test$loc[j+window-1], yplot+yscale/100,col="red",
# border=NA );
} else if ((obs[1] < expecteds[1]) && (obs[3] > expecteds[3]) ) {
sliding.test$clustertype[center] = "blue"; yplot = 0.15;
}
MyBoxC( sliding.test$loc[center],yplot,yscale/100,xscale/200,
col=sliding.test$clustertype[center] );
# rect( sliding.test$loc[j], yplot, sliding.test$loc[j+window-1], yplot+yscale/100,
# border=NA,sliding.test$clustertype[center] );
} else {
sliding.test$clustertype[center] = "green"; #greens are not shown
}
} ## window.bi if loop
}
#sliding.test[1:20,]
#par(new=T);
#plot ( sliding.test$pchisq ~ sliding.test$loc, col='purple', type='l',xlab='',ylab='',log='y', xlim=xlim, axes=F,lty=2);
}#chr loop
dev.off();
q("yes");
Subscribe to:
Posts (Atom)