Examples depicting the visualization and analysis of Retinal neuron data from mouse brain

In [1]:

    

library('curl')
library('R.matlab')
library('gdata')


    

    




Registered S3 method overwritten by 'R.oo':
  method        from       
  throw.default R.methodsS3
R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.

Attaching package: ‘R.matlab’

The following objects are masked from ‘package:base’:

    getOption, isOpen

gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.

gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.

Attaching package: ‘gdata’

The following object is masked from ‘package:stats’:

    nobs

The following object is masked from ‘package:utils’:

    object.size

The following object is masked from ‘package:base’:

    startsWith

In [2]:

    

# Set up URLs for data download..
download_urls <- c('http://neuro.rzg.mpg.de/download/Helmstaedter_et_al_Nature_2013_skeletons_contacts_matrices.mat',
                   'https://media.nature.com/original/nature-assets/nature/journal/v500/n7461/extref/nature12346-s3.zip')


    

In [3]:

    

# Set up Folder name insides 'natverse_examples'..
dataset_name <- 'natverse_mouseconnectome'


    

In [4]:

    

download_filename <- basename(download_urls)
if (is.null(getOption(dataset_name))){
    message("Setting options")
    options(natverse_examples=rappdirs::user_data_dir('R/natverse_examples'))
  }


    

    




Setting options

In [5]:

    

dataset_path <- file.path(getOption('natverse_examples'),dataset_name)
message("Dataset path:", dataset_path)


    

    




Dataset path:/Users/sri/Library/Application Support/R/natverse_examples/natverse_mouseconnectome

In [6]:

    

if (!dir.exists(dataset_path)){
    message("Creating folder: ", basename(dataset_path))
    dir.create(dataset_path, recursive = TRUE, showWarnings = FALSE)
  }


    

In [7]:

    

downloaded_dataset <- ''
if(length(download_urls))
    message("Downloading files to: ", dataset_path)
for (download_fileidx in 1:length(download_urls)){
        localfile <- file.path(dataset_path, download_filename[download_fileidx])
        downloaded_dataset[download_fileidx] <- localfile
        if(file.exists(localfile)) next
        message("Processing URL: ", download_urls[download_fileidx], " (file ", download_fileidx, "/", length(download_urls), ")")
        curl::curl_download(download_urls[download_fileidx], localfile, quiet=FALSE)
      }
message("Downloads done")


    

    




Downloading files to: /Users/sri/Library/Application Support/R/natverse_examples/natverse_mouseconnectome
Downloads done

In [8]:

    

message(paste(c("Dataset downloaded at: ", downloaded_dataset), collapse="\n"))


    

    




Dataset downloaded at: 
/Users/sri/Library/Application Support/R/natverse_examples/natverse_mouseconnectome/Helmstaedter_et_al_Nature_2013_skeletons_contacts_matrices.mat
/Users/sri/Library/Application Support/R/natverse_examples/natverse_mouseconnectome/nature12346-s3.zip

In [9]:

    

raw_data <- R.matlab::readMat(downloaded_dataset[1])


    

In [10]:

    

names(raw_data)


    

    






	
'kn.e2006.ALLSKELETONS.FINAL2012'
	
'kn.e2006.ALLSKELETONS.FINAL2012.cellIDs'
	
'kn.e2006.ALLSKELETONS.FINAL2012.globalTypeIDs.REDOMAR2013'
	
'kn.e2006.ALLSKELETONS.FINAL2012.cellIDs.pure.forSomata'
	
'kn.e2006.ALLSKELETONS.FINAL2012.allSomata'
	
'kn.e2006.ALLSKELETONS.FINAL2012.sortOrderTypes.byDepth.MAR2013'
	
'kn.e2006.ALLSKELETONS.FINAL2012.cellIDs.sortedByType.MAR2013'
	
'kn.e2006.ALLSKELETONS.FINAL2012.cellIDs.sortedByTyp.MAR2013.typ'
	
'kn.allContactData.Interfaces.duplCorr.sizeCorr.APR2013'
	
'kn.contactAreaMx.APR2013'
	
'kn.contactAreaMx.symm.APR2013'
	
