Too few genes being read into adata

[combine-lab]

Too few genes being read into adata

[BW15061999]

Hi,
I upgraded everything and rerun the whole workflow and now is fine!
Thank you!

Best

[rob-p]

Hi @BW15061999,

Could you provide some more information? What reference are you mapping against? How are you creating the splici reference and index? What does your transcript-to-gene mapping file look like? In general, the number of genes is determined by the reference you index and map against, and this, in turn, is determined by the annotation (e.g. GTF file) you use. For human, for example, most annotations have 20-30k genes.

Best,
Rob

[BW15061999]

Hi,
I followed the provided tutorial and used cellranger-mm10-2.1.0 as reference, but the raw data in the original paper is using GRCm38/mm10 as reference. I am not sure if it will influence the result.

Best

[rob-p]

So I can confirm that if this is the number of genes you are getting before filtering, that doesn't seem right. I just downloaded the refdata-cellranger-mm10-2.1.0.tar.gz reference, and build a splici reference from it using pyroe:

pyroe make-splici --filename-prefix splici-ref --flank-trim-length 5 --dedup-seqs refdata-cellranger-mm10-2.1.0/fasta/genome.fa refdata-cellranger-mm10-2.1.0/genes/genes.gtf 95 mm10-2.1.0-splici

the resulting reference has 28692 distinct genes, which is what should be present when you quantify against this index. Does your 3-column transcript-to-gene mapping file (the one you pass to the quant command) have more genes? I will note it seems a bit suspicious that 28692 / 3 = 9564, so it's almost like you are getting a 3rd of the genes — as if you didn't quantify in USA mode but the input is being read as if you did. Since your raw quantification output is in mtx format, can you just do head on your raw count matrix? How many rows and columns does it have?

--Rob

[BW15061999]

Hi,
the mtx file header is 3070 28692 2868327, and the quants_mat_cols.txt file has 28692 lines.
And here is how I load data:

frydir = "quant_res"
e2n_path = "geneid_to_name.txt"
adata = load_fry.load_fry(frydir, output_format = "velocity")
e2n = dict([ l.rstrip().split() for l in open(e2n_path).readlines()])
adata.var_names = [e2n[e] for e in adata.var_names]

and what it generate:

USA mode: True
Using pre-defined output format: velocity
Will populate output field X with sum of counts frorm ['S', 'A'].
Will combine ['S', 'A'] into output layer spliced.
Will combine ['U'] into output layer unspliced.

[rob-p]

Hi @BW15061999 ,

So this is exactly what I was positing. It looks like you have a matrix that is not in USA mode, but it's being loaded as if it is. How did you run the quant command? What transcript-to-gene file did you pass to it?

--Rob