Models aren't attempted for >90% of genes?

[danieljrichard]

Hello,
I've been trying to get SCARlink to work on my datasets. In the Seurat object, I manually set the VariableGenes slot to contain a set of ~1200 marker genes. Preprocessing seems to work, as hvg.txt contains the list of 1200 genes. However, running scarlink itself, the log only appears to attempt ~35 genes, of which only 4 were modelled (it appears for sparsity reasons).

Is there perhaps I'm doing something wrong? Or some silent filter in scarlink that tosses out the majority of variable genes? Any help would be greatly appreciated.

Daniel

[snehamitra]

We don't filter out any genes besides the sparsity threshold. Could you check the number of genes in your coassay_matrix.h5 file?

import h5py
f = h5py.File("output_dir/coassay_matrix.h5", 'r')
print(f.keys())
print(len(f.keys()))
f.close()

[Jia340]

I am not an author, just wonder if you set -np and -p but actually ran the job on only one core? It could happen if your server doesn't have a job scheduler. It appears that the scarlink code split genes across cores using the following code:
gene_names = [rm.gene_names[i] for i in range(p_ix, len(rm.gene_names), total_procs)]

My understanding is that when p_ix (-p) == 1 and total_procs (-np) ==100, genes 1, 101, 201, 301, ... will be run on the same core. With a job scheduler, the p_ix (id of process) will automatically be populated so all genes will be modeled. However, without a job scheduler, only 1/np genes will be run.

The easiest way to fix it is to remove -np and -p. However, it can take a very long time unless you reduce your number of gene upfront. Please correct me if my understanding is wrong - I am also wondering how -np and -p could be used on a server without job scheduler.