Jupyter Notebook Binder

Bulk RNA-seq#

Note

More comprehensive examples are provided for these data types:

Setup#

!lamin init --storage test-bulkrna --schema bionty
Hide code cell output
✅ saved: User(uid='DzTjkKse', handle='testuser1', name='Test User1', updated_at=2024-03-04 17:09:17 UTC)
✅ saved: Storage(uid='pfpl5ptM', root='/home/runner/work/lamin-usecases/lamin-usecases/docs/test-bulkrna', type='local', updated_at=2024-03-04 17:09:17 UTC, created_by_id=1)
💡 loaded instance: testuser1/test-bulkrna
💡 did not register local instance on lamin.ai
import lamindb as ln
from pathlib import Path
import bionty as bt
import pandas as pd
import anndata as ad
💡 lamindb instance: testuser1/test-bulkrna

Ingest data#

Access #

We start by simulating a nf-core RNA-seq run which yields us a count matrix artifact.

(See Nextflow for running this with Nextflow.)

# pretend we're running a bulk RNA-seq pipeline
ln.track(ln.Transform(name="nf-core RNA-seq", reference="https://nf-co.re/rnaseq"))
# create a directory for its output
Path("./test-bulkrna/output_dir").mkdir(exist_ok=True)
# get the count matrix
path = ln.dev.datasets.file_tsv_rnaseq_nfcore_salmon_merged_gene_counts(
    populate_registries=True
)
# move it into the output directory
path = path.rename(f"./test-bulkrna/output_dir/{path.name}")
# register it
ln.Artifact(path, description="Merged Bulk RNA counts").save()
Hide code cell output
💡 saved: Transform(uid='ldQZrHC30Jco8Ern', name='nf-core RNA-seq', type=notebook, reference='https://nf-co.re/rnaseq', updated_at=2024-03-04 17:09:20 UTC, created_by_id=1)
💡 saved: Run(uid='J7W3ylluHXDrdesST1xT', run_at=2024-03-04 17:09:20 UTC, transform_id=1, created_by_id=1)
❗ file has more than one suffix (path.suffixes), using only last suffix: '.tsv' - if you want your file format to be recognized, make an issue: https://github.com/laminlabs/lamindb/issues/new

Transform #

ln.transform.stem_uid = "s5V0dNMVwL9i"
ln.transform.version = "0"
ln.track()
💡 notebook imports: anndata==0.9.2 bionty==0.41.0 lamindb==0.68.0 pandas==1.5.3
💡 saved: Transform(uid='s5V0dNMVwL9i6K79', name='Bulk RNA-seq', short_name='bulkrna', version='0', type=notebook, updated_at=2024-03-04 17:09:25 UTC, created_by_id=1)
💡 saved: Run(uid='w44CFI50r4R7U8jx0Tq1', run_at=2024-03-04 17:09:25 UTC, transform_id=2, created_by_id=1)

Let’s query the artifact:

artifact = ln.Artifact.filter(description="Merged Bulk RNA counts").one()
df = artifact.load()

If we look at it, we realize it deviates far from the tidy data standard Wickham14, conventions of statistics & machine learning Hastie09, Murphy12 and the major Python & R data packages.

Variables are not in columns and observations are not in rows:

df
gene_id gene_name RAP1_IAA_30M_REP1 RAP1_UNINDUCED_REP1 RAP1_UNINDUCED_REP2 WT_REP1 WT_REP2
0 Gfp_transgene_gene Gfp_transgene_gene 0.0 0.000 0.0 0.0 0.0
1 HRA1 HRA1 0.0 8.572 0.0 0.0 0.0
2 snR18 snR18 3.0 8.000 4.0 8.0 8.0
3 tA(UGC)A TGA1 0.0 0.000 0.0 0.0 0.0
4 tL(CAA)A SUP56 0.0 0.000 0.0 0.0 0.0
... ... ... ... ... ... ... ...
120 YAR064W YAR064W 0.0 2.000 0.0 0.0 0.0
121 YAR066W YAR066W 3.0 13.000 8.0 5.0 11.0
122 YAR068W YAR068W 9.0 28.000 24.0 5.0 7.0
123 YAR069C YAR069C 0.0 0.000 0.0 0.0 1.0
124 YAR070C YAR070C 0.0 0.000 0.0 0.0 0.0

125 rows × 7 columns

Let’s change that and move observations into rows:

df = df.T

