####################
#for this code to work you must first insure that limma, marray, and convert are downloaded from bioconducter
#use the following 4 lines if you do not have them
####################

#source("http://bioconductor.org/biocLite.R")
#biocLite("limma")
#biocLite("marray")
#biocLite("convert")
#biocLite("siggenes") #not used until later



###################
#in the working directory you must have
#     all 34 arrays in tab delim txt format titled "array (01).txt","array (02).txt", etc
#     "toReadIn.txt" as specified on my website
###################


setwd("C:/Documents and Settings/Austen/My Documents/~School/~0708Spring/math155/project/tab")


library(marray)
library(limma)
dat.targets<-readTargets("toReadIn.txt")


# to keep the data return 1, to throw it out return 0
dat.filter<-function(x, intensity.threshold=30, percent.threshold=0, no.spotflag.val=0, no.autoflag.val=0){

     # throw out spot if the percent.threshold percent of pixels in both green and red foregrounds
     # is less than the background + 1 standard deviation of each color
     brighter.than.bg<-as.numeric(x[,"% CH1 PIXELS > BG + 1SD"] > percent.threshold ||
                        x[,"% CH2 PIXELS > BG + 1SD"] > percent.threshold)

     # throw out the spot if the scanning software suggests we should
     autoflag<-as.numeric(x[,"Autoflag"] == no.autoflag.val)

     # throw out the spot if both colors' median intensities are less than intensity.threshold
     intensity.median<-as.numeric(x["Ch1 Intensity (Median)"] >= intensity.threshold ||
                                  x["Ch2 Intensity (Median)"] >= intensity.threshold)

     # throw out the spot if it was flagged in the creation of the data
     spotflag<-as.numeric(x[,"Spot Flag"] == no.spotflag.val)

     # throw out controls
     control.names<-c("CONTROL_1", "CONTROL_17", "CONTROL_2",
          "CONTROL_21", "CONTROL_23", "CONTROL_27", "CONTROL_28",
          "CONTROL_29","CONTROL_3", "CONTROL_30", "CONTROL_31",
          "CONTROL_34", "CONTROL_36", "CONTROL_37", "CONTROL_38",
          "CONTROL_4", "CONTROL_44", "CONTROL_45", "CONTROL_48",
          "CONTROL_49", "CONTROL_7", "CONTROL_8")
     controls<-1
     for(i in 1:length(control.names)){
          controls<-controls*as.numeric(x[,"Reporter Name"] != control.names[i])
     }

     return(brighter.than.bg * autoflag * intensity.median * spotflag * controls)

}


full.dat<-read.maimages(dat.targets$Name,source="smd",wt.fun=dat.filter)


############### foreground and background intensities boxplots

win.graph()
par(mfrow=c(2,2))
boxplot(data.frame(log2(full.dat$R)),main="log2 Red foreground",names=(1:34),xlab="Array #")
boxplot(data.frame(log2(full.dat$Rb)),main="log2 Red background",names=(1:34),xlab="Array #")
boxplot(data.frame(log2(full.dat$G)),main="log2 Green foreground",names=(1:34),xlab="Array #")
boxplot(data.frame(log2(full.dat$Gb)),main="log2 Green background",names=(1:34),xlab="Array #")
mtext("Raw Data", outer=T, line=-2)




################################## MA PLOTS

################ raw data

this.method<-"none"
full.dat.ma<-normalizeWithinArrays(full.dat, method=this.method)


win.graph()

temp<-rbind(50:56,cbind(37:42,t(array(1:36,dim=c(6,6)))),c(49,43:48))
layout(temp,widths=c(.5,rep(1,ncol(temp)-1)),
     heights=c(.5,rep(1,nrow(temp)-2),.5))
par(mar=c(0,0,1.2,0))

for(i in 1:34){
     plot(full.dat.ma$A[,i],full.dat.ma$M[,i],
          xlim=c(0,16),ylim=c(-9,9),
          main=NA,xlab=NA,ylab=NA,pch=".",axes=F,frame.plot=T)
     title(main=paste("array",i),line=0.2)
     
}
for(i in 35:36){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F)}
for(i in 37:42){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(2,line=-3.5)}
for(i in 43:48){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(1,line=-3.5)}
mtext("Raw Data", outer=T,line=-2)

################ background subtracted data (median)



this.method<-"median"
full.dat.bgd<-normalizeWithinArrays(full.dat, method=this.method)


win.graph()

