rendered paste body
R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.1 (2006-06-01)
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>
>
>
> ###
> ### Set up output
> ###
>
> ## Use R defaults!
> postscript(file = "Image/plot.3x.ps")
>
> ## Give left margin a bit more space
> #par( mar = c( 5, 4.5, 4, 2 ) + 0.1 )
>
> ## Fix margins (as appropriate)
> par( mex = 1.5 )
>
>
>
>
>
>
> ##
> ## Clear the decks
> ##
>
> ls()
character(0)
> rm(list=ls())
>
>
>
>
>
> ###
> ### Read in datasets...
> ###
>
> ##
> ## dat.or = original data set of oligomers (miller et al. 1987a)
> ## dat.ho = my collected homo-oligomers and subunit group data
> ## dat.he = my collected hetero-oligomers and subunit group data
> ##
>
> dat.or <-
+ read.table( "miller1987.table2.pisa.fix.dat", header=TRUE)
> dat.ho <-
+ read.table( "my.table2.pisa.ho.csv", header=TRUE, sep='\t')
> dat.he <-
+ read.table( "my.table2.pisa.he.csv", header=TRUE, sep='\t')
> dat.he.sub <-
+ read.table( "my.table2.pisa.he.grup.csv", header=TRUE, sep='\t')
>
> ## Test
> head(dat.or, 2)
So Name Id Type SGo Mo Ao Bbo Ao.Bbo
1 2 Avian pancreatic polypeptide 1ppt Ho 1 8493.38 5334.5 1363.3 3971.2
2 2 Subtilisin inhibitor 3ssi Ho 1 23006.20 10147.6 1724.3 8423.3
Ms_1 Ms_2 As Bbs As.Bbs
1 4246.69 4216 3348.9 681.7 2667.2
2 11503.10 11015 5935.9 862.1 5073.8
> head(dat.ho, 2)
So Id.1 Id.2 Mo Ao Bbo Ms As Bbs
1 2 1mxr 1 78802 23670.2 6061.6 39401 14865.9 3030.8
2 2 1h0o 1 67412 23186.2 3107.8 33706 13147.0 1553.9
Name
1 Ribonucleotide reductase R2 from E. coli
2 Ribonucleotide reductase R2 from mouse
> head(dat.he, 2)
Id.1 Id.2 Mo So SGo Ao Bbo
1 1cs0 1 157224 2 2 47427.9 4734.6
2 1dce 1 101420 2 2 37541.1 5598.0
Name
1 Carbamoyl phosphate synthetase
2 Rab geranylgeranyltransferase from rat brain
> head(dat.he.sub, 2)
Id.1 Id.2 Mo So SGo Ao Bbo Name chGroup
1 1cs0 1 157224 2 2 47427.8 4734.7 . A
2 1cs0 1 157224 2 2 47427.8 4734.7 Carbamoyl phosphate synthetase B
Ms Sg As Bbs
1 116188 1 36977.4 2367.3
2 41036 1 15185.1 2367.3
>
> nrow(dat.or)
[1] 21
> nrow(dat.ho)
[1] 169
> nrow(dat.he)
[1] 88
> nrow(dat.he.sub)
[1] 249
>
>
> ## Fix data...
>
> ## Fix col names
> dat.ho$Id <- dat.ho$Id.1
> dat.he$Id <- dat.he$Id.1
>
> dat.he.sub$Id <- dat.he.sub$Id.1
> dat.he.sub$chG <- dat.he.sub$chGroup
>
> ## Hacks
> dat.ho$SGo <- 1
> dat.ho$chG <- "A"
> dat.or$chG <- "A"
>
> ## Very important!
> dat.or$Ms <- dat.or$Ms_2
> dat.or$Mo <- dat.or$Ms * dat.or$So
>
>
>
> ## More tests
> length(unique(dat.or$Id))
[1] 21
> length(unique(dat.ho$Id))
[1] 169
> length(unique(dat.he$Id))
[1] 88
> length(unique(dat.he.sub$Id))
[1] 88
>
>
>
>
>
> ## Prepare the overall OLIGOMER data (with several fields)
> my.fields.o = c( "Id", "So", "SGo", "Mo", "Ao", "Bbo", "Name" )
> #
> dat.to <-
+ data.frame(rbind(cbind( Type = "Or", dat.or[ , my.fields.o ] ),
+ cbind( Type = "Ho", dat.ho[ , my.fields.o ] ),
+ cbind( Type = "He", dat.he[ , my.fields.o ] ) ) )
>
> ## TEST
> head(dat.to, 3)
Type Id So SGo Mo Ao Bbo Name
1 Or 1ppt 2 1 8432 5334.5 1363.3 Avian pancreatic polypeptide
2 Or 3ssi 2 1 22030 10147.6 1724.3 Subtilisin inhibitor
3 Or 2ccy 2 1 26520 12692.3 1622.7 Cytochrome c
> nrow(dat.to); 169 + 88 + 21
[1] 278
[1] 278
>
>
>
> ## Prepare the overall SUBUNIT data (with several fields)
> my.fields.s = c( "Id", "So", "SGo", "chG", "Ms", "As", "Bbs", "Name" )
> ##
> dat.so <-
+ data.frame(rbind(cbind( Type = "Or", dat.or [ , my.fields.s ] ),
+ cbind( Type = "Ho", dat.ho [ , my.fields.s ] ),
+ cbind( Type = "He", dat.he.sub[ , my.fields.s ] ) ) )
>
> ## TEST
> head(dat.so, 3)
Type Id So SGo chG Ms As Bbs Name
1 Or 1ppt 2 1 A 4216 3348.9 681.7 Avian pancreatic polypeptide
2 Or 3ssi 2 1 A 11015 5935.9 862.1 Subtilisin inhibitor
3 Or 2ccy 2 1 A 13260 7158.0 811.4 Cytochrome c
> nrow(dat.so); 169 + 249 + 21
[1] 439
[1] 439
>
>
>
> ## PREPARE
> ho <- which( dat.so$Type == "Ho" )
> he <- which( dat.so$Type == "He" )
> or <- which( dat.so$Type == "Or" )
>
>
>
>
>
> ###
> ### Trend (overview.. Why is this here and not elsewhere?)
> ###
>
> ##
> ## Fit *new* data.Find approximate parameters using the log linear
> ## relationship,
> ##
>
> ##
> ## y = k1 * x ^ k2
> ## log(y) = log(k1) + log(x) * k2
> ##
> ## y = c + x * m
> ##
>
>
> ## Linear regression of log data
> lm.or <- lm( log(Ao) ~ log(Mo), data=dat.or ); summary( lm.or )
Call:
lm(formula = log(Ao) ~ log(Mo), data = dat.or)
Residuals:
Min 1Q Median 3Q Max
-0.10633 -0.03546 -0.02833 0.04692 0.12424
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.48191 0.18006 8.23 1.10e-07 ***
log(Mo) 0.77756 0.01616 48.10 < 2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 0.06557 on 19 degrees of freedom
Multiple R-Squared: 0.9919, Adjusted R-squared: 0.9914
F-statistic: 2314 on 1 and 19 DF, p-value: < 2.2e-16
> lm.ho <- lm( log(Ao) ~ log(Mo), data=dat.ho ); summary( lm.ho )
Call:
lm(formula = log(Ao) ~ log(Mo), data = dat.ho)
Residuals:
Min 1Q Median 3Q Max
-0.364594 -0.063337 -0.002938 0.083001 0.265661
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.10956 0.12964 0.845 0.399
log(Mo) 0.89834 0.01075 83.582 <2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 0.1151 on 167 degrees of freedom
Multiple R-Squared: 0.9767, Adjusted R-squared: 0.9765
F-statistic: 6986 on 1 and 167 DF, p-value: < 2.2e-16
> lm.he <- lm( log(Ao) ~ log(Mo), data=dat.he ); summary( lm.he )
Call:
lm(formula = log(Ao) ~ log(Mo), data = dat.he)
Residuals:
Min 1Q Median 3Q Max
-0.38396 -0.07423 0.01846 0.08754 0.31445
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.55025 0.16891 3.258 0.00161 **
log(Mo) 0.86447 0.01437 60.139 < 2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 0.1293 on 86 degrees of freedom
Multiple R-Squared: 0.9768, Adjusted R-squared: 0.9765
F-statistic: 3617 on 1 and 86 DF, p-value: < 2.2e-16
> ##
> lm.to <- lm(log(Ao) ~ log(Mo), data=dat.to); summary(lm.to)
Call:
lm(formula = log(Ao) ~ log(Mo), data = dat.to)
Residuals:
Min 1Q Median 3Q Max
-0.400863 -0.066219 0.002338 0.084392 0.321525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.423775 0.094645 4.478 1.11e-05 ***
log(Mo) 0.873221 0.007956 109.760 < 2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 0.1209 on 276 degrees of freedom
Multiple R-Squared: 0.9776, Adjusted R-squared: 0.9775
F-statistic: 1.205e+04 on 1 and 276 DF, p-value: < 2.2e-16
>
> ##
> ## Grab parameters
> ##
>
> ## Parameters of slope
> or.k1 <- exp(coef(lm.or)[1]) ## Unlog k1
> or.k2 <- coef(lm.or)[2]
> ho.k1 <- exp(coef(lm.ho)[1]) ## Unlog k1
> ho.k2 <- coef(lm.ho)[2]
> he.k1 <- exp(coef(lm.he)[1]) ## Unlog k1
> he.k2 <- coef(lm.he)[2]
> ##
> to.k1 <- exp(coef(lm.to)[1]) ## Unlog k1
> to.k2 <- coef(lm.to)[2]
>
> ## SD of parameters of slope (from summary)
> or.k1.pm <- exp(summary(lm.or)$coefficients[1,2])
> or.k2.pm <- summary(lm.or)$coefficients[2,2]
> ho.k1.pm <- exp(summary(lm.ho)$coefficients[1,2])
> ho.k2.pm <- summary(lm.ho)$coefficients[2,2]
> he.k1.pm <- exp(summary(lm.he)$coefficients[1,2])
> he.k2.pm <- summary(lm.he)$coefficients[2,2]
> ##
> to.k1.pm <- exp(summary(lm.to)$coefficients[1,2])
> to.k2.pm <- summary(lm.to)$coefficients[2,2]
>
> ## Grab R squared
> or.r <- summary(lm.or)$r.squared
> ho.r <- summary(lm.ho)$r.squared
> he.r <- summary(lm.he)$r.squared
> ##
> to.r <- summary(lm.to)$r.squared
>
>
> ## Print parameters
> my.sigFig = 4
> ##
> my.param <-
+ data.frame(Type=c("Or","Ho","He","To"),
+ k1=c(
+ round(or.k1, my.sigFig),
+ round(ho.k1, my.sigFig),
+ round(he.k1, my.sigFig),
+ ##
+ round(to.k1, my.sigFig)),
+ k1.pm=c(
+ round(or.k1.pm, my.sigFig),
+ round(ho.k1.pm, my.sigFig),
+ round(he.k1.pm, my.sigFig),
+ ##
+ round(to.k1.pm, my.sigFig)),
+ k2=c(
+ round(or.k2, my.sigFig),
+ round(ho.k2, my.sigFig),
+ round(he.k2, my.sigFig),
+ ##
+ round(to.k2, my.sigFig)),
+ k2.pm=c(
+ round(or.k2.pm, my.sigFig),
+ round(ho.k2.pm, my.sigFig),
+ round(he.k2.pm, my.sigFig),
+ ##
+ round(to.k2.pm, my.sigFig)),
+ R2=c(
+ round(or.r, my.sigFig),
+ round(ho.r, my.sigFig),
+ round(he.r, my.sigFig),
+ ##
+ round(to.r, my.sigFig))
+ )
> ##
> my.param
Type k1 k1.pm k2 k2.pm R2
1 Or 4.4014 1.1973 0.7776 0.0162 0.9919
2 Ho 1.1158 1.1384 0.8983 0.0107 0.9767
3 He 1.7337 1.1840 0.8645 0.0144 0.9768
4 To 1.5277 1.0993 0.8732 0.0080 0.9776
>
>
>
>
>
>
>
> ###
> ### Calculate
> ###
>
> ##
> ## At = k3 + Mo * k4 (DEFINED ELSEWHERE)
> ##
>
> or.k3 <- -6.2362829
> or.k4 <- 1.5625581
>
>
>
> ## Calculate Ato, Bto and Bwo (OLIGOMER)
> ##
> dat.to$Ato <- or.k3 + dat.to$Mo * or.k4
> ##
> dat.to$Bto <- dat.to$Ato - dat.to$Ao
> dat.to$Bwo <- dat.to$Bto - dat.to$Bbo
>
> ## Calculate Ats, Bts and Bws (SUBUNITS)
> ##
> dat.so$Ats <- or.k3 + dat.so$Ms * or.k4
> ##
> dat.so$Bws <- dat.so$Ats - dat.so$As
> dat.so$Bts <- dat.so$Bws + dat.so$Bbs
>
>
>
> ## FRAC atack?
>
>
>
> ##
> ## Sample data
> ##
>
> ## Sample some data (OLIGOMER)
> head(cbind( dat.to[ , c( "Type", "Id", "So" ) ],
+ round( dat.to[ , c( "Mo", "Ato", "Ao" ) ] ),
+ round( dat.to[ , c( "Bto", "Bbo", "Bwo" ) ] ) ), 2 )
Type Id So Mo Ato Ao Bto Bbo Bwo
1 Or 1ppt 2 8432 13169 5334 7835 1363 6471
2 Or 3ssi 2 22030 34417 10148 24269 1724 22545
>
> ## Sample some data (SUBUNIT)
> head(cbind( dat.so[ , c( "Id", "So" ) ],
+ round( dat.so[ , c( "Ms", "Ats", "As" ) ] ),
+ round( dat.so[ , c( "Bts", "Bbs", "Bws" ) ] ) ), 2 )
Id So Ms Ats As Bts Bbs Bws
1 1ppt 2 4216 6582 3349 3914 682 3233
2 3ssi 2 11015 17205 5936 12132 862 11269
>
>
>
> ## Checks to hecks?
>
>
>
>
>
> ###
> ### SUMMARY OF BURIED AREA (PER-SUBUNIT)
> ###
>
> ## BIG ONES
> tail(dat.so[order(dat.so$Bbs),
+ c("Type", "Id", "So", "SGo", "chG","As", "Bbs")], 12)
Type Id So SGo chG As Bbs
75 Ho 1jkv 6 1 A 16775.67 10044.80
231 He 1e9y 24 2 B 22251.50 10223.40
27 Ho 1gte 2 1 A 45692.55 11334.90
20 Or 7cat 4 1 A 25789.20 11367.50
1231 He 1hbn 6 3 AD 23331.00 12356.00
197 He 2occ 13 13 A 19678.20 12845.10
1161 He 1iw7 6 6 C 57742.10 13231.00
28 Ho 1kek 2 1 A 51376.20 13505.60
233 He 1ocr 26 13 AN 19671.50 13686.85
1151 He 1iw7 6 6 D 74165.00 16344.80
166 Ho 1c8d 60 1 A 31285.00 16992.50
1212 He 1k28 6 2 A 39351.80 20285.00
>
> ## SIZES (RANGE)
> round(
+ sapply( tapply( dat.so$Bbs, dat.so$Type, range ), c)
+ )
Or Ho He
[1,] 682 916 616
[2,] 11368 16992 20285
>
> ## SIZES (SUMMARY)
> mean(dat.so$Bbs); sd(dat.so$Bbs); median(dat.so$Bbs)
[1] 3353.671
[1] 2489.772
[1] 2708
>
> ## SIZES (SUMMARY - PER GROUP)
> t(cbind(
+ number = tapply( dat.so$Bbs, dat.so$Type, length),
+ ##
+ mean = round( tapply( dat.so$Bbs, dat.so$Type, mean) ),
+ sd = round( tapply( dat.so$Bbs, dat.so$Type, sd) ),
+ median = round( tapply( dat.so$Bbs, dat.so$Type, median) ) ) )
Or Ho He
number 21 169 249
mean 2647 3599 3247
sd 2357 2229 2653
median 1828 3010 2519
>
>
>
> ##
> ## PLOTS
> ##
>
> my.y <- "Interface size per subunit"
> ##
> boxplot( Bbs ~ Type, ylab=my.y, data = dat.so )
> boxplot( Bbs ~ Type, ylab=my.y, data = dat.so, varwidth = T )
> boxplot( Bbs ~ Type, ylab=my.y, data = dat.so, varwidth = T, notch = T )
> ##
> abline( h = mean( log( dat.so$Bbs ) ), col = "grey" )
> abline( h = median( log( dat.so$Bbs ) ), col = "grey" )
>
>
> my.y <- "log ( Interface size per subunit )"
> ##
> boxplot( log( Bbs ) ~ Type, ylab=my.y, data = dat.so )
> boxplot( log( Bbs ) ~ Type, ylab=my.y, data = dat.so, varwidth = T )
> boxplot( log( Bbs ) ~ Type, ylab=my.y, data = dat.so, varwidth = T, notch = T )
> ##
> abline( h = mean( log( dat.so$Bbs ) ), col = "grey" )
> abline( h = median( log( dat.so$Bbs ) ), col = "grey" )
>
>
> plot( sort( dat.or$Bbs ) )
> plot( sort( dat.so$Bbs[ho] ) )
> plot( sort( dat.so$Bbs[he] ) )
>
> plot( sort( dat.or$Bbs ), log = "y" )
> plot( sort( dat.so$Bbs[ho] ), log = "y" )
> plot( sort( dat.so$Bbs[he] ), log = "y" )
>
> plot( ecdf( dat.or$Bbs ) )
> plot( ecdf( dat.so$Bbs[ho] ), add = T )
> plot( ecdf( dat.so$Bbs[he] ), add = T )
>
>
>
> ## HISTO-GRAMMY
>
> ## ONE
> my.mids <- hist(dat.so$Bbs)
> my.brk <- my.mids$breaks
> my.x <- length(my.brk) - 1
>
> abline(v= mean(dat.so$Bbs), col="grey")
> abline(v=median(dat.so$Bbs), col="grey")
>
> ( my.a <- table( c( cut( dat.so$Bbs[ho], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10 11
31 88 32 12 2 2 1 0 1 0 0
> ( my.b <- table( c( cut( dat.so$Bbs[he], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10 11
81 114 28 12 7 1 4 0 1 0 1
>
> barplot( rbind(Ho=my.a,He=my.b), beside=T, legend=T )
>
>
> ## TWO
> my.mids <- hist(dat.so$Bbs[dat.so$Bbs<10000])
> my.brk <- my.mids$breaks
> my.x <- length(my.brk) - 1
>
> abline(v= mean(dat.so$Bbs), col="grey")
> abline(v=median(dat.so$Bbs), col="grey")
>
> ( my.a <- table( c( cut( dat.so$Bbs[ho], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10
2 29 53 35 20 12 5 7 2 0
> ( my.b <- table( c( cut( dat.so$Bbs[he], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10
17 64 74 40 16 12 9 3 4 3
>
> barplot( rbind(Ho=my.a,He=my.b), beside=T, legend=T )
>
>
> ## TRE
> my.mids <- hist(dat.so$Bbs[dat.so$Bbs<5000])
> my.brk <- my.mids$breaks
> my.x <- length(my.brk) - 1
>
> abline(v= mean(dat.so$Bbs), col="grey")
> abline(v=median(dat.so$Bbs), col="grey")
>
> ( my.a <- table( c( cut( dat.so$Bbs[ho], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9
2 16 13 24 29 21 14 12 8
> ( my.b <- table( c( cut( dat.so$Bbs[he], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9
17 36 28 42 32 25 15 6 10
>
> barplot( rbind(Ho=my.a,He=my.b), beside=T, legend=T )
>
>
> ## FOU
> hist(log(dat.so$Bbs),prob=T,breaks=20)
>
> abline(v= mean(log(dat.so$Bbs)), col="grey")
> abline(v=median(log(dat.so$Bbs)), col="grey")
>
> curve( dnorm(x,
+ mean = mean( log( dat.so$Bbs ) ),
+ sd = sd( log( dat.so$Bbs ) ) ), col=2, lwd=2, add=T )
>
>
>
>
>
> ## PLOTTIE FOR SCOTTIE
> my.mids <- hist(dat.so$Bbs)
> my.brk <- my.mids$breaks
> my.x <- length(my.brk) - 1
>
> abline(v= mean(dat.so$Bbs), col="grey")
> abline(v=median(dat.so$Bbs), col="grey")
>
> ( my.a <- table( c( cut( dat.so$Bbs[ho], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10 11
31 88 32 12 2 2 1 0 1 0 0
> ( my.b <- table( c( cut( dat.so$Bbs[he], my.brk, label=F ), 1:my.x ) ) - 1 )
1 2 3 4 5 6 7 8 9 10 11
81 114 28 12 7 1 4 0 1 0 1
>
> barplot( rbind(Ho=my.a,He=my.b), beside=T,
+ xlab="Interface size",
+ ylab="No. Subunits",
+ cex.axis=1.3,
+ cex.lab=1.5,
+ xaxt="n")
>
> my.ticks <- c( 0, 5000, 10000, 15000, 20000 )
>
> my.ticks <- my.ticks/(20000/29)
>
> ##
> axis(side=1,at=my.ticks+1,
+ label=c( 0, 5000, 10000, 15000, 20000 ))
>
>
>
>
>
>
>
> ##
> ## Histogram of the interface area (buried surface area) per subunit
> ##
>
> range( dat.so$Bbs )
[1] 616.5 20285.0
> round( range( dat.so$Bbs ) / 10 ) * 10
[1] 620 20280
>
> hist(dat.so$Bbs,
+ main="",
+ xlab=expression("Interface size / " * ring(A) * phantom()^2 ),
+ ylab="No. Subunits",
+ cex.axis=1.3,
+ cex.lab=1.5
+ )
> ##
> abline(v=median(dat.so$Bbs), lty=2, col=2)
> ##
> split.screen( matrix( c( 0.4, 0.95, 0.4, 0.95 ), ncol=4 ), erase = FALSE )
[1] 1
> ##
> par( mar = c( 2, 2, 0, 2 ) + 0.1 )
> ##
> hist(dat.so$Bbs[dat.so$Bbs<5000],
+ main="",
+ xlim=c(0,5000),
+ cex.axis=1.3,
+ cex.lab=1.5
+ )
> ##
> abline(h=0)
> abline(v=median(dat.so$Bbs), lty=2, col=2)
> ##
> close.screen( all = TRUE )
>
>
> ##q()
>
>
>
>
>
>
>
> ## CHECK INTERFACE AREA AGAINST MOLECULAR WEIGHT
>
> plot(dat.so$Ms,
+ dat.so$Bbs)
>
> summary( my.lm <- lm( Bbs ~ Ms, data = dat.so) )
Call:
lm(formula = Bbs ~ Ms, data = dat.so)
Residuals:
Min 1Q Median 3Q Max
-6377.8 -1158.5 -328.8 756.3 15315.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.718e+03 1.717e+02 10.01 <2e-16 ***
Ms 5.307e-02 4.448e-03 11.93 <2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 2165 on 437 degrees of freedom
Multiple R-Squared: 0.2457, Adjusted R-squared: 0.2439
F-statistic: 142.3 on 1 and 437 DF, p-value: < 2.2e-16
>
> abline( coef( my.lm ), col=2 )
>
>
> plot(dat.so$Ms,
+ dat.so$Bbs, log="xy")
>
> summary( my.lm <- lm( log( Bbs, 10 ) ~ log( Ms, 10 ),
+ data = dat.so) )
Call:
lm(formula = log(Bbs, 10) ~ log(Ms, 10), data = dat.so)
Residuals:
Min 1Q Median 3Q Max
-0.60974 -0.15248 0.01690 0.15550 0.71437
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.80570 0.15001 12.04 <2e-16 ***
log(Ms, 10) 0.37330 0.03421 10.91 <2e-16 ***
---
Signif. codes: 0 $-1òø***òù 0.001 òø**òù 0.01 òø*òù 0.05 òø.òù 0.1 òø òù 1
Residual standard error: 0.2393 on 437 degrees of freedom
Multiple R-Squared: 0.2141, Adjusted R-squared: 0.2123
F-statistic: 119.1 on 1 and 437 DF, p-value: < 2.2e-16
>
> abline( coef( my.lm ), col=2 )
>
>
> ## Multiplot
> my.plot.f <- function( x, lty = 1, ... ) {
+ ##
+ my.i <-
+ dat.so$So == x &
+ dat.so$Type == "Ho"
+ ##
+ points(dat.so$Ms[my.i],
+ dat.so$Bbs[my.i], ... )
+ ##
+ abline( coef( lm( Bbs ~ Ms,
+ data = dat.so[my.i,]) ),
+ lty=lty, ... )
+ }
>
> plot(dat.so$Ms,
+ dat.so$Bbs, type="n")
>
> table(dat.so$So[dat.so$Type=="Ho"])
2 3 4 5 6 7 8 10 11 12 14 16 24 60 180
29 11 23 11 34 3 20 9 2 10 4 1 7 4 1
>
> my.plot.f( 2 )
> my.plot.f( 3, pch = 2, col = 2 )
> my.plot.f( 4, pch = 3, col = 3 )
> my.plot.f( 5, pch = 4, col = 4 )
> my.plot.f( 6, pch = 5, col = 5 )
> my.plot.f( 8, pch = 6, col = 6 )
> my.plot.f( 10, pch = 7, col = 7 )
> my.plot.f( 12, pch = 8, col = 8 )
> my.plot.f( 24, pch = 9, col = 9 )
>
>
>
>
>
>
>
> ## CHECK INTERFACE AREA AGAINST OLIGOMER CLASS
>
> my.plot.f <- function( v, ... ) {
+ ##
+ boxplot( Bbs ~ So,
+ ylab="Interface area (per subunit)",
+ xlab="Subunits per oligomer",
+ data = v, xaxt="n", ... )
+ ##
+ n.out <- sort( unique( v$So ) )
+ ##
+ axis( side = 1, at=c( 1:length( n.out ) ),
+ label = n.out, las = 2 )
+ ##
+ }
>
> ## OK
>
> my.plot.f( dat.so )
> my.plot.f( dat.so, varwidth=T )
> my.plot.f( dat.so, varwidth=T, notch=T )
>
> my.plot.f( dat.so[ dat.so$So <= 12, ] )
> my.plot.f( dat.so[ dat.so$So <= 12, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 12, ], varwidth=T, notch=T )
>
> my.plot.f( dat.so[ dat.so$So <= 8, ] )
> my.plot.f( dat.so[ dat.so$So <= 8, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 8, ], varwidth=T, notch=T )
>
> my.plot.f( dat.so[ dat.so$So <= 6, ] )
> my.plot.f( dat.so[ dat.so$So <= 6, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 6, ], varwidth=T, notch=T )
>
>
>
>
>
>
>
> ## CHECK *FRACTION* INTERFACE AREA AGAINST OLIGOMER CLASS
>
> my.plot.f <- function( v, ... ) {
+ ##
+ ##v <- v[dat.so$Type=="Ho",]
+ ##
+ boxplot( Bbs/As ~ So,
+ ylab = "Fraction interface area (per subunit)",
+ data = v, xaxt="n", ... )
+ ##
+ n.out <- sort( unique( v$So ) )
+ ##
+ axis( side = 1, at=c( 1:length( n.out ) ),
+ label = n.out, las = 2 )
+ ##
+ }
>
> my.f <- function( n ) {
+ 1 - n^( 0.8 - 1 )
+ }
>
> ## OK
>
> my.plot.f( dat.so )
> my.plot.f( dat.so, varwidth=T )
> my.plot.f( dat.so, varwidth=T, notch=T )
>
> #points(x = c(1:15),
> # y = my.f( c(2,3,4,5,6,7,8,10,11,12,14,16,24,60,180) ),
> # pch=3, cex=3, col=3 )
>
> points(x = c(1:24),
+ y = my.f( c(2,3,4,5,6,7,8,9,10,11,12,13,14,
+ 16,18,20,21,22,24,26,28,60,180,360) ),
+ pch=3, cex=3, col=3 )
>
>
> my.plot.f( dat.so[ dat.so$So <= 12, ] )
> my.plot.f( dat.so[ dat.so$So <= 12, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 12, ], varwidth=T, notch=T )
>
> points(x = c(1,2,3,4,5,6,7,8, 9,10,11),
+ y = my.f( c(2,3,4,5,6,7,8,9,10,11,12) ), pch=3, cex=3, col=4 )
>
>
> my.plot.f( dat.so[ dat.so$So <= 8, ] )
> my.plot.f( dat.so[ dat.so$So <= 8, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 8, ], varwidth=T, notch=T )
>
> points(x = c(1,2,3,4,5,6,7),
+ y = my.f( c(2,3,4,5,6,7,8) ), pch=3, cex=3, col=3 )
>
> my.plot.f( dat.so[ dat.so$So <= 6, ] )
> my.plot.f( dat.so[ dat.so$So <= 6, ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So <= 6, ], varwidth=T, notch=T )
>
> points(x = c( 1, 2, 3, 4, 5 ),
+ y = my.f(c(2,3,4,5,6)), pch=3, cex=3, col=4 )
>
> my.plot.f( dat.so[ dat.so$So %in% c(2,3,4,6), ] )
> my.plot.f( dat.so[ dat.so$So %in% c(2,3,4,6), ], varwidth=T )
> my.plot.f( dat.so[ dat.so$So %in% c(2,3,4,6), ], varwidth=T, notch=T )
>
>
>
>
>
>
>
>
> ## PLOT from Ponstingl 2005
>
> my.y <- "Mw of subunits"
>
> boxplot( Ms ~ So, ylab = my.y, data = dat.so )
>
> boxplot( Ms ~ So, ylab = my.y, data = dat.so[dat.so$So <= 12, ] )
> boxplot( Ms ~ So, ylab = my.y, data = dat.so[dat.so$So <= 12, ],
+ varwidth=T )
> boxplot( Ms ~ So, ylab = my.y, data = dat.so[dat.so$So <= 12, ],
+ varwidth=T, notch=T )
>
>
>
>
>
>
>
>
>
> ###
> ### RESIDUALS!!
> ###
>
> to.k1 <- 1.440669
> to.k2 <- 0.879367
>
>
> ## SURFACE AREA (Ao / As)
> dat.to$Ao.p <- to.k1 * dat.to$Mo ^ to.k2
> dat.so$As.p <- to.k1 * dat.so$Ms ^ to.k2
>
> dat.to$Ao.d <- (dat.to$Ao - dat.to$Ao.p) / dat.to$Ao * 100
> dat.so$As.d <- (dat.so$As - dat.so$As.p) / dat.so$As * 100
>
>
> ## BURIED AREA (Bto / Bts)
> dat.to$Bto.p <- dat.to$Ato - dat.to$Ao.p
> dat.so$Bws.p <- dat.so$Ats - dat.so$As.p
>
> dat.to$Bto.d <- (dat.to$Bto - dat.to$Bto.p) / dat.to$Bto * 100
> dat.so$Bws.d <- (dat.so$Bws - dat.so$Bws.p) / dat.so$Bws * 100
>
>
> ## TEST
> mean(dat.to$Ao.d ); mean(abs(dat.to$Ao.d )); fivenum(dat.to$Ao.d )
[1] -2.187613
[1] 9.661773
[1] -50.345076 -8.278252 -1.357469 6.742604 25.645793
> mean(dat.so$As.d ); mean(abs(dat.so$As.d )); fivenum(dat.so$As.d )
[1] 12.46715
[1] 14.92970
[1] -39.654632 1.991580 11.091581 20.111239 54.294523
>
> mean(dat.to$Bto.d); mean(abs(dat.to$Bto.d)); fivenum(dat.to$Bto.d)
[1] 0.0579248
[1] 2.632342
[1] -15.9834611 -2.2174843 0.3383551 2.1419333 10.2631599
> mean(dat.so$Bws.d); mean(abs(dat.so$Bws.d)); fivenum(dat.so$Bws.d)
[1] -11.04762
[1] 11.80701
[1] -372.3731065 -11.3796739 -4.9385843 -0.7172051 12.5076619
>
>
> ## Double check
> plot( dat.to$Ao.d, dat.to$Bto.d );
> abline(h=0,col="grey"); abline(v=0,col="grey")
>
> plot(dat.so$As.d, dat.so$Bws.d);
> abline(h=0,col="grey"); abline(v=0,col="grey")
>
> ## Huh?
> dat.so[ ( x1 <- dat.so$Bws.d < -100 ), ]
Type Id So SGo chG Ms As Bbs Name Ats Bws Bts
1141 He 1ao7 5 5 C 1071 1455.1 1332.6 . 1667.263 212.1634 1544.763
As.p As.d Bws.p Bws.d
1141 665.0604 54.29452 1002.203 -372.3731
> dat.so[ ( x2 <- dat.so$Bws.d > 10 ), ]
Type Id So SGo chG Ms As Bbs Name Ats Bws Bts
249 He 2bbv 360 2 DEF 4461 1669.933 726.6 . 6964.335 5294.402 6021.002
As.p As.d Bws.p Bws.d
249 2332.139 -39.65463 4632.196 12.50766
>
> ## OK
> plot(dat.so$As.d [ ! x1 & ! x2 ], # Remember, x1 and x2
+ dat.so$Bws.d[ ! x1 & ! x2 ]) # are logical vectors.
>
> abline(h=0,col="grey"); abline(v=0,col="grey")
>
>
>
> ## REMOVE
> dat.so <- dat.so[ ! x1 & ! x2, ] #
>
>
> ## BUT DON'T FORGET...
> ho <- which( dat.so$Type == "Ho" )
> he <- which( dat.so$Type == "He" )
> or <- which( dat.so$Type == "Or" )
>
>
>
> ## Try again...
> matplot(x=cbind(dat.to$Ao.d, dat.so$As.d),
+ y=cbind(dat.to$Bto.d, dat.so$Bws.d),
+ pch=1,col=c(1,2));
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 1) in: cbind(1, dat.to$Ao.d, dat.so$As.d)
2: number of rows of result
is not a multiple of vector length (arg 1) in: cbind(1, dat.to$Bto.d, dat.so$Bws.d)
> abline(h=0,col="grey"); abline(v=0,col="grey")
>
> matplot(cbind(dat.so$As.d, dat.to$Ao.d),
+ cbind(dat.so$Bws.d, dat.to$Bto.d),
+ pch=1,col=c(2,1))
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$As.d, dat.to$Ao.d)
2: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Bws.d, dat.to$Bto.d)
> abline(h=0,col="grey"); abline(v=0,col="grey")
>
>
>
>
>
>
>
>
>
>
>
>
>
> ###
> ### Plot
> ###
>
> ## Below we plot the histograms. First we make a plot for Cyrus.
>
> ## PLOT 1
> plot ( dat.so$Ms, dat.so$Bws, ylim=range(dat.so$Bws.p) )
> points( dat.so$Ms, dat.so$Bws.p, col = 2 )
>
> ## PLOT 2
> plot ( dat.so$Ms, dat.so$Bws, ylim=range(dat.so$Bws.p) )
> curve ( (or.k3 + x * or.k4) -
+ (to.k1 * x ^ to.k2), add=T, col = 2, lwd = 2 )
>
> ## PLOT 3
> plot ( dat.so$Ms, dat.so$Bws.d )
> ##
> abline(h=0, col="darkgrey")
>
>
> ## PLOT 4
> matplot(cbind(dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or]),
+ cbind(dat.so$Bws.d[he], dat.so$Bws.d[ho], dat.so$Bws.d[or]),
+ pch=c(3,1,4),
+ col=c(2,1,3),
+ lwd=c(1,1,2)
+ )
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or])
2: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Bws.d[he], dat.so$Bws.d[ho], dat.so$Bws.d[or])
> ##
> abline(h=0, col="darkgrey")
>
> legend(x = "bottomright", inset = 0.05, bg = "white",
+ c("Subunits from the homo-oligomers",
+ "Subunits from the hetero-oligomers",
+ "Subunits from the original (21) oligomrs" ),
+ pch=c(3,1,4),
+ col=c(2,1,3),
+ pt.lwd=c(1,1,2)
+ )
>
>
>
>
> ## PLOT 3.L
> plot ( dat.so$Ms, dat.so$Bws.d, log="x" )
> ##
> abline(h=0, col="darkgrey")
>
>
> ## PLOT 4.L
> matplot(cbind(dat.so$Ms[ho], dat.so$Ms[he], dat.so$Ms[or]),
+ cbind(dat.so$Bws.d[ho], dat.so$Bws.d[he], dat.so$Bws.d[or]),
+ pch=c(1,3,4), log="x"
+ )
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 1) in: cbind(1, dat.so$Ms[ho], dat.so$Ms[he], dat.so$Ms[or])
2: number of rows of result
is not a multiple of vector length (arg 1) in: cbind(1, dat.so$Bws.d[ho], dat.so$Bws.d[he], dat.so$Bws.d[or])
> ##
> abline(h=0, col="darkgrey")
>
> points( dat.so$Ms[ho], dat.so$Bws.d[ho], pch=1 )
> points( dat.so$Ms[or], dat.so$Bws.d[or], pch=4, col=3, lwd=2 )
>
>
>
> ## Nice plot
>
> matplot(cbind(dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or]) / 1000,
+ cbind(dat.so$Bws[he], dat.so$Bws[ho], dat.so$Bws[or]) / 1000,
+ ylim=range(dat.so$Bws.p) / 1000,
+ pch=c(3,1,4),
+ col=c(2,1,3),
+ lwd=c(1,1,2),
+ xlab = expression(
+ "Subunit Molecular Weight " /~kDa ),
+ ylab = expression (
+ "Total Buried Area " / ring(A) * phantom()^2 * 10^-3 ),
+ cex.axis=1.3,
+ cex.lab=1.5
+ )
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or])
2: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Bws[he], dat.so$Bws[ho], dat.so$Bws[or])
> ##
> curve ( ((or.k3 + (x*1000) * or.k4) -
+ (to.k1 * (x*1000) ^ to.k2)) / 1000, add=T, col = 4, lwd = 2, lty = 1 )
> ##
> abline( h = 0, col = "darkgrey" )
> abline( v = 0, col = "darkgrey" )
>
> legend(x = "topleft", inset = 0.05, bg = "white", cex = 1.3,
+ c("Homo-Oligomers",
+ "Hetero-Oligomers",
+ "Original (21) Oligomers"),
+ pch=c(1,3,4),
+ col=c(1,2,3),
+ pt.lwd=c(1,1,2)
+ )
>
>
>
>
>
> ## Make a nice plot with the good bits from abvoe.
>
> postscript(file = "Image/plot.3p.ps", horizontal=FALSE)
>
> ## Reduce upper margin
> par("mar")
[1] 5.1 4.1 4.1 2.1
> par( mar = c( 5, 4, 1, 2 ) + 0.1 )
>
> ## Reset this
> par("mex")
[1] 1
> par( mex = 1.5 )
>
> ## GO
> split.screen( c( 2, 1 ), erase = T )
[1] 1 2
>
> ## PLOT
> screen( 1 )
> ##
> matplot(cbind(dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or]) / 1000,
+ cbind(dat.so$Bws[he], dat.so$Bws[ho], dat.so$Bws[or]) / 1000,
+ ylim=range(dat.so$Bws.p) / 1000,
+ pch=c(3,1,4),
+ col=c(2,1,3),
+ lwd=c(1,1,2),
+ xlab = expression(
+ "Subunit Molecular Weight " /~kDa ),
+ ylab = expression (
+ "Total Buried Area " / ring(A) * phantom()^2 * 10^-3 ),
+ cex.axis=1.3,
+ cex.lab=1.5
+ )
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or])
2: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Bws[he], dat.so$Bws[ho], dat.so$Bws[or])
> ##
> curve ( ((or.k3 + (x*1000) * or.k4) -
+ (to.k1 * (x*1000) ^ to.k2)) / 1000, add=T, col = 4, lwd = 2, lty = 1 )
> ##
> abline( h = 0, col = "darkgrey" )
> abline( v = 0, col = "darkgrey" )
>
> legend(x = "topleft", inset = 0.05, bg = "white", cex = 1.3,
+ c("Homo-Oligomers",
+ "Hetero-Oligomers",
+ "Original (21) Oligomers"),
+ pch=c(1,3,4),
+ col=c(1,2,3),
+ pt.lwd=c(1,1,2)
+ )
>
>
>
> ## PLOT
> screen( 2 )
> ##
>
> matplot(cbind(dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or]) / 1000,
+ cbind(dat.so$Bws.d[he], dat.so$Bws.d[ho], dat.so$Bws.d[or]),
+ pch=c(3,1,4),
+ col=c(2,1,3),
+ lwd=c(1,1,2),
+ xlab = expression ( "Subunit Molecular Weight " /~kDa ),
+ ylab = expression ( "Deviation in the Total Buried Area" ),
+ ##
+ cex.axis=1.3,
+ cex.lab=1.5
+ )
Warning messages:
1: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Ms[he], dat.so$Ms[ho], dat.so$Ms[or])
2: number of rows of result
is not a multiple of vector length (arg 2) in: cbind(1, dat.so$Bws.d[he], dat.so$Bws.d[ho], dat.so$Bws.d[or])
> ##
> abline( h = 0, col = 4, lwd = 2 )
> ##abline( h = 5, col = "darkgrey")
> ##abline( h = -5, col = "darkgrey")
> abline( v = 0, col = "darkgrey")
>
> legend(x = "bottomright", inset = 0.05, bg = "white", cex = 1.3,
+ c("Homo-Oligomers",
+ "Hetero-Oligomers",
+ "Original (21) Oligomers"),
+ pch=c(1,3,4),
+ col=c(1,2,3),
+ pt.lwd=c(1,1,2)
+ )
>
>
>
> ## DONE
> close.screen( all = TRUE )
>
> ## OK
> dev.off()
postscript
2
>
> ## DONE
> ## q()
>
>
>
>
>
>
>
> ##
> ## HISTO-GRAMMY
> ##
>
> ## Plot residuals (OLIGOMER)
> hist(dat.to$Ao.d,prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.to$Ao.d), sd=sd(dat.to$Ao.d)),
+ add=T,col=2,lwd=2)
>
>
> ## Plot residuals (SUBUNITS-HOMO)
> hist(dat.so$As.d[ho],prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.so$As.d[ho]), sd=sd(dat.so$As.d[ho])),
+ add=T,col=2,lwd=2)
>
> ## Plot residuals (SUBUNITS-HETERO)
> hist(
+ dat.so$As.d[he]
+ ,prob=T,breaks=20)
> ##
>
> curve(
+ dnorm(x,
+ mean=
+ mean(dat.so$As.d[he])
+ ,
+ sd= sd(dat.so$As.d[he])
+ ),
+ add=T,col=2,lwd=2)
>
>
>
> ##
> ## Test
> ##
>
> z <- rnorm(10000)
> ##
> my.seq <- seq(min(z),max(z),len=100)
> ##
> hist(z,prob=T,breaks=my.seq)
> ##
> curve(dnorm(x,mean=mean(z),
+ sd=sd(z)),
+ add=T,col=2,lwd=2)
>
>
>
>
>
> ##
> ## Plot the above three histograms on one plot
> ##
>
> hist(dat.to$Ao.d,prob=T,breaks=20, col=2,
+ xlim=c(-60,60),
+ ylim=c(0,0.06),
+ )
>
> hist(dat.so$As.d[ho],prob=T,breaks=20,col=3,add=T)
>
> hist(dat.so$As.d[he],prob=T,breaks=20,col=4,add=T)
>
>
>
>
>
> ##
> ## Plot all three distributions on one plot
> ##
>
> plot(1,1,xlim=c(-60,60),ylim=c(0,0.04),type="n")
> ##
> curve(dnorm(x, mean=mean(dat.to$Ao.d), sd=sd(dat.to$Ao.d)),
+ add=T,col=2,lwd=2)
> ##
> curve(dnorm(x, mean=mean(dat.so$As.d[ho]), sd=sd(dat.so$As.d[ho])),
+ add=T,col=3,lwd=2)
> ##
> curve(dnorm(x, mean=mean(dat.so$As.d[he]), sd=sd(dat.so$As.d[he])),
+ add=T,col=4,lwd=2)
> ##
> abline(v=mean(dat.to$Ao.d), lty=2, col=2)
> abline(v=mean(dat.so$As.d[ho]), lty=2, col=3)
> abline(v=mean(dat.so$As.d[he]), lty=2, col=4)
>
>
>
>
>
> ##
> ## Plot distributions of values for (visual) comparison (Tukey).
> ##
>
> boxplot(list("To (Olig)" = dat.to$Ao.d,
+ "Ho (Subs)" = dat.so$As.d[ho],
+ "He (Subs)" = dat.so$As.d[he]),
+ notch=T, varwidth=F,col=c(2,3,4),
+ ylab="Deviation in Ao")
> ##
> abline(h=mean(dat.to$Ao.d), lty=2, col=2)
> abline(h=mean(dat.so$As.d[ho]), lty=2, col=3)
> abline(h=mean(dat.so$As.d[he]), lty=2, col=4)
>
>
>
>
>
> ##
> ## Statistics of comparison
> ##
>
> (t.to.vs.ho <- t.test(dat.to$Ao.d, dat.so$As.d[ho] ))$p.value
[1] 6.646485e-17
> (t.to.vs.he <- t.test(dat.to$Ao.d, dat.so$As.d[he]))$p.value
[1] 1.413913e-39
> (t.ho.vs.he <- t.test(dat.so$As.d[ho], dat.so$As.d[he]))$p.value
[1] 3.537439e-11
>
>
>
> ##
> ## Test above statistics
> ##
>
> ## Tweek...
> x <- rnorm( 1000, 0.0, 1 )
> y <- rnorm( 1000, 0.3, 4 )
> ##
> plot( density(x,bw=0.5))
> lines(density(y,bw=0.5),col=2)
> ##
> t.test(x,y)$p.value
[1] 0.004663707
> ##
> abline(v=mean(x), lty=2)
> abline(v=mean(y), lty=2, col=2)
>
>
>
>
>
>
>
> ###
> ### Prepare an overall plot
> ###
>
> ## Use some R defaults!
> postscript(file = "Image/plot.3y.ps",
+ horizontal=FALSE)
>
> ## Store default par (for later)
> old.par <- par(no.readonly = TRUE)
>
> ## Give left margin a bit more space and remove top margin
> par( mar = c( 5, 4.5, 0, 2 ) + 0.1 )
>
> ## Divide up page...
> par(mfrow=c(4,1))
>
> ## Make labels big..
> par(cex.lab=1.5)
>
> ## Standardize plots
> my.xlim = c(-60,60)
>
>
>
>
>
> ## Plot histograms
>
> ## ONE
> hist(dat.to$Ao.d,prob=T,breaks=20,xlim=my.xlim,
+ main="",xlab="Deviation for oligomer ASA")
> ##
> curve(dnorm(x,
+ mean=mean(dat.to$Ao.d),
+ sd= sd(dat.to$Ao.d)),
+ add=T,col=2,lwd=2)
> ##
> abline(v=mean(dat.to$Ao.d), lty=2, col=2)
>
>
>
> ## TWO
> hist(dat.so$As.d[ho],prob=T,breaks=20,xlim=my.xlim,
+ main="",xlab="Deviation for homo-subunit ASA")
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$As.d[ho]),
+ sd= sd(dat.so$As.d[ho])),
+ add=T,col=3,lwd=2)
> ##
> abline(v=mean(dat.so$As.d[ho]), lty=2, col=3)
>
>
>
> ## THREE
> hist(dat.so$As.d[he],prob=T,breaks=20,xlim=my.xlim,
+ main="",xlab="Deviation for hetero-subunit ASA")
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$As.d[he]),
+ sd= sd(dat.so$As.d[he])),
+ add=T,col=4,lwd=2)
> ##
> abline(v=mean(dat.so$As.d[he]), lty=2, col=4)
>
>
>
> ## FOUR (plot the curves togehter)
> plot(1,1,ylim=c(0,0.04),type="n",xlim=my.xlim,
+ main="", xlab="Deviation",ylab="Density")
> ##
> curve(dnorm(x,
+ mean=mean(dat.to$Ao.d),
+ sd= sd(dat.to$Ao.d)),
+ add=T,col=2,lwd=2)
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$As.d[ho]),
+ sd= sd(dat.so$As.d[ho])),
+ add=T,col=3,lwd=2)
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$As.d[he]),
+ sd = sd(dat.so$As.d[he])),
+ add=T,col=4,lwd=2)
> ##
> abline(v=mean(dat.to$Ao.d ), lty=2, col=2)
> abline(v=mean(dat.so$As.d[ho]), lty=2, col=3)
> abline(v=mean(dat.so$As.d[he]), lty=2, col=4)
>
> ## OK
> par(old.par)
> dev.off()
postscript
2
> ## q()
>
>
>
>
>
>
>
> ###
> ### Investigate size of extreem hetero-subunits (larege deviation)
> ###
>
> ## Plot residuals (SUBUNITS-HETERO)
> hist(dat.so$As.d[he],prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.so$As.d[he]), sd=sd(dat.so$As.d[he])),
+ add=T, col=2, lwd=2)
> ##
> abline(v=30, col="grey", lwd=3)
>
> ## UPPER deviations
> tail( dat.so[he,][ order( dat.so$As.d[he] ), ], 6 )
Type Id So SGo chG Ms As Bbs
591 He 1gp2 3 3 G 5899 5675.3 2156.0
421 He 1n0w 2 2 B 3678 3789.8 1055.8
441 He 1dev 2 2 B 4110 4238.5 1296.4
226 He 1l0l 22 11 G 9289 8711.3 4006.2
891 He 1j1d 3 3 C 13406 12074.9 4399.4
1571 He 1ub4 6 2 C 8411 8075.6 3886.2
Name Ats Bws Bts As.p
591 G protein heterotrimer 9211.294 3535.994 5691.994 2981.689
421 . 5740.852 1951.052 3006.852 1968.092
441 . 6415.878 2177.378 3473.778 2169.989
226 . 14508.366 5797.066 9803.266 4444.933
891 . 20941.418 8866.518 13265.918 6137.266
1571 MazE / MazF antidote-toxin complex 13136.440 5060.840 8947.040 4073.293
As.d Bws.p Bws.d
591 47.46201 6229.605 -76.17692
421 48.06870 3772.760 -93.37051
441 48.80291 4245.889 -95.00013
226 48.97510 10063.433 -73.59529
891 49.17336 14804.152 -66.96692
1571 49.56049 9063.146 -79.08384
>
> plot(dat.so$As.d[he],
+ dat.so$Ms [he])
>
> abline(v=40, col="grey", lwd=3)
> abline(v=30,col="grey")
>
> abline(h=30000,col="grey")
> abline(h=20000,col="grey")
>
> ## Check two particular outliers...
> dat.so[he &
+ dat.so$As.d[he] > 30 &
+ dat.so$Ms [he] > 30000,]
Type Id So SGo chG Ms As Bbs
95 Ho 1hwx 6 1 A 55480 22538.50 5485.50
99 Ho 1dku 6 1 A 32299 13743.05 3899.85
174 He 1ohh 8 3 H 4007 3626.30 1682.00
178 He 1dd3 8 2 CD 3584 3070.70 1911.60
Name Ats
95 Glutamate dehydrogenase from Bos taurus 86684.487
99 Phosphoribosylpyrophosphate synthetase from B. subtilis 50462.828
174 . 6254.934
178 Ribosomal protein L12 from T. maritima 5593.972
Bws Bts As.p As.d Bws.p Bws.d
95 64145.987 69631.487 21399.383 5.054094 65285.10 -1.775820
99 36719.778 40619.628 13298.305 3.236146 37164.52 -1.211187
174 2628.634 4310.634 2122.094 41.480458 4132.84 -57.223860
178 2523.272 4434.872 1923.792 37.350047 3670.18 -45.453202
>
>
> ## HETERO-SUBUNIT OUTLIERS MW
>
> ## The mean molecular weight of subunits within a range of deviations
> mean( dat.so$Ms[ he & dat.so$As.d[he] > 30 ] )
[1] 32681.81
>
> ## The mean molecular weight of similar subunit blocks
> my.blocks <-
+ cut( dat.so$As.d[he], breaks=8, label=F )
>
> ## Check the size of each block
> table(my.blocks)
my.blocks
1 2 3 4 5 6 7 8
12 25 45 52 51 18 16 28
>
> ## Check average Mw of each block
> tapply(dat.so$Ms[he], my.blocks, mean)
1 2 3 4 5 6 7 8
55935.750 39676.200 35185.658 29557.035 24009.245 22543.111 15531.688 6540.393
>
>
>
>
>
>
>
> ##
> ## Is the difference in mean between the As.d of homo- and
> ## hetero-oligomers significantly different once we have removed the
> ## extreem cases? -- Only just.
> ##
>
> dat.so.trim <- dat.so$As.d < 30
>
>
>
> ###
> ### Statistics of comparison
> ###
>
> t.test( dat.to$Ao.d, dat.so$As.d[ he & dat.so.trim[he] ] )$p.value
[1] 5.596112e-39
> t.test( dat.so$As.d[ho], dat.so$As.d[ he & dat.so.trim[he] ] )$p.value
[1] 1.351140e-06
>
>
>
> ##
> ## And after removing all small hetero-subunits??
> ## -- Trimming at 10 KDa = Sig diff
> ## -- Trimming at 20 KDa = Sig diff (just)
> ##
>
> my.trim <- 0
> ##
> dat.so.trim <- dat.so$Ms > my.trim
> ##
> (t.test(dat.so$As.d[ ho & dat.so.trim[ho] ],
+ dat.so$As.d[ he & dat.so.trim[he] ]))$p.value
[1] 1
>
>
> my.trim <- 10000
> ##
> dat.so.trim <- dat.so$Ms > my.trim
> ##
> (t.test(dat.so$As.d[ ho & dat.so.trim[ho] ],
+ dat.so$As.d[ he & dat.so.trim[he] ]))$p.value
[1] 0.4662064
>
>
> my.trim <- 20000
> ##
> dat.so.trim <- dat.so$Ms > my.trim
> ##
> (t.test(dat.so$As.d[ ho & dat.so.trim[ho] ],
+ dat.so$As.d[ he & dat.so.trim[he] ]))$p.value
[1] 0.752284
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> ##
> ## Investigate size of extreem homo-subunits (larege deviation)
> ##
>
> ## Plot residuals (SUBUNITS-HETERO)
> hist(dat.so$As.d[ho],prob=T,breaks=20)
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$As.d[ho]),
+ sd= sd(dat.so$As.d[ho])
+ ),
+ add=T,col=2,lwd=2)
>
> plot(dat.so$As.d[ho],
+ dat.so$Ms [ho])
>
>
>
>
>
>
>
>
>
>
>
>
>
> ###
> ### Overall residuals (BSA)
> ###
>
> ## Predict values (OLIGOMER)
> dat.to$Btop <-
+ (or.k3 + dat.to$Mo * or.k4) -
+ (to.k1 * dat.to$Mo ^ to.k2)
>
> ## Predict values (SUBUNITS)
> dat.so$Bwsp <-
+ (or.k3 + dat.so$Ms * or.k4) -
+ (to.k1 * dat.so$Ms ^ to.k2)
>
> ## Turn predictions into percentage deviation (OLIGOMER)
> dat.to$Bto.d <-
+ (dat.to$Bto - dat.to$Btop) / dat.to$Bto * 100
>
> ## Turn predictions into percentages (SUBUNITS)
> dat.so$Bws.d <-
+ (dat.so$Bws - dat.so$Bwsp) / dat.so$Bws * 100
>
>
>
>
>
> ###
> ### PLOT
> ###
>
> ## Plot we have already seen...
> plot(dat.to$Ao.d, dat.to$Bto.d)
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
>
> ## New plots (SUBUNITS-HOMO)
> plot(dat.so$As.d[ho], dat.so$Bws.d[ho])
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
> ## New plots (SUBUNITS-HETERO)
> plot(dat.so$As.d[he], dat.so$Bws.d[he])
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
>
> ## Who is in extremis
> head( dat.so[he,][ order( dat.so$As.d[he] ), ], 1 )
Type Id So SGo chG Ms As Bbs
371 He 1cse 2 2 E 27240 9811.3 743.2
Name Ats Bws
371 Subtilisin Carlsberg / elastase inhibitor (eglin) complex 42557.85 32746.55
Bts As.p As.d Bws.p Bws.d Bwsp
371 33489.75 11448.25 -16.6843 31109.6 4.998837 31109.6
> tail( dat.so[he,][ order( dat.so$As.d[he] ), ], 1 )
Type Id So SGo chG Ms As Bbs Name
1571 He 1ub4 6 2 C 8411 8075.6 3886.2 MazE / MazF antidote-toxin complex
Ats Bws Bts As.p As.d Bws.p Bws.d Bwsp
1571 13136.44 5060.84 8947.04 4073.293 49.56049 9063.146 -79.08384 9063.146
>
> dat.to[dat.to$Id=="1ao7",]
Type Id So SGo Mo Ao Bbo
491 He 1ao7 5 5 80201 33161.4 8475.3
Name Ato
491 T-cell receptor / viral peptide / HLA-A complex from human 125312.5
Bto Bwo Ao.p Ao.d Bto.p Bto.d Btop
491 92151.09 83675.79 29589.53 10.77117 95722.96 -3.876101 95722.96
>
>
> ## OK, remove in extremis
> dat.he.trim <-
+ dat.so$As.d [he] > -20 &
+ dat.so$Bws.d[he] > -100
>
> plot(dat.so$As.d [ he & dat.he.trim ],
+ dat.so$Bws.d[ he & dat.he.trim ] )
>
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
>
>
> ## Plot both subunits together
>
> xmin <- min( dat.so$As.d [ho], dat.so$As.d [ he & dat.he.trim ] )
> xmax <- max( dat.so$As.d [ho], dat.so$As.d [ he & dat.he.trim ] )
> ymin <- min( dat.so$Bws.d[ho], dat.so$Bws.d[ he & dat.he.trim ] )
> ymax <- max( dat.so$Bws.d[ho], dat.so$Bws.d[ he & dat.he.trim ] )
>
> plot(dat.so$As.d [ho],
+ dat.so$Bws.d[ho],col=3,
+ xlab="Deviation in predicted Ao",
+ ylab="Deviation in predicted Bt",
+ xlim=c(xmin,xmax),
+ ylim=c(ymin,ymax),
+ cex.lab=1.5
+ )
>
> points(dat.so$As.d [ he & dat.he.trim ],
+ dat.so$Bws.d[ he & dat.he.trim ],col=4)
> points(dat.so$As.d [ho],
+ dat.so$Bws.d[ho],col=3)
>
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
>
>
> ## Plot all together
>
> xran <- range(dat.so$As.d[ho], dat.so$As.d [ he & dat.he.trim ], dat.to$Ao.d )
> yran <- range(dat.so$Bws.d[ho], dat.so$Bws.d[ he & dat.he.trim ], dat.to$Bto.d)
>
> plot(dat.so$As.d[ho], dat.so$Bws.d[ho],col=3,
+ xlab="Deviation in predicted Ao",
+ ylab="Deviation in predicted Bt",
+ xlim=xran,
+ ylim=yran,
+ cex.lab=1.5
+ )
>
> points(dat.so$As.d [ he & dat.he.trim ],
+ dat.so$Bws.d[ he & dat.he.trim ], col=4)
> points(dat.so$As.d [ho],
+ dat.so$Bws.d[ho], col=3)
> points(dat.to$Ao.d,
+ dat.to$Bto.d, col=2)
>
>
> abline(h=0,col="grey")
> abline(v=0,col="grey")
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> ###
> ### Duplicate A_S plots above in terms of A_B (deviations).
> ###
>
> ## Plot residuals (OLIGOMER)
> hist(dat.to$Bto.d,prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.to$Bto.d), sd=sd(dat.to$Bto.d)),
+ add=T,col=2,lwd=2)
>
> ## Plot residuals (SUBUNITS-HOMO)
> hist(dat.so$Bws.d[ho],prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.so$Bws.d[ho]), sd=sd(dat.so$Bws.d[ho])),
+ add=T,col=2,lwd=2)
>
> ## Plot residuals (SUBUNITS-HETERO)
> hist(dat.so$Bws.d[he],prob=T,breaks=20)
> ##
> curve(dnorm(x, mean=mean(dat.so$Bws.d[he]), sd=sd(dat.so$Bws.d[he])),
+ add=T,col=2,lwd=2)
>
>
>
>
>
> ## REMOVE 2 (SILLY) EXTREMES (FROM CONSIDERATION)
> my.i <-
+ dat.so$Bws.d > -200 &
+ dat.so$Bws.d < 12
>
>
>
>
>
> ## Plot residuals (SUBUNITS-HETERO) II
> hist( dat.so$Bws.d[ he & my.i ], prob=T, breaks=20 )
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$Bws.d[ he & my.i ]),
+ sd= sd(dat.so$Bws.d[ he & my.i ])),
+ add=T,col=2,lwd=2)
Warning messages:
1: longer object length
is not a multiple of shorter object length in: he & my.i
2: longer object length
is not a multiple of shorter object length in: he & my.i
>
>
>
>
>
>
> ###
> ### Plot distributions of values for (visual) comparison (Tukey).
> ###
>
> boxplot(list("To (Olig)" = dat.to$Bto.d,
+ "Ho (Subs)" = dat.so$Bws.d[ho],
+ "He (Subs)" = dat.so$Bws.d[he & my.i]),
+ notch=T, varwidth=F,col=c(2,3,4),
+ ylab="Deviation in Ao")
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
> ##
> abline(h=mean(dat.to$Bto.d), lty=2, col=2)
> abline(h=mean(dat.so$Bws.d[ho]), lty=2, col=3)
> abline(h=mean(dat.so$Bws.d[he & my.i]), lty=2, col=4)
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
>
>
>
> ## TRIMULA
> my.j <-
+ dat.so$Bws.d[he] > -20 &
+ dat.so$Bws.d[he] < 12
>
>
>
>
> ## REPLOT
> boxplot(list("To (Olig)" = dat.to$Bto.d,
+ "Ho (Subs)" = dat.so$Bws.d[ho],
+ "He (Subs)" = dat.so$Bws.d[he][my.j]),
+ notch=T, varwidth=F,col=c(2,3,4),
+ ylab="Deviation in Ao")
> ##
> abline(h=mean(dat.to$Bto.d), lty=2, col=2)
> abline(h=mean(dat.so$Bws.d[ho]), lty=2, col=3)
> abline(h=mean(dat.so$Bws.d[he][my.j]), lty=2, col=4)
>
>
>
>
> ###
> ### Statistics of comparison
> ###
>
> (t.to.vs.ho <- t.test(dat.to$Bto.d, dat.so$Bws.d[ho] ))$p.value
[1] 1.747294e-14
> (t.to.vs.he <- t.test(dat.to$Bto.d, dat.so$Bws.d[he & my.i]))$p.value
[1] 7.754711e-30
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
> (t.ho.vs.he <- t.test(dat.so$Bws.d[ho], dat.so$Bws.d[he & my.i]))$p.value
[1] 4.935567e-10
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
> (t.to.vs.he <- t.test(dat.to$Bto.d, dat.so$Bws.d[he][my.j]))$p.value
[1] 1.040179e-26
> (t.ho.vs.he <- t.test(dat.so$Bws.d[ho], dat.so$Bws.d[he][my.j]))$p.value
[1] 0.0933224
>
> ## Test
> (t.ho.vs.he <- t.test(dat.so$Bws.d[he],
+ dat.so$Bws.d[he & my.i]))$p.value
[1] 0.005095323
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
> ## Test
> (t.ho.vs.he <- t.test(dat.so$Bws.d[he & my.i],
+ dat.so$Bws.d[he][my.j]))$p.value
[1] 8.059372e-08
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
>
>
>
>
>
>
>
>
> ## TEST
> range(dat.to$Bto.d)
[1] -15.98346 10.26316
> range(dat.so$Bws.d[ho])
[1] -34.03630 3.46382
> range(dat.so$Bws.d[he & my.i])
[1] -95.000127 4.998837
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
> ## TEST II (NUMBER OF SUBUNITS BELOW THE RANGE OF THE OLIGOMERS)
> nrow(dat.so[ho & dat.so$Bws.d < -15.82116, ])
[1] 72
Warning message:
longer object length
is not a multiple of shorter object length in: ho & dat.so$Bws.d < -15.82116
> nrow(dat.so[he & dat.so$Bws.d < -15.82116, ])
[1] 72
Warning message:
longer object length
is not a multiple of shorter object length in: he & dat.so$Bws.d < -15.82116
>
> ## And with tiny trimmy?
> nrow(dat.so[he & dat.so$Bws.d < -15.82116 & my.i, ])
[1] 72
Warning message:
longer object length
is not a multiple of shorter object length in: he & dat.so$Bws.d < -15.82116
>
> ## And with big trimmy?
> nrow(dat.so[he & dat.so$Bws.d < -15.82116 & my.j, ])
[1] 57
Warning messages:
1: longer object length
is not a multiple of shorter object length in: he & dat.so$Bws.d < -15.82116
2: longer object length
is not a multiple of shorter object length in: he & dat.so$Bws.d < -15.82116 & my.j
>
>
>
>
>
>
>
>
>
>
>
>
>
> ###
> ### Prepare an overall plot
> ###
>
> ##
> ## BEGIN: Plot!
> ##
>
> ## Use some R defaults!
> postscript( file = "Image/plot.3z.ps", horizontal=FALSE )
>
> ## Store default par (for later)
> old.par <- par(no.readonly = TRUE)
>
> ## Give left margin a bit more space and remove top margin
> par( mar = c( 5, 4.5, 0.5, 2 ) + 0.1 )
>
> ## Divide up page and make labels / axis big..
> par(mfrow=c(3,1))
>
> par(cex.axis=1.5)
> par(cex.lab =1.7)
>
>
> ## Plot histograms
>
> source("legendLabel.R")
[1] "loading the legendLabel() function"
>
> ## ONE
> hist( dat.to$Bto.d, breaks=20,
+ main="",
+ xlab=expression( paste("Deviation in ", A[B] )),
+ ylab="Count"
+ )
>
> legendLabel(x="topleft", inset=0.2, cex=2.2, "Oligomers")
>
>
> ## TWO
> hist( dat.so$Bws.d[ho], breaks=20,
+ main="",
+ xlab=expression( paste("Deviation in ", A[B] )),
+ ylab="Count"
+ )
>
> legendLabel(x="topleft", inset=0.2, cex=2.2, "Homo-Subunits")
>
>
> ## THREE
> hist(dat.so$Bws.d[he & my.i], breaks=20,
+ main="",
+ xlab=expression( paste("Deviation in ", A[B] )),
+ ylab="Count"
+ )
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
> legendLabel(x="topleft", inset=0.2, cex=2.2, "Hetero-Subunits")
>
>
> ## OK
> par(old.par)
> dev.off()
postscript
2
> ## q()
>
>
>
>
>
>
>
>
>
> ###
> ### Investigate extreem hetero-subunits (larege deviation)
> ###
>
> ## Plot residuals (SUBUNITS-HETERO)
> hist(dat.so$Bws.d[he & my.i],prob=T,breaks=20)
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
> ##
> curve(dnorm(x,
+ mean=mean(dat.so$Bws.d[he & my.i]),
+ sd= sd(dat.so$Bws.d[he & my.i])),
+ add=T, col=2, lwd=2)
Warning messages:
1: longer object length
is not a multiple of shorter object length in: he & my.i
2: longer object length
is not a multiple of shorter object length in: he & my.i
> ##
> abline(v=-25, col="grey", lwd=3)
>
> ## LOWER deviations
> tail( dat.so[he,][ order( dat.so$Bws.d[he]), ], 6 )
Type Id So SGo chG Ms As Bbs
451 He 1ct4 2 2 E 18616 7446.0 616.5
1261 He 1n62 6 3 BE 87626 28194.6 5107.6
1811 He 1atn 2 2 D 28799 10764.2 885.3
601 He 1ehk 3 3 A 60305 20088.3 3990.0
971 He 1jnr 4 2 AC 73023 23442.9 6741.6
371 He 1cse 2 2 E 27240 9811.3 743.2
Name Ats
451 Proteinase B / ovomucoid inhibitor (third domain) 29082.35
1261 . 136914.48
1811 Actin / DNAse I complex 44993.87
601 . 94223.83
971 . 114096.44
371 Subtilisin Carlsberg / elastase inhibitor (eglin) complex 42557.85
Bws Bts As.p As.d Bws.p Bws.d Bwsp
451 21636.35 22252.85 8191.461 -10.01157 20890.88 3.445412 20890.88
1261 108719.88 113827.48 31985.456 -13.44533 104929.02 3.486811 104929.02
1811 34229.67 35114.97 12022.466 -11.68936 32971.41 3.675950 32971.41
601 74135.53 78125.53 23027.626 -14.63203 71196.20 3.964801 71196.20
971 90653.54 97395.14 27247.719 -16.23015 86848.72 4.197099 86848.72
371 32746.55 33489.75 11448.247 -16.68430 31109.60 4.998837 31109.60
>
> plot( dat.so$Bws.d[he & my.i],dat.so$Ms[he & my.i])
Warning messages:
1: longer object length
is not a multiple of shorter object length in: he & my.i
2: longer object length
is not a multiple of shorter object length in: he & my.i
>
> abline(v=-25, col="grey", lwd=3)
> abline(v=-20, col="grey")
>
> abline(h=30000,col="grey")
> abline(h=20000,col="grey")
>
> ## Check two particular outliers...
> dat.so[dat.so$Bws.d < -20 &
+ dat.so$Ms > 30000,]
Type Id So SGo chG Ms As Bbs Name Ats Bws Bts
1171 He 1iw7 6 6 F 39554 24323.5 5776.3 . 61799.19 37475.69 43251.99
1212 He 1k28 6 2 A 61269 39351.8 20285.0 . 95730.14 56378.34 76663.34
As.p As.d Bws.p Bws.d Bwsp
1171 15892.12 34.66353 45907.07 -22.49828 45907.07
1212 23351.02 40.66087 72379.12 -28.38109 72379.12
>
>
> ## HETERO-SUBUNIT OUTLIERS MW
>
> ## The mean molecular weight of subunits within a range of deviations
> mean(dat.so$Ms[he][dat.so$Bws.d[he]< -20])
[1] 10270.13
>
> ## The mean molecular weight of similar subunit blocks
> my.blocks <-
+ cut( dat.so$Bws.d[he & my.i], breaks=8, label=F )
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
> ## Check the size of each block
> table(my.blocks)
my.blocks
1 2 3 4 5 6 7 8
3 5 12 12 8 20 102 275
>
> ## Check average Mw of each block
> tapply(dat.so$Ms[he & my.i], my.blocks, mean)
1 2 3 4 5 6 7 8
3457.667 7245.800 6130.125 6868.125 9593.000 16263.475 22816.969 38513.139
Warning message:
longer object length
is not a multiple of shorter object length in: he & my.i
>
>
>
>
>
>
>
> ##
> ## Is the difference in mean between the Bws.d of homo- and
> ## hetero-oligomers significantly different once we have removed the
> ## extreem cases? -- Only just.
> ##
>
> dat.he.trim <-
+ dat.he.sub[my.i,][dat.so$Bws.d[he] < -20,]
>
>
>
> ###
> ### Statistics of comparison
> ###
>
> t.test(dat.to$Bto.d, dat.so$Bws.d[ he & dat.he.trim ])$p.value
[1] 1.724498e-18
Warning messages:
1: & not meaningful for factors in: Ops.factor(left, right)
2: & not meaningful for factors in: Ops.factor(left, right)
3: & not meaningful for factors in: Ops.factor(left, right)
4: & not meaningful for factors in: Ops.factor(left, right)
5: & not meaningful for factors in: Ops.factor(left, right)
6: & not meaningful for factors in: Ops.factor(left, right)
> t.test(dat.so$Bws.d[ho], dat.so$Bws.d[ he & dat.he.trim ])$p.value
[1] 6.925334e-10
Warning messages:
1: & not meaningful for factors in: Ops.factor(left, right)
2: & not meaningful for factors in: Ops.factor(left, right)
3: & not meaningful for factors in: Ops.factor(left, right)
4: & not meaningful for factors in: Ops.factor(left, right)
5: & not meaningful for factors in: Ops.factor(left, right)
6: & not meaningful for factors in: Ops.factor(left, right)
>
>
>
> ##
> ## And after removing all small hetero-subunits??
> ## -- Sig diff after trimming at 10 KDa
> ## -- Barly sig diff after trimming at 20 KDa
> ##
>
> my.trim <- 0
> my.trim <- 10000
> my.trim <- 20000
>
> dat.he.sub.trim <- dat.he.sub[dat.so$Ms[he]>my.trim,]
> dat.ho.trim <- dat.ho [dat.ho $Ms>my.trim,]
> ##
> nrow(dat.he.sub.trim)
[1] 138
> nrow(dat.ho.trim)
[1] 120
> ##
> (t.test(dat.ho.trim$Bws.d, dat.so$Bws.d[ he & dat.he.trim ]))$p.value
Error in t.test.default(dat.ho.trim$Bws.d, dat.so$Bws.d[he & dat.he.trim]) :
not enough 'x' observations
In addition: Warning messages:
1: & not meaningful for factors in: Ops.factor(left, right)
2: & not meaningful for factors in: Ops.factor(left, right)
3: & not meaningful for factors in: Ops.factor(left, right)
4: & not meaningful for factors in: Ops.factor(left, right)
5: & not meaningful for factors in: Ops.factor(left, right)
6: & not meaningful for factors in: Ops.factor(left, right)
7: is.na() applied to non-(list or vector) in: is.na(x)
Execution halted