IsoTools2 integration with Iso-Seq data

[RoyPollak]

Hi IsoTools team,
Thank you very much for developing and maintaining this valuable tool!

I am trying to integrate PacBio data with IsoTools but am running through some difficulties and be thankful if you could share some of your experience with me.

After getting PacBio output, I implemented Iso-Seq and Pigeon pipelines to create a Seurat object. I filtered out low quality cells based on UMIs and gene expression characteristics. At this point I wish to preform pseudo bulking of the single cells and implement IsoTools in search for differentially spliced genes.

I have previously managed to filter mapped reads from pbmm2 .bam files, and have created an IsoTools object. I have then preformed the relevant steps as you guide in your tutorials to find deferentially spliced genes. Thank you!

Yet, as part of the process I redefined the transcripts with IsoTools algorithm. Currently, I wish to integrate Iso-Seq defined transcripts with IsoTools. Thus, I retrieved from the Seurat object an isoform abundance table (.csv file) and from Pigeon a .gff file. As it is thoroughly explained in your online Transcriptome Import tutorial, I tried executing IsoTools add_sample_from_csv() function. Yet, the function does not fully run.
Here is an exemplary code with the error messages:

isotools_obj = Transcriptome.from_reference("../Mouse/gencode.vM36.annotation.gtf", file_format='auto')
id_map=isotools_obj.add_sample_from_csv(
f”…/aggregated_counts_by_sample.csv", # A file based on a Seurat object
transcripts_file=”Batch001.collapsed.gff”, # A file generated by the Iso-Seq pipeline - https://isoseq.how/classification/isoseq-collapse.html
transcript_id_col='transcript_id',
reconstruct_genes=False,
sep=','
)
TypeError: _read_gff_file() got an unexpected keyword argument 'infer_genes'

I then converted the .gff file into a .gtf file using the gffread Batch001.collapsed.gff -T -o output.collapsed.gtf function and re-run the add_sample_from_csv() using the .gtf file and received this error:
AttributeError: 'Series' object has no attribute 'chr'

I would be very thankful if you could help me figure out how to integrate these tools together. Below you may find my .csv, .gff and .gtf file heads.
Once again, thank you very much for developing and maintaining this valuable tool.

Roy

B.t.w., the link to the files used for the Transcriptome Import tutorial does not work for me.

CSV file:
transcript_id,LPS_Female_sum_coverage,LPS_Male_sum_coverage,Salmonella_Female_sum_coverage,Salmonella_Male_sum_coverage,Steady_State_Female_sum_coverage,Steady_State_Male_sum_coverage
PB.100005.1,5,4,2,10,4,0
PB.100005.100,73,78,243,418,148,252
PB.100005.103,0,1,1,1,1,2
PB.100005.104,0,1,0,2,0,1
PB.100005.107,0,0,0,0,1,2
PB.100005.108,1,0,0,2,0,0
PB.100005.115,0,1,0,0,0,0
PB.100005.116,0,0,0,1,0,0
PB.100005.120,0,0,0,0,0,0

GFF file:
##pacbio-collapse-version 1.0
##date Mon Sep 9 06:22:33 2024 UTC
chr1 PacBio transcript 3905348 3905724 . - . gene_id "PB.1"; transcript_id "PB.1.1";
chr1 PacBio exon 3905348 3905724 . - . gene_id "PB.1"; transcript_id "PB.1.1";
chr1 PacBio transcript 3911074 3911560 . - . gene_id "PB.2"; transcript_id "PB.2.1";
chr1 PacBio exon 3911074 3911560 . - . gene_id "PB.2"; transcript_id "PB.2.1";
chr1 PacBio transcript 3966446 3966898 . - . gene_id "PB.3"; transcript_id "PB.3.1";
chr1 PacBio exon 3966446 3966898 . - . gene_id "PB.3"; transcript_id "PB.3.1";
chr1 PacBio transcript 3971921 3972757 . - . gene_id "PB.4"; transcript_id "PB.4.1";
chr1 PacBio exon 3971921 3972757 . - . gene_id "PB.4"; transcript_id "PB.4.1";