temp<-rbind(50:56,cbind(37:42,t(array(1:36,dim=c(6,6)))),c(49,43:48))
layout(temp,widths=c(.5,rep(1,ncol(temp)-1)),
     heights=c(.5,rep(1,nrow(temp)-2),.5))
par(mar=c(0,0,1.2,0))

for(i in 1:34){
     plot(full.dat.bgd$A[,i],full.dat.bgd$M[,i],
          xlim=c(0,16),ylim=c(-9,9),
          main=NA,xlab=NA,ylab=NA,pch=".",axes=F,frame.plot=T)
     title(main=paste("array",i),line=0.2)
     
}
for(i in 35:36){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F)}
for(i in 37:42){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(2,line=-3.5)}
for(i in 43:48){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(1,line=-3.5)}
mtext("Background Subtracted Data", outer=T,line=-2)






################ loess normalized data



this.method<-"loess"
full.dat.norm<-normalizeWithinArrays(full.dat, method=this.method)


win.graph()

temp<-rbind(50:56,cbind(37:42,t(array(1:36,dim=c(6,6)))),c(49,43:48))
layout(temp,widths=c(.5,rep(1,ncol(temp)-1)),
     heights=c(.5,rep(1,nrow(temp)-2),.5))
par(mar=c(0,0,1.2,0))

for(i in 1:34){
     plot(full.dat.norm$A[,i],full.dat.norm$M[,i],
          xlim=c(0,16),ylim=c(-9,9),
          main=NA,xlab=NA,ylab=NA,pch=".",axes=F,frame.plot=T)
     title(main=paste("array",i),line=0.2)
     
}
for(i in 35:36){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F)}
for(i in 37:42){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(2,line=-3.5)}
for(i in 43:48){plot(NA,xlim=c(0,16),ylim=c(-9,9),axes=F);axis(1,line=-3.5)}
mtext("Loess Normalized Data", outer=T,line=-2)





win.graph()
boxplot(data.frame(full.dat.norm$M),main="Normalized Data",names=(1:34),xlab="Array #")




############################################# TEST FOR NORMALIZED DATA


##### ks.test

ks.test.pvals<-c()
for(i in 1:nrow(full.dat.norm$M)){
     ks.test.pval<-ks.test(full.dat.norm$M[i,],
          "pnorm",mean(full.dat.norm$M[i,],na.rm=T),sd(full.dat.norm$M[i,],na.rm=T))$p.value
     ks.test.pvals<-c(ks.test.pvals,ks.test.pval)
}
hist(ks.test.pvals)



###### looney test

####################### this function was written by John Kloke \/ \/ \/ \/

looney<-function(y,nbs=2000,alpha=0.05) {

############################################################
#
#  looney - get the correlation of the data and their 
#    expected value under normality 
#    estimate the pval and critical value via parametric bootstrap
#    for testing 
#      H0: Data are Normal
#      HA: Data are not Normal
#    ie no need for the Looney Table
#
#  input 
#    y - data you would like tested for normality 
#    nbs - number of bootstrap samples to take
#    alpha - approximates a cv for this size test
#
#  output
#    cor - the correlation between the data and there 
#      expected value under normality
#    pval - a bootstrap estimate of the pvalue for th
#    ecv - the estimated critical value 
#      (should be quite close to the value in the Looney Table)
#
############################################################

	n<-length(y)
	ey<-qnorm((rank(y) - 0.375)/(n+0.25))
	est<-cor(ey,y)

	cors<-rep(0,nbs)

	for( i in 1:nbs ) {

		z<-rnorm(n)
		ez<-qnorm((rank(z) - 0.375)/(n+0.25))

		cors[i]<-cor(ez,z)

	}

	pval=sum( cors < est )/nbs
	ecv=sort(cors)[nbs*alpha]

	list(cor=est, pval=pval, ecv=ecv)

}


####################################### /\ /\ /\ /\ end Kloke's function



####### NOTE: this is commented out only because it takes a long time to run
## if you want looney test p vals then uncomment the following 6 lines

#looney.test.pvals<-c()
#for(i in 1:nrow(full.dat.norm$M)){
#     looney.test.pval<-looney(rm.na(full.dat.norm$M[i,]),nbs=1000,alpha=0.1)$pval
#     looney.test.pvals<-c(looney.test.pvals,looney.test.pval)
#}
#hist(looney.test.pvals)



##############################


