重画boxplot去掉自交 "2016-10-20 22:14:09 CST"

因为要去掉数据中的自交,所以之前画的boxplot需要重新做,代码如下:

setwd("D:/Sync/R/8801/fullsib/160620")
Myfile<-c("8801fs_05_14.csv")
df<-read.csv(Myfile)
str(df)
# 'data.frame': 2226 obs. of  20 variables:
#   $ TreeId : int  3141 3142 3143 3144 3145 3146 3147 3148 3149 3150 ...
# $ Mum    : Factor w/ 23 levels "CF29","R1","R1_OP",..: 18 18 18 18 18 18 18 18 18 18 ...
# $ Dad    : Factor w/ 13 levels "","CF29","R1",..: 8 8 8 8 8 8 8 8 8 8 ...
# $ Com    : Factor w/ 66 levels "CF29RC402","CF29Y11",..: 52 52 52 52 52 52 52 52 52 52 ...
# $ Tree   : int  NA 123 NA 107 NA 120 104 122 118 121 ...
# $ Fam    : int  84 84 84 84 84 84 84 84 84 84 ...
# $ Rep    : int  1 1 1 1 1 1 1 1 1 1 ...
# $ Row    : int  41 41 41 42 42 42 43 43 44 44 ...
# $ Col    : int  52 53 55 51 54 55 52 54 51 53 ...
# $ Recipro: Factor w/ 18 levels "",".","RC13RC82",..: 12 12 12 12 12 12 12 12 12 12 ...
# $ Plot   : Factor w/ 250 levels "","1CF29Y11",..: 42 42 42 42 42 42 42 42 42 42 ...
# $ h05    : num  15.5 15 15.5 16.5 17.5 15.8 17.2 17.2 15.5 16.1 ...
# $ d05    : num  14.5 16 13 13.3 13.4 14.7 18.3 17.3 14.5 14.8 ...
# $ v05    : num  0.326 0.384 0.262 0.292 0.314 ...
# $ h14    : num  NA 19.9 NA 20.7 NA 20.8 23.2 23.1 19.2 21.6 ...
# $ d14    : num  NA 20.2 NA 15.9 NA 19.1 23.1 24 17.1 17.4 ...
# $ v14    : num  NA 0.812 NA 0.523 NA ...
# $ delta.h: num  NA 4.9 NA 4.2 NA 5 6 5.9 3.7 5.5 ...
# $ delta.d: num  NA 4.2 NA 2.6 NA 4.4 4.8 6.7 2.6 2.6 ...
# $ delta.v: num  NA 0.428 NA 0.231 NA ...
names(df)
# [1] "TreeId"  "Mum"     "Dad"     "Com"     "Tree"    "Fam"     "Rep"     "Row"    
# [9] "Col"     "Recipro" "Plot"    "h05"     "d05"     "v05"     "h14"     "d14"    
# [17] "v14"     "delta.h" "delta.d" "delta.v"
names(df)[12:length(names(df))]<-c("H05","DBH05","V05","H14","DBH14","V14","HI", "DI","VI")
df$Mum<-as.character(df$Mum)
df$Dad<-as.character(df$Dad)
# df$Mum<-factor(df$Mum)
# df$Dad<-factor(df$Dad)
df.rm.self<-df[df$Mum!=df$Dad,]
levels(df.rm.self$Com)
str(df.rm.self)
df.rm.self[df.rm.self$Com=="Y9Y9",]
df.rm.self$Com<-as.character(df.rm.self$Com)
df.rm.self$Com<-factor(df.rm.self$Com)
# par(las=2)
# plot(df.rm.self$H05~df.rm.self$Com)

