How to transform neurons to a list of dotprops?

89 views
Skip to first unread message

Jingpeng Wu

unread,
Oct 10, 2018, 10:35:13 PM10/10/18
to nat-user
Hi,

I am trying to read a serials of swc files and do hierarchical clustering. 

here is my code example. It works to create a list of dotprops, but can not be used to perform hierarchical clustering. 
The data structure is probably different with the example of kcs20.

```
neuronList = c()
for (neuronId in neuronIdList) {
  fileName = paste("../01_data/0703/postprocessed/atlas_space/swc/", neuronId, ".swc", sep="")
  neuron = read.neuron(fileName)
  neuron = dotprops(neuron)
  neuronList <- c(neuronList, neuron)
}
```

Could you give me some hints of how to do it correctly?

Jingpeng

Gregory Jefferis

unread,
Oct 11, 2018, 4:41:01 AM10/11/18
to Jingpeng Wu, nat-user
Dear Jingpeng,

Thanks for your interest in the nat package/NBLAST.

To read in multiple neurons, there is a function called http://jefferis.github.io/nat/reference/read.neurons.html which will simplify things considerably. If you have a suggestion to make this easier to find in the documentation, please raise a GitHub issue. If you have not looked at http://jefferis.github.io/nat/reference/index.html, you may find this helpful.

Here's some sample code (untested) based on yours.

# assuming this directory contains your neurons
# you probably don't need the pattern if there are only swc files in the directory
mynlist=read.neurons("../01_data/0703/postprocessed/atlas_space/swc/", pattern='swc$')
# See http://jefferis.github.io/nat/reference/dotprops.html for the meaning of the additional arguments
# resample=1 => 1µm resampling and assumes your neurons are calibrated in µm
mydotprops=dotprops(mynlist, resample=1, k=5)

library(nat.nblast)
aba=nblast_allbyall(mydotprops)
hc=nhclust(scoremat=aba)
plot(hc)
plot3d(hc, k=3, db=mynlist)

All the best,

Greg.

On 11 Oct 2018, at 03:35, Jingpeng Wu <jingp...@gmail.com> wrote:

neuronList = c()
for (neuronId in neuronIdList) {
  fileName = paste("../01_data/0703/postprocessed/atlas_space/swc/", neuronId, ".swc", sep="")
  neuron = read.neuron(fileName)
  neuron = dotprops(neuron)
  neuronList <- c(neuronList, neuron)
}


--
Gregory Jefferis, PhD
Division of Neurobiology
MRC Laboratory of Molecular Biology
Francis Crick Avenue
Cambridge Biomedical Campus
Cambridge, CB2 OQH, UK




Jingpeng Wu

unread,
Oct 11, 2018, 2:57:59 PM10/11/18
to nat-user
Dear Gregory,

Thanks for response. It is really helpful. I have read the neuron list and clustered the cells.

BTW, the method of nhclust should be specified in my environment of latest R 3.5. the default of "ward" is not working, it should be "ward.D" or "ward.D2". I am using "ward.D" for my test.

I am selecting some neurons from the directory, and my node coordinates are in nm rather than micron.
so my code is like this:
```
fileNameList = c()
for (neuronId in neuronIdList) {
  fileName = paste("../01_data/0703/postprocessed/atlas_space/swc/", neuronId, ".swc", sep="")
  fileNameList <-c(fileNameList, fileName)
}

neurons = read.neurons(fileNameList)
mydotprops = dotprops(neurons, resmaple=1000, k=20)
# plot3d(mydotprops)

aba=nblast_allbyall(mydotprops)
hc=nhclust(scoremat=aba, method="ward.D")
plot(hc)
plot3d(hc, k=3, db=mynlist)
```

This code works and gives me classes of neurons, but he result is not as expected. The clusters looks random. I am probably doing something wrong. Do you have any idea? Any hints would be appreciated. 

Best,
Jingpeng

Jingpeng Wu

unread,
Oct 11, 2018, 3:02:35 PM10/11/18
to nat-user
According to the documentation, I might need to set my own sd using algorithm version 1. I am working on a zebra fish dataset from EM. 

Gregory Jefferis

unread,
Oct 11, 2018, 4:37:43 PM10/11/18
to Jingpeng Wu, nat-user
I’m glad that was helpful. As I mentioned in my last message you need your units to be in microns. This is very important for nblast. If your neurons are calibrated in nm then you need to convert ie divide by 1000. You can just take the object and divide by 1e3. 

mydotprops = dotprops(neurons/1e3, resample=1, k=5)

Also I suggest you construct your file list more like this:

fileNameList = paste0("../01_data/0703/postprocessed/atlas_space/swc/",neuronIdList, ".swc")

Most base functions in R are vectorised and recycle arguments so that you can avoid for loops. 

Finally there is a warning about the clustering method but it is not an error. If you do specify the clustering method, use the same one as it would use (I can’t remember which of the two that is). 

Best,

Greg. 

Sent from my iPhone
--
You received this message because you are subscribed to the Google Groups "nat-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nat-user+u...@googlegroups.com.
Visit this group at https://groups.google.com/group/nat-user.
To view this discussion on the web, visit https://groups.google.com/d/msgid/nat-user/92a0da43-6c3c-48e5-a317-cebcac8242ec%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Jingpeng Wu

unread,
Oct 12, 2018, 2:18:46 PM10/12/18
to nat-user
Dear Greg,

your suggestions are super helpful! I got reasonable results after transforming to micron unit. 
The vectorized construction is much faster.

Best,
Jingpeng

Gregory Jefferis

unread,
Oct 13, 2018, 3:19:04 AM10/13/18
to Jingpeng Wu, nat-user
Great! Happy to try and answer additional queries or to give comments offline if you want to get input on your results. All the best, Greg.

Sent from my iPhone
Reply all
Reply to author
Forward
0 new messages