GTF file:
chr1 PacBio transcript 3250396 3250740 . + . transcript_id "PB.7.1"; gene_id "PB.7"
chr1 PacBio exon 3250396 3250740 . + . transcript_id "PB.7.1"; gene_id "PB.7";
chr1 PacBio transcript 3905348 3905724 . - . transcript_id "PB.1.1"; gene_id "PB.1"
chr1 PacBio exon 3905348 3905724 . - . transcript_id "PB.1.1"; gene_id "PB.1";
chr1 PacBio transcript 3911074 3911560 . - . transcript_id "PB.2.1"; gene_id "PB.2"
chr1 PacBio exon 3911074 3911560 . - . transcript_id "PB.2.1"; gene_id "PB.2";
chr1 PacBio transcript 3966446 3966898 . - . transcript_id "PB.3.1"; gene_id "PB.3"
chr1 PacBio exon 3966446 3966898 . - . transcript_id "PB.3.1"; gene_id "PB.3";
chr1 PacBio transcript 3971921 3972757 . - . transcript_id "PB.4.1"; gene_id "PB.4"
chr1 PacBio exon 3971921 3972757 . - . transcript_id "PB.4.1"; gene_id "PB.4";

[mpizzagalli777]

Hi @RoyPollak I have run into the same issue. Where you able to resolve this?

[YalanBi]

Dear @RoyPollak and @mpizzagalli777,

Thank you for reporting the issue.
If you would like to import transcripts identified using external packages in gff/gff3 format, you need to have an index .tbi file for your gff3 file. This can be generated using Tabix or other packages, and should be saved in the same directory as your gff3 file.
If you use a gtf file for importing, it should work. I could not reproduce the error reported by @RoyPollak using the gtf files generated by some external algorithms that I have. If you could provide me with a subset of your CSV and GTF files containing matching transcripts, I will be able to look into it.
We are currently working on implementing the gff/gff3 import.

Regards,
Yalan

[mpizzagalli777]

Thanks so much for the quick reply. Just to clarify the workflow for analyzing a transcriptome generated from an external package (from here https://isotools.readthedocs.io/en/latest/notebooks/03b_transcriptome_import.html):

1. import the transcriptome using from_reference (this is the gencode gff3 with the .tbi in the same folder).
2. add_sample_from_csv
   coverage_csv_file is the counts table, containing a column with TxIDs and then other columns corresponding to each sample
   transcripts_file is the gtf that is output from the external tool and contains the identified transcripts (both novel and known)

When we run this, we are getting a warning for each transcript (both novel and known) that there is missing infos in the gtf:

WARNING:skipping transcript BambuTx97 from gene ENSG00000159592.11, missing infos in gtf: ENSG00000159592.11
WARNING:skipping transcript BambuTx974 from gene ENSG00000132424.17, missing infos in gtf: ENSG00000132424.17
WARNING:skipping transcript ENST00000000233.10 from gene ENSG00000004059.11, missing infos in gtf: ENSG00000004059.11
WARNING:skipping transcript ENST00000000412.8 from gene ENSG00000003056.8, missing infos in gtf: ENSG00000003056.8
WARNING:skipping transcript ENST00000001008.6 from gene ENSG00000004478.8, missing infos in gtf: ENSG00000004478.8

Do you know what might be causing this?

And to provide the full error message above, this is what is output:

---

AttributeError Traceback (most recent call last)
/tmp/ipykernel_2313414/100262762.py in ?()
35 "H9_2": 'H9_2',
36 "H9_3": 'H9_3'
37 }
38
---> 39 id_map = isoseq.add_sample_from_csv(
40 coverage_csv_file = counts_csv,
41 transcripts_file=transcripts_gtf,
42 sample_cov_cols=read_count_cols,

~/.conda/envs/medaka/lib/python3.10/site-packages/isotools/_transcriptome_io.py in ?(self, coverage_csv_file, transcripts_file, transcript_id_col, sample_cov_cols, sample_properties, add_chromosomes, infer_genes, reconstruct_genes, fuzzy_junction, min_exonic_ref_coverage, sep)
329 row.transcript_id,
330 row.gene_id,
331 e.args[0],
332 )
--> 333 except AssertionError as e:
334 logger.warning(
335 "skipping transcript %s from gene %s: duplicate transcript id; Error: %s",
336 row.transcript_id,

~/.local/lib/python3.10/site-packages/pandas/core/generic.py in ?(self, name)
6295 and name not in self._accessors
6296 and self._info_axis._can_hold_identifiers_and_holds_name(name)
6297 ):
6298 return self[name]
-> 6299 return object.getattribute(self, name)

AttributeError: 'Series' object has no attribute 'chr'

[mpizzagalli777]

Hi, Sorry for all of the messages but I thought it might be helpful to mention that adding the infer_gene parameter fixed this error and we were able to import the data.

However, this brought up a different question. How do the infer_gene and reconstruct_gene parameters differ in terms of output and when should each of those be implemented? Sorry for the basic question but wanted to clarify their use cases. We are importing a transcriptome generated using Bambu.

infer_genes – If True, gene structure is inferred from the transcripts. Useful for gtf files without gene information.
reconstruct_genes – If True, transcript gene assignment from gtf is ignored, and transcripts are grouped to genes from scratch.