source("D:/Sync/R/rcodes/d.box.order.r")
names(df.rm.self)[4]<-c("com")
#pdf
pdf("grow.all.crosses.order.pdf", width = 7, height = 4,family="Times",pointsize = 8)
par(mar = c(6, 4, .5, .5),las=2)
X<-names(df.rm.self) 
for (j in 12:length(X)){
  data.i<-df.rm.self[, X[j]]
  d.box.order(data.i, df.rm.self, X, j)
}
dev.off()
# win.metafile
win.metafile("grow.all.crosses.H05.wmf", width = 7, height = 4,family="serif",pointsize = 8)
par(mar = c(6, 4, .5, .5),las=2)
X<-names(df.rm.self) 
  data.i<-df.rm.self[, X[12]]
  d.box.order(data.i, df.rm.self, X, 12)
dev.off()

com<-unique(df$Com)
com<-as.character(com)
comcom<-read.csv("comcom.growth.csv")
# comcom<-read.csv("comcom.growth.rmself.csv")
comcom<-as.matrix(comcom)
names(df)[4]<-c("com")
pdf("grow.performance.rmself.pdf",paper="a4", width = 0, height = 9,family="Times",pointsize = 10)
par(mfrow=c(3,3),mar = c(6, 4, .5, .5),las=2)
for (i in 1:nrow(comcom)){
  df.com<-df[df$com %in% com[comcom[i,]],]
  df.com$com<-factor(df.com$com)
  X<-names(df.com) 
  for (j in 12:length(X)){
    data.i<-df.com[, X[j]]
    model=lm(data.i ~ df.com$com )
    ANOVA=aov(model)
    TUKEY <- TukeyHSD(x=ANOVA, 'df.com$com', conf.level=0.99)
    # a<-TUKEY$`df.com$com`[,c(1,4)]
    # colnames(a)<-c(paste(X[j],"diff",sep = "_"),"p adj")
    # write.table(a,file="growth_p_value.csv",sep=",",quote=F,append=T)
    LABELS=generate_label_df(TUKEY , "df.com$com")
    a=boxplot(data.i~df.com$com,
              ylim=c(min(data.i,na.rm=T) ,
                     1.01*max(data.i,na.rm=T)) , 
              #xlab="Crosses",
              ylab=X[j] , main="")
    abline(h=mean(data.i,na.rm=T),lty=2)
    # I want to write the letter over each box. Over is how high I want to write it.
    over=0.02*max( a$stats[nrow(a$stats),] )
    #Add the labels
    text( c(1:nlevels(df.com$com)) , 
          a$stats[nrow(a$stats),]+over , LABELS[,1])
  }}
dev.off()

# win.metafile
win.metafile("grow.performance.rmself.wmf", width = 8.3, height = 11.7,family="serif",pointsize = 10)
par(mfrow=c(3,3),mar = c(6, 4, .5, .5),las=2)
  df.com<-df[df$com %in% com[comcom[1,]],]
  df.com$com<-factor(df.com$com)
  X<-names(df.com) 
  for (j in 12:length(X)){
    data.i<-df.com[, X[j]]
    model=lm(data.i ~ df.com$com )
    ANOVA=aov(model)
    TUKEY <- TukeyHSD(x=ANOVA, 'df.com$com', conf.level=0.99)
    # a<-TUKEY$`df.com$com`[,c(1,4)]
    # colnames(a)<-c(paste(X[j],"diff",sep = "_"),"p adj")
    # write.table(a,file="growth_p_value.csv",sep=",",quote=F,append=T)
    LABELS=generate_label_df(TUKEY , "df.com$com")
    a=boxplot(data.i~df.com$com,
              ylim=c(min(data.i,na.rm=T) ,
                     1.01*max(data.i,na.rm=T)) , 
              #xlab="Crosses",
              ylab=X[j] , main="")
    abline(h=mean(data.i,na.rm=T),lty=2)
    # I want to write the letter over each box. Over is how high I want to write it.
    over=0.02*max( a$stats[nrow(a$stats),] )
    #Add the labels
    text( c(1:nlevels(df.com$com)) , 
          a$stats[nrow(a$stats),]+over , LABELS[,1])
  }
dev.off()
save(df,df.rm.self,file = "8801.未替换家系名及去掉自交.RData")
Sys.time()
# [1] "2016-10-20 22:14:09 CST"

推荐阅读更多精彩内容