'kn.contactAreaMx.forDisplay.APR2013'
	
'kn.e2006.typeMx.REDOAPR2013'
	
'kn.e2006.typeMx.REDOAPR2013.selfNorm'

In [11]:

    

skeleton_metadata <- raw_data[2:8]
saveRDS(skeleton_metadata, file=file.path(dataset_path,'skeleton_metadata.rds'))


    

In [12]:

    

skall <- raw_data$kn.e2006.ALLSKELETONS.FINAL2012


    

In [13]:

    

rownames(skall[[1]][[1]])


    

    






	
'parameters'
	
'commentsString'
	
'nodes'
	
'nodesNumDataAll'
	
'edges'
	
'thingID'
	
'branchpointsString'
	
'branchpoints'
	
'sanityReport'
	
'edgeSel'

In [14]:

    

#' Parse a single skeleton object cleaning up the rather messy structure that comes out of readMat
parse.moritz.skel<-function(x, simple_fields_only=TRUE){
  # delist
  x=x[[1]]
  stopifnot(inherits(x,'array'))
  
  vars=rownames(x)
  stopifnot(all(c('nodes','edges') %in% vars))
  
  #simplevars=intersect(c('nodes','edges', 'edgeSel'), vars)
  simplevars=intersect(c('nodes','edges'), vars)
  r=sapply(simplevars, function(v) x[v,,][[1]], simplify=FALSE)
  
  othervars=setdiff(vars, simplevars)
  process_var<-function(y){
    if(inherits(y,'list')) y=y[[1]]
    if(inherits(y,'array')) apply(y, 1, unlist) else {
      if(is.numeric(y)) drop(y) else y
    }
  }
  if(simple_fields_only)
    structure(r, class=c('skel','list'))
  else {
    r2=sapply(othervars, function(v) process_var(x[v,,]), simplify = FALSE)
    structure(c(r, r2), class=c('skel','list'))
  }
}


    

In [15]:

    

# Convert all the neurons to intermediate skeleton format (here just a neuronlist in 'nat')
skallp=nat::nlapply(skall, parse.moritz.skel, .progress='text')
# give the neurons names
names(skallp)=sprintf("sk%04d", seq_along(skallp))


    

    




Registered S3 method overwritten by 'nat':
  method             from
  as.mesh3d.ashape3d rgl 






    




  |======================================================================| 100%

In [16]:

    

skallp_unique <- skallp


    

In [17]:

    

attr(skallp_unique,'df')=data.frame(
  row.names=names(skallp_unique),
  skid=seq.int(skallp_unique),
  cellid=skeleton_metadata$kn.e2006.ALLSKELETONS.FINAL2012.cellIDs[1,],
  typeid=skeleton_metadata$kn.e2006.ALLSKELETONS.FINAL2012.globalTypeIDs.REDOMAR2013[1,],
  cellid.soma=skeleton_metadata$kn.e2006.ALLSKELETONS.FINAL2012.cellIDs.pure.forSomata[1,]
  )


    

In [18]:

    

head(attr(skallp_unique,'df'))


    

    






A data.frame: 6 × 4

	
skid
cellid
typeid
cellid.soma
	
<int>
<dbl>
<dbl>
<dbl>


	
sk0001
1
1
10204
1
	
sk0002
2
1
10204
1
	
sk0003
3
1
10204
1
	
sk0004
4
1
10204
1
	
sk0005
5
1
10204
1
	
sk0006
6
1
10204
1

In [19]:

    

attr(skallp_unique,'df')$nedges=sapply(skallp_unique,function(x) nrow(x$edges))


    

In [20]:

    

message("Number of unique cells(neurons): ", length(unique(attr(skallp_unique,'df')$cellid)))


    

    




Number of unique cells(neurons): 1185

In [21]:

    

message("Number of tracings: ", length(skallp_unique))


    

    




Number of tracings: 4188

In [22]:

    

sk.uniq_temp=skallp_unique[by(attr(skallp_unique,'df'),attr(skallp_unique,'df')$cellid,function(x) x$skid[which.max(x$nedges)])]


    

In [23]:

    

message("Number of tracings used: ", length(sk.uniq_temp))


    

    




