Smart-seq vars in 0.25

[JBreunig]

Describe the issue
Trying to align paired Smart-seq data with the following command (or similar with alternate Human references), I get variable var types in my loom file but never the standard "gid" as with 10X data alignment to facilitate easy conversion back to the gene symbol with t2g.txt or other means.

EDIT: aligning a human 10X dataset (changing -x 10xv3) with 0.25 & same reference works fine and gives appropriate vars:
HuGBMv3_10X2.var
Out[3]:
Empty DataFrame
Columns: []
Index: [ENSG00000277400.1, ENSG00000274847.1, ENSG00000276256.1, ENSG00000278198.1

HuGBMv3_10X2
Out[4]:
AnnData object with n_obs × n_vars = 98545 × 58367

Using this reference or other (all human) references from your repositories or other with smart-seq fastqs, I get:
adata
Out[56]: AnnData object with n_obs × n_vars = 384 × 227368

adata.var
Out[55]:
Empty DataFrame
Columns: []
Index: [ENST00000003583.12, ENST00000003912.7, ENST00000008440.9,....

adata
Out[50]:
AnnData object with n_obs × n_vars = 384 × 188753

adata.var
Out[33]:
Empty DataFrame
Columns: []
Index: [ENST00000631435.1, ENST00000434970.2, ENST00000448914.1, ENST00000415118.1, ENST00000632684.1,...

adata
Out[27]: AnnData object with n_obs × n_vars = 384 × 845338

adata.var
Out[28]:
Empty DataFrame
Columns: []
Index: [ENSG00000277400.1.A14056, ENSG00000277400.1.A32841, ENSG00000277400.1.A35311,
What is the exact command that was run?
kb count -i /mnt/WD2TBunderside/WorkingHumanRefKBtools081520/Patcher096release/index.idx -g /mnt/WD2TBunderside/WorkingHumanRefKBtools081520/Patcher096release/t2g.txt -x SMARTSEQ -o BT1030e BT1030*.gz --loom --verbose -t 56 -m 196G

or are these outputs expected?

[Lioscro]

Hi, @JBreunig,
It looks like kallisto is outputting transcript IDs instead of gene IDs, probably because processing smartseq fastqs require a slightly different procedure than standard UMI-based single cell fastqs.

I agree that gene IDs are easier to work with. I'll take a look at how we can convert those to gene IDs instead.

Regarding the last adata.var you show, I'm not sure what is happening with the genes there. Perhaps you are using the index and t2g from this (https://github.com/pachterlab/kallisto-transcriptome-indices) repository? If so, my recommendation would be to manually generating the transcriptome index following the procedure described in the FAQ 1 in our tutorials because there is some additional columns in the t2g that kb uses, while the t2g in the repository do not include those.

[JBreunig]

Ok, that all makes sense. Thanks in advance for looking into converting to gene IDs!

[Lioscro]

Hi, @JBreunig,

I've pushed an update to the devel branch that should now output gene IDs instead of transcript IDs as the columns.
Could you try running and verify? Once I get the confirmation, we'll release a new version of kb with these changes.
You can install the devel branch with the following command:

pip install git+https://github.com/pachterlab/kb_python@devel --upgrade

[JBreunig]

Hi Joseph,
I think it's close but do you have a suggestion to remove the .XXXXXX places and condense duplicates from the format "ENSG00000277400.1.A14056" to "ENSG00000277400.1" as is shown in the tutorials and typical with droplet data?

Here are the outputs using the same reference validated with human 10X datasets:

BT1030f = sc.read('/mnt/1TBretrySamsung/HuEpendymoma/BT1030/BT1030rd/adata.loom', cache=True, validate=False)
writing an h5ad cache file to speedup reading next time

adata = BT1030f

t2g = pd.read_csv("/mnt/WD2TBunderside/WorkingHumanRefKBtools081520/t2g.txt", header=None, sep="\t", names=["tid", "gid", "gene"]) #Human

t2g = t2g.drop_duplicates(["gid", "gene"])
t2g = t2g.set_index("gid")
t2g.head()
Out[8]:
tid gene
gid
ENSG00000223972.5 ENST00000456328.2 DDX11L1
ENSG00000227232.5 ENST00000488147.1 WASH7P
ENSG00000278267.1 ENST00000619216.1 MIR6859-1
ENSG00000243485.5 ENST00000473358.1 MIR1302-2HG
ENSG00000284332.1 ENST00000607096.1 MIR1302-2

adata.var.head()#
Out[9]:
Empty DataFrame
Columns: []
Index: [ENSG00000277400.1.A14056, ENSG00000277400.1.A32841, ENSG00000277400.1.A35311, ENSG00000277400.1.A36796, ENSG00000277400.1.A44325]

adata
Out[10]: AnnData object with n_obs × n_vars = 384 × 845338

adata.var["Gene"] = adata.var.index.map(t2g["gene"])

adata.var.head()#
Out[12]:
Gene
ENSG00000277400.1.A14056 NaN
ENSG00000277400.1.A32841 NaN
ENSG00000277400.1.A35311 NaN
ENSG00000277400.1.A36796 NaN
ENSG00000277400.1.A44325 NaN

[Lioscro]

I actually have no idea where the AXXXX suffixes are coming from. I have never seen that happen on my test data. Are you using a cached version of the output matrix?

[JBreunig]

Just to be clear, cached from the loom?

When I look at the genes.txt and transcripts.txt in the output directory from today, they are there as well:
genes.txt
ENSG00000277400.1.A14056
ENSG00000277400.1.A32841
ENSG00000277400.1.A35311
ENSG00000277400.1.A36796

transcripts.txt
ENSG00000277400.1.A14056
ENSG00000277400.1.A32841
ENSG00000277400.1.A35311
ENSG00000277400.1.A36796

Retrying now with a new, rebuilt index...will post when done.

[Lioscro]

Yes, cached from the loom. I mentioned that because I noticed the cache=True option in your sc.read command, which may be reading a cached (older) version of the loom file.

I re-ran using the transcriptome provided here (https://github.com/pachterlab/kallisto-transcriptome-indices/releases) and am still unable to reproduce the outputs you are getting, and none of the AXXXX suffixes appear in the GTF, CDNA fasta, nor the t2g (so I'm not sure where these are coming from).

Please let me know how it looks when running with a manually built index.

[JBreunig]

Ok, seems to work using that index. Now I get a "gene_name" column that is easy to work with...thanks!!

adata
Out[6]:
AnnData object with n_obs × n_vars = 384 × 53292
var: 'gene_name'

adata.var.head()#
Out[8]:
gene_name
ENSG00000000003.14 TSPAN6
ENSG00000000005.6 TNMD
ENSG00000000419.12 DPM1
ENSG00000000457.14 SCYL3
ENSG00000000460.17 C1orf112