spot.miRNA<-c()
for(i in 1:32){spot.miRNA<-c(spot.miRNA,rep(((24*(i-1)+1):(24*i)),2))}
mean.2<-function(x){mean(x,na.rm=T)}
norm.M<-array(NA,dim=c(nrow(full.dat.norm$M)/2,ncol(full.dat.norm$M)))
for(i in 1:ncol(full.dat.norm$M)){
     norm.M[,i]<-tapply(full.dat.norm$M[,i],spot.miRNA,mean.2)
}


reporter.name.1536<-abbreviate(read.delim("array (01).txt",skip=15)$Reporter.Name,25)
reporter.name<-c()
for(i in 1:32){reporter.name<-c(reporter.name,reporter.name.1536[(1:24)+(24*(i-1))])}






design<-modelMatrix(dat.targets,ref="UHR")

fit0<-lmFit(norm.M, design)
cont.matrix1<-makeContrasts(ERMS-ARMS, levels=design)
cont.matrix2<-makeContrasts(PRMS-GIST, levels=design)

fit1<-contrasts.fit(fit0, cont.matrix1)
fit2<-contrasts.fit(fit0, cont.matrix2)

fit1B<-eBayes(fit1)
fit2B<-eBayes(fit2)

dat.topTable1<-topTable(fit1B,genelist=reporter.name,number=25)
dat.topTable2<-topTable(fit2B,genelist=reporter.name,number=25)

dat.topTable1
dat.topTable2



##### the following 4 functions write the topTable to an excel file and a txt file in the working directory

#write.xls(dat.topTable1,file="topTenGenes1.xls")
#write.xls(dat.topTable1,file="topTenGenes1.txt")
#write.xls(dat.topTable2,file="topTenGenes2.xls")
#write.xls(dat.topTable2,file="topTenGenes2.txt")


#these are the significant points as noted from the two topTables

volcanoplot(fit1B,main="ARMS vs ERMS")
text(-3.4,4.5,"a")
text(-3.4,4.7,"a")
text(5.5,3.5,"b")
text(4.9,1.7,"b")
text(4.4,3,"d")
text(4.2,2.3,"d")


volcanoplot(fit2B,main="PRMS vs GIST")
text(4.7,9.8,"a")
text(4.6,8.4,"a")
text(-5.6,4.8,"b")
text(-5.6,4.6,"b")
text(-4.3,3.6,"c")
text(-4.2,3.6,"c")
text(4.6,3.6,"d")
text(4.5,3.6,"d")






#########################
# code to put spots to miRNA correctly (as on hw 7):





target.normals<-which(dat.targets$DiseaseState=="normal")
target.tumors<-which(dat.targets$DiseaseState=="tumor")
p.vals<-c()
t.vals<-c()
for(i in 1:nrow(norm.M)){
ifelse(min(length(rm.na(norm.M[i,target.normals])),(length(rm.na(norm.M[i,target.tumors]))))>1,
     {p.vals<-c(p.vals,t.test(norm.M[i,target.normals],norm.M[i,target.tumors])$p.value);
       t.vals<-c(t.vals,t.test(norm.M[i,target.normals],norm.M[i,target.tumors])$statistic)},
     {p.vals<-c(p.vals,NA);t.vals<-c(t.vals,NA)})
}




#################### project 7
library(siggenes)





#cl.list<-c(1,1,2,3,1,1,2,4,4,3,5,5,6,5,7,7,4,3,3,3,1,8,1,1,1,4,3,3,4,9,9,9,9,9)

cl.list.all<-c(1,1,2,3,1,1,2,4,4,3,5,5,5,6,6,4,3,3,3,1,1,1,1,4,3,3,4,7,7,7,7,7)
norm.sam.all<-sam(norm.M[,-c(13,22)],cl.list.all,gene.names=reporter.name)
plot(norm.sam.all,2.01)
summary(norm.sam.all,2.01)



cl.list.tumor<-replace(rep(2,34),target.normals,1)
norm.sam.tumor<-sam(norm.M,cl.list.tumor,gene.names=reporter.name)
plot(norm.sam.tumor,1.3)
mtext("normal vs tumor tissue",line=.2)
summary(norm.sam.tumor,1.3)







##################### project 8




design.tumor<-cbind(as.numeric(dat.targets$DiseaseState=="normal"),
                   as.numeric(dat.targets$DiseaseState=="tumor"))
colnames(design.tumor)<-c("normal","tumor")