Number of tracings used: 1185

In [24]:

    

message("Reading SI 4 spreadsheet sheet 3 ...")
utils::unzip(downloaded_dataset[2],exdir = file.path(dataset_path, sub('\\.zip$', '', basename(downloaded_dataset[2]))))
SI4.s3=gdata::read.xls(file.path(dataset_path, sub('\\.zip$', '', basename(downloaded_dataset[2])),'Helmstaedter_et_al_SUPPLinformation4.xlsx'), sheet=3)
message("Reading done")


    

    




Reading SI 4 spreadsheet sheet 3 ...
Reading done

In [25]:

    

head(SI4.s3)


    

    






A data.frame: 6 × 4

	
ID.in.cell2cell.Mx
cell.ID.in.skeleton.db
internal.type.ID
sorted.type.ID..as.in.Type.Mx.
	
<int>
<int>
<int>
<int>


	
1
77
201
1
	
2
75
201
1
	
3
16
231
2
	
4
47
231
2
	
5
88
231
2
	
6
23
202
3

In [26]:

    

saveRDS(SI4.s3,file=file.path(dataset_path,"SI4.s3.rds"))


    

In [27]:

    

# assign sorted typed id as used in the paper
attr(sk.uniq_temp,'df')$stypeid=0
attr(sk.uniq_temp,'df')$stypeid[match(SI4.s3$cell.ID.in.skeleton.db,attr(sk.uniq_temp,'df')$cellid)]=SI4.s3$sorted.type.ID..as.in.Type.Mx.


    

In [28]:

    

# Make classes for 
# 1:12 ganglion
# 13:57 amacrine
# 58:71 bipolar
df=attr(sk.uniq_temp,'df')
df$class=''
df<-within(df,class[stypeid%in%1:12]<-'ganglion')
df<-within(df,class[stypeid%in%13:57]<-'amacrine')
# narow 13-24
# wide: 25, 28, 30–32, 37, 39, 41, 47, 53 and 57
# medium: types 26, 27, 29, 33-36, 38, 40, 42–46, 48–52 and 54-57
# typo for that last 57?

df<-within(df,class[stypeid%in%58:71]<-'bipolar')
# I assume that these must be Muller glial cells
df<-within(df,class[stypeid%in%77]<-'glial')
df$class=factor(df$class)

df$subclass=''
df$subclass[df$stypeid%in%13:24]='narrow field'
df$subclass[df$stypeid%in%c(25, 28, 30:32, 37, 39, 41, 47, 53, 57)]='wide field'
df$subclass[df$stypeid%in%c(26, 27, 29, 33:36, 38, 40, 42:46, 48:52,54:57)]='medium field'
df$subclass=factor(df$subclass)

df$ntype=''
df$ntype[df$stypeid%in%c(33,51)]='starburst amacrine'
attr(sk.uniq_temp,'df')=df


    

In [29]:

    

saveRDS(sk.uniq_temp, file=file.path(dataset_path,'sk.uniq.rds'))


    

In [30]:

    

file.path(dataset_path,'sk.uniq.rds')


    

    





'/Users/sri/Library/Application Support/R/natverse_examples/natverse_mouseconnectome/sk.uniq.rds'

In [ ]:

    




    

In [33]:

    

load('/Users/sri/Documents/R/dev/nat.examples/03-helmstaedter2013/sk.uniq.rda')


    

In [34]:

    

str(sk.uniq[[1]])


    

    




List of 2
 $ nodes: num [1:2316, 1:4] 124658 121126 121242 121390 121605 ...
 $ edges: num [1:2295, 1:2] 3 4 5 6 7 8 9 10 11 12 ...
 - attr(*, "header")=List of 3
  ..$ description: chr "MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Sat Aug 17 09:49:30 2013                                                 "
  ..$ version    : chr "5"
  ..$ endian     : chr "little"
 - attr(*, "class")= chr [1:2] "skel" "list"

In [35]:

    

str(sk.uniq_temp[[1]])


    

    




List of 2
 $ nodes: num [1:2316, 1:4] 124658 121126 121242 121390 121605 ...
 $ edges: num [1:2295, 1:2] 3 4 5 6 7 8 9 10 11 12 ...
 - attr(*, "class")= chr [1:2] "skel" "list"