df
0 1 2 3 4 5 6 7 8 9 ... 115 116 117 118 119 120 121 122 123 124
gene_id Gfp_transgene_gene HRA1 snR18 tA(UGC)A tL(CAA)A tP(UGG)A tS(AGA)A YAL001C YAL002W YAL003W ... YAR050W YAR053W YAR060C YAR061W YAR062W YAR064W YAR066W YAR068W YAR069C YAR070C
gene_name Gfp_transgene_gene HRA1 snR18 TGA1 SUP56 TRN1 tS(AGA)A TFC3 VPS8 EFB1 ... FLO1 YAR053W YAR060C YAR061W YAR062W YAR064W YAR066W YAR068W YAR069C YAR070C
RAP1_IAA_30M_REP1 0.0 0.0 3.0 0.0 0.0 0.0 1.0 55.0 36.0 632.0 ... 4.357 0.0 1.0 0.0 1.0 0.0 3.0 9.0 0.0 0.0
RAP1_UNINDUCED_REP1 0.0 8.572 8.0 0.0 0.0 0.0 0.0 72.0 33.0 810.0 ... 15.72 0.0 0.0 0.0 3.0 2.0 13.0 28.0 0.0 0.0
RAP1_UNINDUCED_REP2 0.0 0.0 4.0 0.0 0.0 0.0 0.0 115.0 82.0 1693.0 ... 13.772 0.0 4.0 0.0 2.0 0.0 8.0 24.0 0.0 0.0
WT_REP1 0.0 0.0 8.0 0.0 0.0 1.0 0.0 60.0 63.0 1115.0 ... 13.465 0.0 0.0 0.0 1.0 0.0 5.0 5.0 0.0 0.0
WT_REP2 0.0 0.0 8.0 0.0 0.0 0.0 0.0 30.0 25.0 704.0 ... 6.891 0.0 1.0 0.0 0.0 0.0 11.0 7.0 1.0 0.0

7 rows × 125 columns

Now, it’s clear that the first two rows are in fact no observations, but descriptions of the variables (or features) themselves.

Let’s create an AnnData object to model this. First, create a dataframe for the variables:

var = pd.DataFrame({"gene_name": df.loc["gene_name"].values}, index=df.loc["gene_id"])
var.head()
gene_name
gene_id
Gfp_transgene_gene Gfp_transgene_gene
HRA1 HRA1
snR18 snR18
tA(UGC)A TGA1
tL(CAA)A SUP56

Now, let’s create an AnnData:

# we're also fixing the datatype here, which was string in the tsv
adata = ad.AnnData(df.iloc[2:].astype("float32"), var=var)

adata
AnnData object with n_obs × n_vars = 5 × 125
    var: 'gene_name'

The AnnData object is in tidy form and complies with conventions of statistics and machine learning:

adata.to_df()
gene_id Gfp_transgene_gene HRA1 snR18 tA(UGC)A tL(CAA)A tP(UGG)A tS(AGA)A YAL001C YAL002W YAL003W ... YAR050W YAR053W YAR060C YAR061W YAR062W YAR064W YAR066W YAR068W YAR069C YAR070C
RAP1_IAA_30M_REP1 0.0 0.000 3.0 0.0 0.0 0.0 1.0 55.0 36.0 632.0 ... 4.357 0.0 1.0 0.0 1.0 0.0 3.0 9.0 0.0 0.0
RAP1_UNINDUCED_REP1 0.0 8.572 8.0 0.0 0.0 0.0 0.0 72.0 33.0 810.0 ... 15.720 0.0 0.0 0.0 3.0 2.0 13.0 28.0 0.0 0.0
RAP1_UNINDUCED_REP2 0.0 0.000 4.0 0.0 0.0 0.0 0.0 115.0 82.0 1693.0 ... 13.772 0.0 4.0 0.0 2.0 0.0 8.0 24.0 0.0 0.0
WT_REP1 0.0 0.000 8.0 0.0 0.0 1.0 0.0 60.0 63.0 1115.0 ... 13.465 0.0 0.0 0.0 1.0 0.0 5.0 5.0 0.0 0.0
WT_REP2 0.0 0.000 8.0 0.0 0.0 0.0 0.0 30.0 25.0 704.0 ... 6.891 0.0 1.0 0.0 0.0 0.0 11.0 7.0 1.0 0.0

5 rows × 125 columns

Validate #

Let’s create a Artifact object from this AnnData.

Almost all gene IDs are validated:

genes = bt.Gene.from_values(
    adata.var.index,
    bt.Gene.stable_id,
    organism="saccharomyces cerevisiae",  # or set globally with bt.settings.organism
)
Hide code cell output
did not create Gene records for 2 non-validated stable_ids: 'Gfp_transgene_gene', 'YAR062W'
# also register the 2 non-validated genes obtained from Bionty
ln.save(genes)

Register #

efs = bt.ExperimentalFactor.lookup()
organism = bt.Organism.lookup()
features = ln.Feature.lookup()
curated_file = ln.Artifact.from_anndata(
    adata,
    description="Curated bulk RNA counts"
)

Hence, let’s save this artifact:

curated_file.save()

Link to validated metadata records:

curated_file.features.add_from_anndata(var_field=bt.Gene.stable_id, organism="saccharomyces cerevisiae")
2 terms (1.60%) are not validated for stable_id: Gfp_transgene_gene, YAR062W
curated_file.labels.add(efs.rna_seq, features.assay)
curated_file.labels.add(organism.saccharomyces_cerevisiae, features.organism)
curated_file.describe()
Artifact(uid='jOUGOQEZkvNarw3k6EhQ', suffix='.h5ad', accessor='AnnData', description='Curated bulk RNA counts', size=28180, hash='6bieh8XjOCCz6bJToN4u1g', hash_type='md5', visibility=1, key_is_virtual=True, updated_at=2024-03-04 17:09:27 UTC)