fit0.tumor<-lmFit(norm.M, design.tumor)
cont.matrix.tumor<-makeContrasts(normal-tumor, levels=design.tumor)
fit1.tumor<-contrasts.fit(fit0.tumor, cont.matrix.tumor)
fit1B.tumor<-eBayes(fit1.tumor)
dat.topTable.tumor<-topTable(fit1B.tumor,genelist=reporter.name,number=20)
volcanoplot(fit1B.tumor,main="normal - tumor")
abline(h=dat.topTable.tumor$B[20],col="grey")
abline(v=0,col="grey")

rna.sig<-as.numeric(rownames(dat.topTable.tumor))
rna.rand<-sample(1:768,20,rep=F)

rna.sig.k2.true<-as.numeric(dat.topTable.tumor$logFC>0)+1 #1 is logFC neg, 2 is logFC pos

rna.sig.k2.true

################ dendrograms 1: 1-abs(cor), complete

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<-"complete"
this.dist<-"1-abs(cor)"


####significant miRNA
rna.sig.dists<-as.dist(1-abs(cor(t(norm.M[rna.sig,]),use="complete.obs")))
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)


####random miRNA
rna.rand.dists<-as.dist(1-abs(cor(t(norm.M[rna.rand,]),use="complete.obs")))
rna.rand.clusts<-hclust(rna.rand.dists,method=this.method)
rna.rand.k2<-cutree(rna.rand.clusts,k=2)

plot(rna.rand.clusts,
     labels=rna.rand,
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for random miRNA",
     sub="")




################ dendrograms 2: 1-abs(cor), average

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<-"average"
this.dist<-"1-abs(cor)"


####significant miRNA
rna.sig.dists<-as.dist(1-abs(cor(t(norm.M[rna.sig,]),use="complete.obs")))
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)


####random miRNA
rna.rand.dists<-as.dist(1-abs(cor(t(norm.M[rna.rand,]),use="complete.obs")))
rna.rand.clusts<-hclust(rna.rand.dists,method=this.method)
rna.rand.k2<-cutree(rna.rand.clusts,k=2)

plot(rna.rand.clusts,
     labels=rna.rand,
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for random miRNA",
     sub="")



################ dendrograms 3: 1-(cor)^2, complete

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<-"complete"
this.dist<-"1-(cor)^2"


####significant miRNA
rna.sig.dists<-as.dist(1-(cor(t(norm.M[rna.sig,]),use="complete.obs"))^2)
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)


####random miRNA
rna.rand.dists<-as.dist(1-(cor(t(norm.M[rna.rand,]),use="complete.obs"))^2)
rna.rand.clusts<-hclust(rna.rand.dists,method=this.method)
rna.rand.k2<-cutree(rna.rand.clusts,k=2)

plot(rna.rand.clusts,
     labels=rna.rand,
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for random miRNA",
     sub="")





################ dendrograms 4: 1-(cor)^2, average

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<-"average"
this.dist<-"1-(cor)^2"


####significant miRNA
rna.sig.dists<-as.dist(1-(cor(t(norm.M[rna.sig,]),use="complete.obs"))^2)
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)


####random miRNA
rna.rand.dists<-as.dist(1-(cor(t(norm.M[rna.rand,]),use="complete.obs"))^2)
rna.rand.clusts<-hclust(rna.rand.dists,method=this.method)
rna.rand.k2<-cutree(rna.rand.clusts,k=2)

plot(rna.rand.clusts,
     labels=rna.rand,
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for random miRNA",
     sub="")




################ dendrograms 5: 1-(cor)^2, average

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<- "mcquitty"
this.dist<-"1-(cor)^2"


####significant miRNA
rna.sig.dists<-as.dist(1-(cor(t(norm.M[rna.sig,]),use="complete.obs"))^2)
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)



################ dendrograms 6: 1-abs(cor), mcquitty

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<- "mcquitty"
this.dist<-"1-abs(cor)"


####significant miRNA
rna.sig.dists<-as.dist(1-abs(cor(t(norm.M[rna.sig,]),use="complete.obs")))
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)



####random miRNA
rna.rand.dists<-as.dist(1-(cor(t(norm.M[rna.rand,]),use="complete.obs"))^2)
rna.rand.clusts<-hclust(rna.rand.dists,method=this.method)
rna.rand.k2<-cutree(rna.rand.clusts,k=2)

plot(rna.rand.clusts,
     labels=rna.rand,
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for random miRNA",
     sub="")


####################### end project 8

####### begin project 9 ###




dat.topTable.tumor100<-topTable(fit1B.tumor,genelist=reporter.name,number=100)
rna.sig<-as.numeric(rownames(dat.topTable.tumor100))
rna.sig.k2.true<-as.numeric(dat.topTable.tumor100$logFC>0)+1 #1 is logFC neg, 2 is logFC pos
rna.sig.k2.true





