Gaussian Mixture Model to determine co-localization

library(rstan)

gaussianMM<-'
data {  
  int K;
  int N;
  real x[N];
}

parameters {
  simplex[K] theta;
  real mu[K];
}
model {
  real logps[K];
  for (k in 1:K) {
    mu[k] ~ normal( 0, 10 );
  }
  for (n in 1:N) {
    for (k in 1:K) {
      logps[k] <- log(theta[k]) + normal_log( x[n], mu[k], .3);
    }
    lp__ <- lp__ + log_sum_exp( logps );
  }
}'



quartz(width= 47.658/10, height= 47.658/10)
par(mfrow=c(2,2), yaxs ='i', xaxs='i', mar=c(4,4,1,1))
alpha<-30
limit<-2
ch.size<-0.5
color<-c(rep(paste('#008837', alpha, sep=''),nrow(dataset)), rep(paste('#7B3294', alpha, sep=''), nrow(datasetPvalb)), rep(paste('#CA0020', alpha, sep=''), nrow(datasetNPY) ) )
x<-sqrt(mydata[,1])
y<-sqrt(mydata[,2])
rnd.index<-sample(1:length(x))
plot(x[rnd.index],y[rnd.index], pch=16, cex=ch.size, xlim=c(0.5,2),ylim=c(0.5,2), col=color[rnd.index], xlab='Cy3 Channel (NPY)', ylab='FITC Channel (Lhx6)', asp=1, axes=F)

right.quad<-round(length(which((x > 16000)&(y>30000)))/nrow(mydata)*100, 2)

axis(1, at=c(0,limit/2, limit) )
axis(2, at=c(0,limit/2, limit) , las=1)
abline(v=30000, lwd=4)
abline(h=16000, lwd=4)
abline(v= 30000, lwd=2, col='orange')
abline(h= 16000, lwd=2, col='orange')
text( (limit-30000)*0.5+ 30000, (limit-16000)*0.5+16000, paste('Lhx6 & NPY\n',right.quad,'%'), cex=1)

par(mar=c(4,2,1,1), yaxs='i', xaxs='i',)
x<-sqrt(mydata[,3])
y<-sqrt(mydata[,1])
plot(x[rnd.index],y[rnd.index], pch=16, cex=ch.size, xlim=c(0.5,2),ylim=c(0.5,2), col= color[rnd.index], xlab='', ylab='',  asp=1, axes=F)
right.quad<-round(length(which((x > 16000)&(y> 16000)))/nrow(mydata)*100, 2)

axis(1, at=c(0,limit/2, limit) )
axis(2, at=c(0,limit/2, limit) , las=1)
abline(v= 16000, lwd=4)
abline(h= 16000, lwd=4)
abline(v= 16000, lwd=2, col='orange')
abline(h= 16000, lwd=2, col='orange')
text( (limit-16000)*0.5+ 16000, (limit-16000)*0.5+16000, paste('Lhx6 & Pvalb\n',right.quad,'%'), cex=1)


par(mar=c(4,2,1,1), yaxs ='i', xaxs='i',)
plot(0, col=0, axes=F, ylab='', xlab='')
x<-sqrt(mydata[,3])
y<-sqrt(mydata[,2])
plot(x[rnd.index],y[rnd.index], pch=16, cex=ch.size, xlim=c(0.5,2),ylim=c(0.5,2), col=color[rnd.index], xlab='Cy5 Channel (Pvalb)', ylab='Cy3 Channel (NPY)',  asp=1, axes=F)
x.max<-10.09543
y.max<-12.29063
right.quad<-round(length(which((x > x.max)&(y> y.max)))/nrow(mydata)*100, 2)

axis(1, at=c(0,limit/2, limit) )
axis(2, at=c(0,limit/2, limit) , las=1)
abline(v= x.max, lwd=4)
abline(h= y.max, lwd=4)
abline(v= x.max, lwd=2, col='orange')
abline(h= y.max, lwd=2, col='orange')


# Ward Hierarchical Clustering
par(mfrow=c(1,2))
d <- dist(sqrt(mydata), method = "euclidean") # distance matrix
fit <- hclust(d, method="ward") 
plot(fit,  col='gray', main='hierarchical clustering of fluorescence') # display dendogram
groups <- cutree(fit, k=7) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters 
rect.hclust(fit, k=7, border= unique(groups))
plot(mydata2$ML, mydata2$DV, col=groups, pch=16, cex=0.3, ylim=c(-8,2), xlim=c(-5,9), xlab='Mediolateral', ylab='Dorsoventral')
legend('topright', paste('Cluster', c(1:8)), pch=16, col=unique(groups), cex=0.85)

library(mclust)
fit2 <- Mclust(mydata)
plot(fit2) # plot results 
summary(fit2) # display the best model

plot(mydata2$ML, mydata2$DV, col=fit2$classification, pch=16, cex=0.3)

text( (16-x.max)*0.5+ x.max, (16-y.max)*0.5+ y.max, paste('Pvalb & NPY\n',right.quad,'%'), cex=1)

mydata[which(is.na(mydata), arr.ind=TRUE)]<-0

km <- kmeans(sqrt(mydata),4,10000)
# run principle component analysis
pc<-prcomp(sqrt(mydata))
# plot dots
plot(pc$x[,1], pc$x[,2],col=km$cluster,pch=16)

tsne_data <- tsne(cbind(sqrt(mydata), mydata2$DV, mydata2$ML), k=2, max_iter=500, epoch=500)
plot(tsne_data[,1], tsne_data[,2], col=km$cluster, pch=16)