Provenance:
  🗃️ storage: Storage(uid='pfpl5ptM', root='/home/runner/work/lamin-usecases/lamin-usecases/docs/test-bulkrna', type='local', updated_at=2024-03-04 17:09:17 UTC, created_by_id=1)
  💫 transform: Transform(uid='s5V0dNMVwL9i6K79', name='Bulk RNA-seq', short_name='bulkrna', version='0', type=notebook, updated_at=2024-03-04 17:09:25 UTC, created_by_id=1)
  👣 run: Run(uid='w44CFI50r4R7U8jx0Tq1', run_at=2024-03-04 17:09:25 UTC, transform_id=2, created_by_id=1)
  👤 created_by: User(uid='DzTjkKse', handle='testuser1', name='Test User1', updated_at=2024-03-04 17:09:17 UTC)
Features:
  var: FeatureSet(uid='dyHwHBVm3vMLca8sFVvy', n=123, type='number', registry='bionty.Gene', hash='Se1pVG-akmzgoXOqCxhW', updated_at=2024-03-04 17:09:27 UTC, created_by_id=1)
    'None', 'None', 'TGA1', 'SUP56', 'TRN1', 'None', 'TFC3', 'VPS8', 'EFB1', 'None', 'SSA1', 'ERP2', 'FUN14', 'SPO7', 'MDM10', 'SWC3', 'CYS3', 'DEP1', 'SYN8', 'NTG1', ...
  external: FeatureSet(uid='cjAKfNHXDWeP4r8X3hQF', n=2, registry='core.Feature', hash='9d2OAIOoOUVicxSdf_tN', updated_at=2024-03-04 17:09:27 UTC, created_by_id=1)
    🔗 assay (1, bionty.ExperimentalFactor): 'RNA-Seq'
    🔗 organism (1, bionty.Organism): 'saccharomyces cerevisiae'
Labels:
  🏷️ organism (1, bionty.Organism): 'saccharomyces cerevisiae'
  🏷️ experimental_factors (1, bionty.ExperimentalFactor): 'RNA-Seq'

Query data#

We have two files in the artifact registry:

ln.Artifact.df()
uid storage_id key suffix accessor description version size hash hash_type n_objects n_observations transform_id run_id visibility key_is_virtual created_at updated_at created_by_id
id
2 jOUGOQEZkvNarw3k6EhQ 1 None .h5ad AnnData Curated bulk RNA counts None 28180 6bieh8XjOCCz6bJToN4u1g md5 None None 2 2 1 True 2024-03-04 17:09:27.103607+00:00 2024-03-04 17:09:27.136903+00:00 1
1 coFAm7s8limgb0LxPKw2 1 output_dir/salmon.merged.gene_counts.tsv .tsv None Merged Bulk RNA counts None 3787 xxw0k3au3KtxFcgtbEr4eQ md5 None None 1 1 1 False 2024-03-04 17:09:24.712136+00:00 2024-03-04 17:09:24.712160+00:00 1
curated_file.view_lineage()
_images/7375681436beb326011c6d22d71b61594f8306c53fe32a425b02574d190af371.svg

Let’s by query by gene:

genes = bt.Gene.lookup()
genes.spo7
Gene(uid='2pkcLeMEB6aS', symbol='SPO7', stable_id='YAL009W', ncbi_gene_ids='851224', biotype='protein_coding', description='Putative regulatory subunit of Nem1p-Spo7p phosphatase holoenzyme; regulates nuclear growth by controlling phospholipid biosynthesis, required for normal nuclear envelope morphology, premeiotic replication, and sporulation ', synonyms='', updated_at=2024-03-04 17:09:27 UTC, organism_id=1, public_source_id=13, created_by_id=1)
# a gene set containing SPO7
feature_set = ln.FeatureSet.filter(genes=genes.spo7).first()
# artifacts that link to this feature set
ln.Artifact.filter(feature_sets=feature_set).df()
uid storage_id key suffix accessor description version size hash hash_type n_objects n_observations transform_id run_id visibility key_is_virtual created_at updated_at created_by_id
id
2 jOUGOQEZkvNarw3k6EhQ 1 None .h5ad AnnData Curated bulk RNA counts None 28180 6bieh8XjOCCz6bJToN4u1g md5 None None 2 2 1 True 2024-03-04 17:09:27.103607+00:00 2024-03-04 17:09:27.136903+00:00 1
# clean up test instance
!lamin delete --force test-bulkrna
!rm -r test-bulkrna
Hide code cell output
💡 deleting instance testuser1/test-bulkrna
✅     deleted instance settings file: /home/runner/.lamin/instance--testuser1--test-bulkrna.env
✅     instance cache deleted
✅     deleted '.lndb' sqlite file
❗     consider manually deleting your stored data: /home/runner/work/lamin-usecases/lamin-usecases/docs/test-bulkrna