################ dendrogram top 100: 1-abs(cor), complete

#method can be: "ward", "single", "complete", "average", "mcquitty", "median" or "centroid".
this.method<-"complete"
this.dist<-"1-abs(cor)"


####significant miRNA
rna.sig.dists<-as.dist(1-abs(cor(t(norm.M[rna.sig,]),use="complete.obs"))) #'complete' for cor
rna.sig.clusts<-hclust(rna.sig.dists,method=this.method)
rna.sig.k2<-cutree(rna.sig.clusts,k=2)

plot(rna.sig.clusts,
     labels=paste(ifelse(rna.sig.k2.true==1,"(-)","(+)"),rna.sig),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")
mtext("(+): logFC(normal-tumor)>0, (-): logFC(normal-tumor)<0",side=1,line=1)



################

require(cluster)

max.disses<-c()
for(i in 2:20){
     max.disses<-c(max.disses,mean(pam(rna.sig.dists,i)$clusinfo[,"max_diss"]))
}

plot(2:20,max.disses,xlab="# clusters",ylab="max dissimilarity",main="MAX diss vs # Clusters")

avg.disses<-c()
for(i in 2:20){
     avg.disses<-c(avg.disses,mean(pam(rna.sig.dists,i)$clusinfo[,"av_diss"]))
}


plot(2:20,avg.disses,xlab="# clusters",ylab="average dissimilarity",main="AVERAGE diss vs # Clusters")


avg.sils<-c()
for(i in 2:20){
     avg.sils<-c(avg.sils,pam(rna.sig.dists,i)$silinfo$avg.width)
}

plot(2:20,avg.sils,xlab="# clusters",ylab="average silhouette width",main="Average Silhouette Width vs # Clusters")
points(8,avg.sils[7],lwd=3)



clus8<-pam(rna.sig.dists,8)$clustering
clus9<-pam(rna.sig.dists,9)$clustering
clus13<-pam(rna.sig.dists,13)$clustering

table(clus8,clus9)
table(clus8,clus13)
table(clus9,clus13)

cor.rna.sig<-cor(norm.M[rna.sig,],use="complete.obs")



################### code from hardin for adjrand func###############
adjrand<-function(clust1,clust2){

 

group.tab<-table(clust1,clust2)

col.sum<-apply(group.tab,2,sum)

row.sum<-apply(group.tab,1,sum)

tot.sum<-sum(group.tab)

 

num1<-sum(choose(group.tab,2))

num2<-sum(choose(row.sum,2))*sum(choose(col.sum,2))/choose(tot.sum,2)

den1<-0.5*(sum(choose(row.sum,2))+sum(choose(col.sum,2)))

 

arand<- (num1-num2)/(den1-num2)

 

arand}

#######################end code from hardin#############


adjrand(clus8,clus9)
adjrand(clus8,clus13)
adjrand(clus9,clus13)


cutree(rna.sig.clusts,k=9)
plot(rna.sig.clusts,
#     labels=paste(cutree(rna.sig.clusts,k=9),"-",rna.sig),
     labels=cutree(rna.sig.clusts,k=9),
     xlab=paste("dist =",this.dist,",  method =",this.method),
     main="Dendrogram for significant miRNA",
     sub="")


########## end project 9########




library(pamr)

single.arrays<-c(which(dat.targets$TissueType=="ERMS"),which(dat.targets$TissueType=="DDLPS"))

dat.norm2<-c()
dat.norm2$x<-norm.M[,-single.arrays]
rownames(dat.norm2$x)<-reporter.name
dat.norm2$y<-dat.targets$TissueType[-single.arrays]
dat.norm2$geneid<-as.character(c(1:nrow(dat.norm2$x)))
dat.norm2$genenames<-rownames(dat.norm2$x)

dat.norm<-pamr.knnimpute(dat.norm2)


pamr.menu(dat.norm)






dat.norm.NT2<-c()
dat.norm.NT2$x<-norm.M
rownames(dat.norm.NT2$x)<-reporter.name
dat.norm.NT2$y<-dat.targets$DiseaseState
dat.norm.NT2$geneid<-as.character(c(1:nrow(dat.norm.NT2$x)))
dat.norm.NT2$genenames<-rownames(dat.norm.NT2$x)

dat.norm.NT<-pamr.knnimpute(dat.norm.NT2)


pamr.menu(dat.norm.NT)



################################################################################################ end #