In [36]:

    

head(attr(sk.uniq_temp,'df'))


    

    






A data.frame: 6 × 9

	
skid
cellid
typeid
cellid.soma
nedges
stypeid
class
subclass
ntype
	
<int>
<dbl>
<dbl>
<dbl>
<int>
<dbl>
<fct>
<fct>
<chr>


	
sk0003
 3
1
10204
1
2295
57
amacrine
medium field

	
sk0012
12
2
10990
2
1259
72
        
            

	
sk0018
18
3
10999
3
  57
 0
        
            

	
sk0024
24
4
  181
4
4123
43
amacrine
medium field

	
sk0027
27
5
  216
5
3865
 9
ganglion
            

	
sk0035
35
6
  214
6
4352
10
ganglion
            

In [37]:

    

head(attr(sk.uniq,'df'))


    

    






A data.frame: 6 × 9

	
skid
cellid
typeid
cellid.soma
nedges
stypeid
class
subclass
ntype
	
<int>
<dbl>
<dbl>
<dbl>
<int>
<dbl>
<fct>
<fct>
<chr>


	
sk0003
 3
1
10204
1
2295
57
amacrine
medium field

	
sk0012
12
2
10990
2
1259
72
        
            

	
sk0018
18
3
10999
3
  57
 0
        
            

	
sk0024
24
4
  181
4
4123
43
amacrine
medium field

	
sk0027
27
5
  216
5
3865
 9
ganglion
            

	
sk0035
35
6
  214
6
4352
10
ganglion
            

In [38]:

    

all.equal(attr(sk.uniq,'df'),attr(sk.uniq_temp,'df'))


    

    





TRUE

In [39]:

    

all.equal(sk.uniq,sk.uniq_temp, check.attributes = F)


    

    





TRUE

In [ ]:

    




    

In [40]:

    

#install.packages("/Users/sri/Documents/R/dev/nat",repos = NULL, 
#type = "source",force=T)


    

In [41]:

    

library(nat)


    

    




Loading required package: rgl

Attaching package: ‘nat’

The following object is masked from ‘package:gdata’:

    resample

The following objects are masked from ‘package:base’:

    intersect, setdiff, union

In [42]:

    

#' Convert Helmstaedter's matlab skel format into nat::neuron objects
#' 
#' @description skel objects are my direct R translation of the matlab data
#'   provided by Briggman, Helmstaedter and Denk.
#' @param x A skel format neuron to convert
#' @param ... arguments passed to as.neuron
#' @return An object of class \code{neuron}
#' @seealso \code{\link{as.neuron}}
as.neuron.skel<-function(x, ...) {
  as.neuron(as.ngraph(x), ...)
}


    

In [43]:

    

#' Convert Helmstaedter's matlab skel format into nat::ngraph objects
#' 
#' @details \code{ngraph} objects are thin wrappers for \code{igraph::graph} 
#'   objects.
#'   
#'   This function always removes self edges (i.e. when a node is connected to 
#'   itself), which seem to be used as a placeholder in the Helmstaedter dataset
#'   when there are no valid edges.
#'   
#'   Isolated nodes are also removed by default (i.e. the nodes that are not
#'   connected by any edges.)
#' @param remove.isolated Whether or not to remove isolated nodes from the graph
#'   (see Details).
#' @inheritParams as.neuron.skel
#' @seealso \code{\link[nat]{ngraph}}, \code{\link{as.neuron.skel}}
as.ngraph.skel<-function(x, remove.isolated=TRUE, ...) {
  self_edges=x$edges[,1]==x$edges[,2]
  if(sum(self_edges)>0){
    if(sum(self_edges)==nrow(x$edges)){
      # there are only self edges - make a dummy empty edge matrix
      x$edges=matrix(nrow=0, ncol=2)
    } else {
      x$edges[!self_edges, , drop=FALSE]
    }
  }
  ng=ngraph(x$edges, vertexnames = seq_len(nrow(x$nodes)), xyz=x$nodes, ...)
  if(remove.isolated){
    isolated_vertices=igraph::V(ng)[igraph::degree(ng)==0]
    g=igraph::delete.vertices(graph=ng,isolated_vertices)
    ng=as.ngraph(g)
  }
}


    

