NMF on RORB data

#
import argparse
import math
from intern.remote.boss import BossRemote
from intern.resource.boss.resource import *
import configparser
#import grequests # for async requests, conflicts with requests somehow
import requests
import numpy as np
from numpy import genfromtxt
import h5py
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import shutil
import blosc
from IPython.core.debugger import set_trace
import sys
import os
import itertools
from functools import partial
from multiprocessing import Pool
from multiprocessing.dummy import Pool as ThreadPool
from multiprocessing import cpu_count
import csv
import datetime

def getCube(di):
    data = di['rem'].get_cutout(di['ch_rsc'], di['res'], di['xrng'],
          di['yrng'] ,di['zrng'])
    print(di['loc']) #DEBUG
    sys.stdout.flush() #DEBUG
    ### OUTPUT will be a numpy array in Z Y X order 
    #avg = np.mean(data)
    return(data)

CONFIG_FILE = "neurodata.conf"
COLL_NAME = "collman"
EXP_NAME = "M247514_Rorb_1_Site3Align2_LENS_Session1_CROP"
COORD_FRAME = COLL_NAME + "_" + EXP_NAME

CHAN_NAMES = ['DAPI1', 'DAPI2', 'DAPI3', 'GABA', 'GAD2', 'Gephyrin', 'GluN1',
                'MBP', 'PSD95', 'synapsin', 'TdTomato', 'VGlut1']


config = configparser.ConfigParser()
config.read(CONFIG_FILE)
TOKEN = config['Default']['token']
boss_url = ''.join( ( config['Default']['protocol'],'://',config['Default']['host'],'/v1/' ) )
#print(boss_url)
#'https://api.boss.neurodata.io/v1/'

#intern
rem = BossRemote(CONFIG_FILE)

chres = [ChannelResource(ch,COLL_NAME,EXP_NAME,'image',datatype='uint8') for ch in CHAN_NAMES]

data = []
for i in range(len(chres)):
  data.append(rem.get_cutout(chres[i],0, [100,300],[100,300],[12,25]))

colMaj = [data[i].flatten() for i in range(len(data))]

out = np.asarray(colMaj)

head = "DAPI1, DAPI2, DAPI3, GABA, GAD2, Gephyrin, GluN1, MBP, PSD95, synapsin, TdTomato, VGlut1"

np.savetxt("test1Cube.csv", np.transpose(out), delimiter = ',', fmt = "%d", header = head)

0.1 A test case

Using a chunk from the RORB data at x = [100,300], y = [100,300], z = [12,25]

cu <- fread("test1Cube.csv")
colnames(cu) <- cnames
ind <- apply(cu, 1, function(x) all(x ==0))

tmp <- as.numeric(cu[,8][[1]])
b1 <- array(tmp, dim = c(200,200,13))

plot(raster::brick(b1, transpose = TRUE), col = gray.colors(255))

tmp <- matrix(NA, 520000, 3)
tmp[ind,] <- 0

nam <- cu[!ind,]
tmp <- matrix(NA, nrow(cu), 3)
tmp[ind,] <- 0

N <- foreach(x = c(1:ncol(nam))) %dopar% {
    nmf(nam, rank = x)
}

err <- 
foreach(x = 1:ncol(nam), .combine = 'c') %do% {
  norm(as.matrix(nam) - (N[[x]]@fit@W %*% N[[x]]@fit@H), "F")^2 / 
      norm(as.matrix(nam), "F")^2
}

getElbows(err)

## [1] 4 6 8

for(i in 2:5){
  if(i < 3){
    tmp[!ind, 1:i] <- N[[i]]@fit@W
  } else {
    tmp[!ind,1:3] <- N[[i]]@fit@W[,1:3]
  }
  am <- array(tmp, dim = c(200,200,i))
  show(plotRGB(raster::brick(am), r = 1, g = 2, b = 3))
  show(plot(raster::brick(am, transpose = TRUE)))
  show(heatmap(N[[i]]@fit@H, col = gray.colors(255)))
}

## Warning in .local(x, ...): layer was changed to 2

## Warning in couldBeLonLat(x): CRS is NA. Assuming it is longitude/latitude

## NULL

## NULL

## $rowInd
## [1] 2 1
## 
## $colInd
##  [1] 12 11  4 10  1  6  2  7  3  9  8  5
## 
## $Rowv
## NULL
## 
## $Colv
## NULL

## Warning in couldBeLonLat(x): CRS is NA. Assuming it is longitude/latitude

## NULL

## NULL

## $rowInd
## [1] 1 2 3
## 
## $colInd
##  [1] 11  4 10  8  5 12  1  6  9  3  2  7
## 
## $Rowv
## NULL
## 
## $Colv
## NULL

## Warning in couldBeLonLat(x): CRS is NA. Assuming it is longitude/latitude

## NULL

## NULL

## $rowInd
## [1] 1 4 3 2
## 
## $colInd
##  [1] 11  5 10 12  8  4  7  1  6  9  2  3
## 
## $Rowv
## NULL
## 
## $Colv
## NULL

## Warning in couldBeLonLat(x): CRS is NA. Assuming it is longitude/latitude

## NULL

## NULL

## $rowInd
## [1] 1 2 5 4 3
## 
## $colInd
##  [1] 11 12  8  4 10  5  9  7  6  1  2  3
## 
## $Rowv
## NULL
## 
## $Colv
## NULL