In [44]:

    

options(nat.progress="none")


    

In [45]:

    

options(warn=-1)
skn=nlapply(sk.uniq_temp, as.neuron, OmitFailures = T, Verbose = FALSE)
options(warn=0)


    

In [46]:

    

class(skn[[1]])


    

    






	
'neuron'
	
'list'

In [47]:

    

#options(nat.plotengine = 'plotly')
#getOption('nat.plotengine')


    

In [48]:

    

library(htmlwidgets)
library(plotly)
library(IRdisplay)


    

    




Loading required package: ggplot2

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

In [54]:

    

clearplotlyscene()
tempval <- plot3d(skn[1:5], plotengine ='plotly', soma = TRUE)
#tempval$plotlyscenehandle


    

In [55]:

    

htmlwidgets::saveWidget(plotly::as_widget(tempval$plotlyscenehandle), 'helmstaedter2013_01.html')
IRdisplay::display_html(paste0('<iframe src="','helmstaedter2013_01.html" width=100%, height=500> frameborder="0" </iframe>'))


    

    





 frameborder="0" 

In [60]:

    

clearplotlyscene()
tempval <- plot3d(skn,ntype=='starburst amacrine',col=stypeid, plotengine = 'plotly', soma = TRUE, alpha = 0.5)
#tempval$plotlyscenehandle


    

In [61]:

    

htmlwidgets::saveWidget(plotly::as_widget(tempval$plotlyscenehandle), 'helmstaedter2013_02.html')
IRdisplay::display_html(paste0('<iframe src="','helmstaedter2013_02.html" width=100%, height=500> frameborder="0" </iframe>'))


    

    





 frameborder="0" 

In [ ]:

    




    

In [62]:

    

sessionInfo()


    

    






R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] IRdisplay_0.7.0   plotly_4.9.0.9000 ggplot2_3.2.1     htmlwidgets_1.3  
 [5] nat_1.9.1.9000    rgl_0.100.19      Matrix_1.2-17     gdata_2.18.0     
 [9] R.matlab_3.6.2    curl_4.0         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2              lattice_0.20-38         tidyr_0.8.3            
 [4] gtools_3.8.1            assertthat_0.2.1        zeallot_0.1.0          
 [7] digest_0.6.20           mime_0.7                R6_2.4.0               
[10] plyr_1.8.4              repr_1.0.1              backports_1.1.4        
[13] evaluate_0.14           httr_1.4.0              pillar_1.4.2           
[16] rlang_0.4.0             lazyeval_0.2.2          uuid_0.1-2             
[19] data.table_1.12.2       miniUI_0.1.1.1          R.utils_2.9.0          
[22] R.oo_1.22.0             webshot_0.5.1           igraph_1.2.4.1         
[25] munsell_0.5.0           shiny_1.3.2             compiler_3.6.0         
[28] httpuv_1.5.1            xfun_0.8.1              pkgconfig_2.0.2        
[31] base64enc_0.1-3         htmltools_0.3.6         tidyselect_0.2.5       
[34] tibble_2.1.3            codetools_0.2-16        viridisLite_0.3.0      
[37] withr_2.1.2             crayon_1.3.4            dplyr_0.8.1            
[40] later_0.8.0             nat.utils_0.5.1         R.methodsS3_1.7.1      
[43] rappdirs_0.3.1          grid_3.6.0              jsonlite_1.6           
[46] xtable_1.8-4            gtable_0.3.0            magrittr_1.5           
[49] scales_1.0.0            nabor_0.5.0             promises_1.0.1         
[52] vctrs_0.2.0             IRkernel_1.0.2.9000     tools_3.6.0            
[55] manipulateWidget_0.10.0 glue_1.3.1              purrr_0.3.2            
[58] crosstalk_1.0.0         yaml_2.2.0              colorspace_1.4-1       
[61] filehash_2.4-2          pbdZMQ_0.3-3            knitr_1.23             

In [ ]: