3.5. Vignette for extracting Long Ids#
Reference Guided Isoform Discovery (Bulk)
This vignettes can be used to extract reference gen and transcript Ids from Gffcompare results and map them back to transcript counts and reconstructed gtf for further analysis. It takes the Gffcompare tracking file, trasncript counts matrix and reconstructed strngtie GTF as an input and outputs the matrix and GTF replacing Stringtie Ids with Long Ids as below.
Long Id format:
Stringtie.Gene.Id^Stringtie.Transcript.Id^Reference.Gene.Id^Reference.Transcript.Id^GffcompareCode
eg:
# Uncomment the commands below to run this cell if using collab
#!git clone https://github.com/MethodsDev/kinnex-documentation-external.git
#!pip install -r kinnex-documentation-external/requirements.txt
# wd="kinnex-documentation-external/test_data/"
import os
import re
import glob
import pandas as pd
import numpy as np
import pyranges as pr
import seaborn as sns
import matplotlib.pyplot as plt
/home/docs/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/__init__.py:8: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', 100)
pd.set_option('display.max_colwidth', 100)
pd.options.display.max_colwidth = 100
Matching gene and transcript ids from Gffcompare outs with Isoquant outs#
!which conda
looking at the tracking file from gffcompare outs for matching ids
#Uncomment if using Binder Jupyter Lab Interface
wd="../test_data/"
colnames=["query_transfrag_id", "query_locus_id", "reference_gene_id","class_code","qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len"]
tracking_df = pd.read_csv(wd+'gffcmp.tracking', sep='\t',names=colnames, header=None)
print(tracking_df.shape)
tracking_df.head()
(25855, 5)
query_transfrag_id
query_locus_id
reference_gene_id
class_code
qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len
0
TCONS_00000001
XLOC_000001
ENSG00000228794.11|ENST00000670780.1
j
q1:MSTRG.5|MSTRG.5.1|7|0.000000|0.000000|0.000000|3219
1
TCONS_00000002
XLOC_000001
ENSG00000228794.11|ENST00000445118.7
j
q1:MSTRG.5|MSTRG.5.2|6|0.000000|0.000000|0.000000|3060
2
TCONS_00000003
XLOC_000001
ENSG00000228794.11|ENST00000445118.7
j
q1:MSTRG.5|MSTRG.5.3|6|0.000000|0.000000|0.000000|3060
3
TCONS_00000004
XLOC_000002
ENSG00000187608.10|ENST00000649529.1
=
q1:MSTRG.6|MSTRG.6.1|2|0.000000|0.000000|0.000000|700
4
TCONS_00000005
XLOC_000003
ENSG00000217801.11|ENST00000456409.6
=
q1:MSTRG.8|MSTRG.8.1|4|0.000000|0.000000|0.000000|937
#col0 gene_ids qJ:gene_id
tracking_df["gene_id"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[0].split(':')[1])
#col1 transcript_id
tracking_df["transcript_id"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[1])
#col2 num_exons
tracking_df["num_exons"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[2])
#col3 FPKM
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[3]).unique()
#col4 TPM
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[4]).unique()
#col5 cov
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[5]).unique()
#col6 len
tracking_df["len"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[6])
tracking_df=tracking_df[['query_transfrag_id','query_locus_id','reference_gene_id','class_code','gene_id','transcript_id','num_exons','len']]
tracking_df.transcript_id.unique().shape
(25855,)
tracking_df
#check - gffcompare matching at gene vs transcript level
query_transfrag_id
query_locus_id
reference_gene_id
class_code
gene_id
transcript_id
num_exons
len
0
TCONS_00000001
XLOC_000001
ENSG00000228794.11|ENST00000670780.1
j
MSTRG.5
MSTRG.5.1
7
3219
1
TCONS_00000002
XLOC_000001
ENSG00000228794.11|ENST00000445118.7
j
MSTRG.5
MSTRG.5.2
6
3060
2
TCONS_00000003
XLOC_000001
ENSG00000228794.11|ENST00000445118.7
j
MSTRG.5
MSTRG.5.3
6
3060
3
TCONS_00000004
XLOC_000002
ENSG00000187608.10|ENST00000649529.1
=
MSTRG.6
MSTRG.6.1
2
700
4
TCONS_00000005
XLOC_000003
ENSG00000217801.11|ENST00000456409.6
=
MSTRG.8
MSTRG.8.1
4
937
...
...
...
...
...
...
...
...
...
25850
TCONS_00025851
XLOC_011477
ENSG00000131002.13|ENST00000693214.1
=
MSTRG.12167
MSTRG.12167.2
3
2533
25851
TCONS_00025852
XLOC_011477
ENSG00000131002.13|ENST00000253320.8
c
MSTRG.12166
MSTRG.12166.1
2
1820
25852
TCONS_00025853
XLOC_011478
ENSG00000198692.10|ENST00000361365.7
=
MSTRG.12169
MSTRG.12169.1
7
1392
25853
TCONS_00025854
XLOC_011479
ENSG00000012817.16|ENST00000541639.5
j
MSTRG.12168
MSTRG.12168.1
8
1527
25854
TCONS_00025855
XLOC_011480
ENSG00000187191.16|ENST00000382365.7
j
MSTRG.12170
MSTRG.12170.1
12
2761
25855 rows × 8 columns
Looking at the distribution of class codes (similar to SQANTI)#
Code Relationship to reference transcript = complete, exact intron chain match c contained in reference transcript (intron compatible) k contains reference transcript (reverse containment) m retained intron(s) compared to reference, full intron chain match
everywhere else n completely overlaps intron from reference transcript, partial or no intron
chain match everywhere else j multi-exon with at least one junction match e single exon that partially covers an intron from reference o other same strand overlap with reference exons s intron match on the opposite strand (likely a mapping error) x exonic overlap on the opposite strand i fully contained within a reference intron y contains a reference within its intron(s) p possible polymerase run-on (close to reference but no overlap) r repeat (at least 50% bases are soft-masked) u none of the above (unknown, intergenic)
sns.histplot(tracking_df.class_code).set(title='Distribution of class codes - Gffcompare \n Ref Guided track Bulk Kinnex ')
[Text(0.5, 1.0, 'Distribution of class codes - Gffcompare \n Ref Guided track Bulk Kinnex ')]
We have combined count matrix, generated with Isoquant in Quantification mode, quantifying against a reconstrutced gtf - stringtie_merged.gtf
isoquant_counts_df = pd.read_csv(wd+'/refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv', sep='\t')
isoquant_counts_df.head(5)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[13], line 1
----> 1 isoquant_counts_df = pd.read_csv(wd+'/refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv', sep='\t')
2 isoquant_counts_df.head(5)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1026, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)
1013 kwds_defaults = _refine_defaults_read(
1014 dialect,
1015 delimiter,
(...) 1022 dtype_backend=dtype_backend,
1023 )
1024 kwds.update(kwds_defaults)
-> 1026 return _read(filepath_or_buffer, kwds)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:620, in _read(filepath_or_buffer, kwds)
617 _validate_names(kwds.get("names", None))
619 # Create the parser.
--> 620 parser = TextFileReader(filepath_or_buffer, **kwds)
622 if chunksize or iterator:
623 return parser
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1620, in TextFileReader.__init__(self, f, engine, **kwds)
1617 self.options["has_index_names"] = kwds["has_index_names"]
1619 self.handles: IOHandles | None = None
-> 1620 self._engine = self._make_engine(f, self.engine)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1880, in TextFileReader._make_engine(self, f, engine)
1878 if "b" not in mode:
1879 mode += "b"
-> 1880 self.handles = get_handle(
1881 f,
1882 mode,
1883 encoding=self.options.get("encoding", None),
1884 compression=self.options.get("compression", None),
1885 memory_map=self.options.get("memory_map", False),
1886 is_text=is_text,
1887 errors=self.options.get("encoding_errors", "strict"),
1888 storage_options=self.options.get("storage_options", None),
1889 )
1890 assert self.handles is not None
1891 f = self.handles.handle
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/common.py:873, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
868 elif isinstance(handle, str):
869 # Check whether the filename is to be opened in binary mode.
870 # Binary mode does not support 'encoding' and 'newline'.
871 if ioargs.encoding and "b" not in ioargs.mode:
872 # Encoding
--> 873 handle = open(
874 handle,
875 ioargs.mode,
876 encoding=ioargs.encoding,
877 errors=errors,
878 newline="",
879 )
880 else:
881 # Binary mode
882 handle = open(handle, ioargs.mode)
FileNotFoundError: [Errno 2] No such file or directory: '../test_data//refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv'
isoquant_counts_df.shape
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 isoquant_counts_df.shape
NameError: name 'isoquant_counts_df' is not defined
isoquant_counts_df.rename(columns={'#feature_id': 'isoform_id'}, inplace=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 isoquant_counts_df.rename(columns={'#feature_id': 'isoform_id'}, inplace=True)
NameError: name 'isoquant_counts_df' is not defined
tracking_df[['reference_gene_id','class_code','transcript_id']]
reference_gene_id
class_code
transcript_id
0
ENSG00000228794.11|ENST00000670780.1
j
MSTRG.5.1
1
ENSG00000228794.11|ENST00000445118.7
j
MSTRG.5.2
2
ENSG00000228794.11|ENST00000445118.7
j
MSTRG.5.3
3
ENSG00000187608.10|ENST00000649529.1
=
MSTRG.6.1
4
ENSG00000217801.11|ENST00000456409.6
=
MSTRG.8.1
...
...
...
...
25850
ENSG00000131002.13|ENST00000693214.1
=
MSTRG.12167.2
25851
ENSG00000131002.13|ENST00000253320.8
c
MSTRG.12166.1
25852
ENSG00000198692.10|ENST00000361365.7
=
MSTRG.12169.1
25853
ENSG00000012817.16|ENST00000541639.5
j
MSTRG.12168.1
25854
ENSG00000187191.16|ENST00000382365.7
j
MSTRG.12170.1
25855 rows × 3 columns
print(isoquant_counts_df.shape)
# 59734
isoquant_counts_df[isoquant_counts_df.isoform_id.str.contains('_')]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 print(isoquant_counts_df.shape)
2 # 59734
3 isoquant_counts_df[isoquant_counts_df.isoform_id.str.contains('_')]
NameError: name 'isoquant_counts_df' is not defined
isoquant_wRefids=isoquant_counts_df[~isoquant_counts_df.isoform_id.str.contains('_')].merge(tracking_df[['reference_gene_id','class_code','transcript_id']], \
how='left', left_on='isoform_id', right_on='transcript_id')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 isoquant_wRefids=isoquant_counts_df[~isoquant_counts_df.isoform_id.str.contains('_')].merge(tracking_df[['reference_gene_id','class_code','transcript_id']], \
2 how='left', left_on='isoform_id', right_on='transcript_id')
NameError: name 'isoquant_counts_df' is not defined
isoquant_wRefids.loc[isoquant_wRefids.isoform_id.isin(['MSTRG.5415.3','MSTRG.12284.4',
'MSTRG.6514.2','MSTRG.6514.1',
'MSTRG.18675.1','MSTRG.19411.3'])]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.isoform_id.isin(['MSTRG.5415.3','MSTRG.12284.4',
2 'MSTRG.6514.2','MSTRG.6514.1',
3 'MSTRG.18675.1','MSTRG.19411.3'])]
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.shape)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 print(isoquant_wRefids.shape)
NameError: name 'isoquant_wRefids' is not defined
#print(isoquant_wRefids.reference_gene_id.head(10))
#print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[-1]).head(10))
#print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[0]).head(10))
isoquant_wRefids[['ref_gene_id','ref_transcript_id']]=\
isoquant_wRefids['reference_gene_id'].str.split('|',expand=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 6
1 #print(isoquant_wRefids.reference_gene_id.head(10))
2 #print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[-1]).head(10))
3 #print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[0]).head(10))
5 isoquant_wRefids[['ref_gene_id','ref_transcript_id']]=\
----> 6 isoquant_wRefids['reference_gene_id'].str.split('|',expand=True)
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',].shape)
isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 print(isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',].shape)
2 isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',]
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids=isoquant_wRefids[['transcript_id','ref_gene_id','ref_transcript_id','class_code',\
'HG002_E0','HG002_E0','HG002_E2','K562_E0','K562_E1','K562_E2']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 isoquant_wRefids=isoquant_wRefids[['transcript_id','ref_gene_id','ref_transcript_id','class_code',\
2 'HG002_E0','HG002_E0','HG002_E2','K562_E0','K562_E1','K562_E2']]
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.to_csv(wd+'isoquant_counts_mapped.tsv', sep='\t')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 isoquant_wRefids.to_csv(wd+'isoquant_counts_mapped.tsv', sep='\t')
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.shape)
isoquant_wRefids.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 1
----> 1 print(isoquant_wRefids.shape)
2 isoquant_wRefids.head()
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['strgt_gene_ids']=\
isoquant_wRefids.transcript_id.apply(lambda x: ('.').join(x.split('.')[0:2]))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 2
1 isoquant_wRefids['strgt_gene_ids']=\
----> 2 isoquant_wRefids.transcript_id.apply(lambda x: ('.').join(x.split('.')[0:2]))
NameError: name 'isoquant_wRefids' is not defined
(isoquant_wRefids[['ref_gene_id','strgt_gene_ids']]\
.groupby('strgt_gene_ids')['ref_gene_id'].nunique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 1
----> 1 (isoquant_wRefids[['ref_gene_id','strgt_gene_ids']]\
2 .groupby('strgt_gene_ids')['ref_gene_id'].nunique())
NameError: name 'isoquant_wRefids' is not defined
Gffcompare may assign distinct Reference Genes and Transcripts to the same Stringtie Gene Ids. We would want to preserve that assignemnt.
Here in this sub-section we explore these cases.
stringtie_gene_id_counts=isoquant_wRefids.groupby('strgt_gene_ids').agg({'ref_gene_id':lambda x: len(pd.unique(x))})
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 1
----> 1 stringtie_gene_id_counts=isoquant_wRefids.groupby('strgt_gene_ids').agg({'ref_gene_id':lambda x: len(pd.unique(x))})
NameError: name 'isoquant_wRefids' is not defined
stringtie_gene_id_counts=stringtie_gene_id_counts.add_suffix('_Count').reset_index()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 1
----> 1 stringtie_gene_id_counts=stringtie_gene_id_counts.add_suffix('_Count').reset_index()
NameError: name 'stringtie_gene_id_counts' is not defined
stringtie_gene_id_counts
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 1
----> 1 stringtie_gene_id_counts
NameError: name 'stringtie_gene_id_counts' is not defined
len(isoquant_wRefids.loc\
[isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids)]\
.strgt_gene_ids.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[31], line 1
----> 1 len(isoquant_wRefids.loc\
2 [isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids)]\
3 .strgt_gene_ids.unique())
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc\
[isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids),\
['strgt_gene_ids','transcript_id','ref_gene_id','ref_transcript_id','class_code']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[32], line 1
----> 1 isoquant_wRefids.loc\
2 [isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids),\
3 ['strgt_gene_ids','transcript_id','ref_gene_id','ref_transcript_id','class_code']]
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[33], line 1
----> 1 len(isoquant_wRefids.transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.ref_transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 1
----> 1 len(isoquant_wRefids.ref_transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].ref_transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[35], line 1
----> 1 len(isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].ref_transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].sort_values(by=['ref_transcript_id'])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[36], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].sort_values(by=['ref_transcript_id'])
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[37], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].groupby(['class_code']).size()
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids=isoquant_wRefids.fillna(value="")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[38], line 1
----> 1 isoquant_wRefids=isoquant_wRefids.fillna(value="")
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['isoform_id']=isoquant_wRefids[['transcript_id', 'ref_gene_id', 'ref_transcript_id','class_code']]\
.agg('^'.join, axis=1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[39], line 1
----> 1 isoquant_wRefids['isoform_id']=isoquant_wRefids[['transcript_id', 'ref_gene_id', 'ref_transcript_id','class_code']]\
2 .agg('^'.join, axis=1)
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.head(5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[40], line 1
----> 1 isoquant_wRefids.head(5)
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['isoform_id']=isoquant_wRefids.isoform_id.str.replace("=", "eq")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[41], line 1
----> 1 isoquant_wRefids['isoform_id']=isoquant_wRefids.isoform_id.str.replace("=", "eq")
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids['class_code']=="=",].head(3)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[42], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids['class_code']=="=",].head(3)
NameError: name 'isoquant_wRefids' is not defined
isoquant_counts_combined_wRefids=isoquant_wRefids[['isoform_id',\
'HG002_E0','HG002_E1','HG002_E2','K562_E0','K562_E1','K562_E2']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[43], line 1
----> 1 isoquant_counts_combined_wRefids=isoquant_wRefids[['isoform_id',\
2 'HG002_E0','HG002_E1','HG002_E2','K562_E0','K562_E1','K562_E2']]
NameError: name 'isoquant_wRefids' is not defined
isoquant_counts_combined_wRefids.head(5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[44], line 1
----> 1 isoquant_counts_combined_wRefids.head(5)
NameError: name 'isoquant_counts_combined_wRefids' is not defined
isoquant_counts_combined_wRefids.to_csv(\
wd+'/isoquant_counts_combined_wRefids.tsv', \
sep='\t')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[45], line 1
----> 1 isoquant_counts_combined_wRefids.to_csv(\
2 wd+'/isoquant_counts_combined_wRefids.tsv', \
3 sep='\t')
NameError: name 'isoquant_counts_combined_wRefids' is not defined
Updating the gtf with Ref Ensemble ids
gtf_f=wd+'stringtie_merged.gtf'
# as PyRanges-object
gr_obj = pr.read_gtf(gtf_f)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:230, in skiprows(f)
228 import gzip
--> 230 fh = gzip.open(f)
231 for i, l in enumerate(fh):
File ~/.asdf/installs/python/3.12.10/lib/python3.12/gzip.py:65, in open(filename, mode, compresslevel, encoding, errors, newline)
64 if isinstance(filename, (str, bytes, os.PathLike)):
---> 65 binary_file = GzipFile(filename, gz_mode, compresslevel)
66 elif hasattr(filename, "read") or hasattr(filename, "write"):
File ~/.asdf/installs/python/3.12.10/lib/python3.12/gzip.py:201, in GzipFile.__init__(self, filename, mode, compresslevel, fileobj, mtime)
200 if fileobj is None:
--> 201 fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
202 if filename is None:
FileNotFoundError: [Errno 2] No such file or directory: '../test_data/stringtie_merged.gtf'
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[46], line 3
1 gtf_f=wd+'stringtie_merged.gtf'
2 # as PyRanges-object
----> 3 gr_obj = pr.read_gtf(gtf_f)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:318, in read_gtf(f, full, as_df, nrows, duplicate_attr, ignore_bad)
245 def read_gtf(
246 f,
247 full=True,
(...) 251 ignore_bad: bool = False,
252 ):
253 """Read files in the Gene Transfer Format.
254
255 Parameters
(...) 315 >>> # 18 hidden columns: gene_name, gene_source, gene_biotype, transcript_id, transcript_version, transcript_name, transcript_source, transcript_biotype, tag, transcript_support_level, ... (+ 8 more.)
316 """
--> 318 _skiprows = skiprows(f)
320 if full:
321 gr = read_gtf_full(f, as_df, nrows, _skiprows, duplicate_attr, ignore_bad=ignore_bad)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:235, in skiprows(f)
233 break
234 except (OSError, TypeError): # not a gzipped file, or StringIO
--> 235 fh = open(f)
236 for i, l in enumerate(fh):
237 if l[0] != "#":
FileNotFoundError: [Errno 2] No such file or directory: '../test_data/stringtie_merged.gtf'
gr_df=gr_obj.df
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[47], line 1
----> 1 gr_df=gr_obj.df
NameError: name 'gr_obj' is not defined
print(gr_df.shape)
gr_df.head()
# 244000 number of unique transcripts from gffcompare tracking file
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[48], line 1
----> 1 print(gr_df.shape)
2 gr_df.head()
3 # 244000 number of unique transcripts from gffcompare tracking file
NameError: name 'gr_df' is not defined
gr_df[gr_df.Feature=='transcript'].transcript_id.unique().shape
# 25855 transcripts in gtf
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[49], line 1
----> 1 gr_df[gr_df.Feature=='transcript'].transcript_id.unique().shape
2 # 25855 transcripts in gtf
NameError: name 'gr_df' is not defined
#MSTRG.5
gr_df.loc[gr_df['gene_id']=='MSTRG.5',].transcript_id.unique()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[50], line 2
1 #MSTRG.5
----> 2 gr_df.loc[gr_df['gene_id']=='MSTRG.5',].transcript_id.unique()
NameError: name 'gr_df' is not defined
print(isoquant_wRefids.shape)
isoquant_wRefids.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[51], line 1
----> 1 print(isoquant_wRefids.shape)
2 isoquant_wRefids.head()
NameError: name 'isoquant_wRefids' is not defined
gr_wRefids=gr_df.merge(\
isoquant_wRefids[['transcript_id','strgt_gene_ids','ref_gene_id','ref_transcript_id','isoform_id','class_code']],\
how='left', on='transcript_id')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[52], line 1
----> 1 gr_wRefids=gr_df.merge(\
2 isoquant_wRefids[['transcript_id','strgt_gene_ids','ref_gene_id','ref_transcript_id','isoform_id','class_code']],\
3 how='left', on='transcript_id')
NameError: name 'gr_df' is not defined
print(gr_wRefids.shape)
gr_wRefids.head(5)
# 7618 transcript ids in the gtf do not exist in isoquant results.
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[53], line 1
----> 1 print(gr_wRefids.shape)
2 gr_wRefids.head(5)
3 # 7618 transcript ids in the gtf do not exist in isoquant results.
NameError: name 'gr_wRefids' is not defined
gr_wRefids[['ref_transcript_id','class_code']].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[54], line 1
----> 1 gr_wRefids[['ref_transcript_id','class_code']].groupby(['class_code']).size()
NameError: name 'gr_wRefids' is not defined
gr_wRefids[gr_wRefids['ref_transcript_id']==""].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[55], line 1
----> 1 gr_wRefids[gr_wRefids['ref_transcript_id']==""].groupby(['class_code']).size()
NameError: name 'gr_wRefids' is not defined
gr_wRefids=gr_wRefids.rename(columns={
'strgt_gene_ids':'mstrg_gene',
'transcript_id':'mstrg_transcript',
'isoform_id':'transcript_id',
'ref_gene_id': 'ensmb_gene',
'ref_transcript_id':'ensmb_transcript'})
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[56], line 1
----> 1 gr_wRefids=gr_wRefids.rename(columns={
2 'strgt_gene_ids':'mstrg_gene',
3 'transcript_id':'mstrg_transcript',
4 'isoform_id':'transcript_id',
5 'ref_gene_id': 'ensmb_gene',
6 'ref_transcript_id':'ensmb_transcript'})
NameError: name 'gr_wRefids' is not defined
len(gr_wRefids[gr_wRefids['transcript_id'].isnull()])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[57], line 1
----> 1 len(gr_wRefids[gr_wRefids['transcript_id'].isnull()])
NameError: name 'gr_wRefids' is not defined
gr_wRefids['transcript_id']=gr_wRefids['transcript_id'].fillna("")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[58], line 1
----> 1 gr_wRefids['transcript_id']=gr_wRefids['transcript_id'].fillna("")
NameError: name 'gr_wRefids' is not defined
len(gr_wRefids.loc[(gr_wRefids['class_code']=='u') & (gr_wRefids['Feature']=='transcript'),])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[59], line 1
----> 1 len(gr_wRefids.loc[(gr_wRefids['class_code']=='u') & (gr_wRefids['Feature']=='transcript'),])
NameError: name 'gr_wRefids' is not defined
gr_wRefids=gr_wRefids.drop(['mstrg_gene','mstrg_transcript','ensmb_gene',\
'ensmb_transcript','class_code'], axis=1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[60], line 1
----> 1 gr_wRefids=gr_wRefids.drop(['mstrg_gene','mstrg_transcript','ensmb_gene',\
2 'ensmb_transcript','class_code'], axis=1)
NameError: name 'gr_wRefids' is not defined
gr_wRefids.head(10)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[61], line 1
----> 1 gr_wRefids.head(10)
NameError: name 'gr_wRefids' is not defined
Writing gtf - replacing transcript Ids in strtingtie GTF with Long Ids#
gr_wRefids_pr= pr.PyRanges(gr_wRefids)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[62], line 1
----> 1 gr_wRefids_pr= pr.PyRanges(gr_wRefids)
NameError: name 'gr_wRefids' is not defined
gr_wRefids_pr.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[63], line 1
----> 1 gr_wRefids_pr.head()
NameError: name 'gr_wRefids_pr' is not defined
gr_wRefids_pr.to_gtf(\
wd+"stringtie_merged.wRefIds.gtf")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[64], line 1
----> 1 gr_wRefids_pr.to_gtf(\
2 wd+"stringtie_merged.wRefIds.gtf")
NameError: name 'gr_wRefids_pr' is not defined
#! cat stringtie_merged.wRefIds.gtf | head -5
#Uncomment if using Binder Jupyter Lab Interface
wd="../test_data/"
colnames=["query_transfrag_id", "query_locus_id", "reference_gene_id","class_code","qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len"]
tracking_df = pd.read_csv(wd+'gffcmp.tracking', sep='\t',names=colnames, header=None)
print(tracking_df.shape)
tracking_df.head()
(25855, 5)
| query_transfrag_id | query_locus_id | reference_gene_id | class_code | qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len | |
|---|---|---|---|---|---|
| 0 | TCONS_00000001 | XLOC_000001 | ENSG00000228794.11|ENST00000670780.1 | j | q1:MSTRG.5|MSTRG.5.1|7|0.000000|0.000000|0.000000|3219 |
| 1 | TCONS_00000002 | XLOC_000001 | ENSG00000228794.11|ENST00000445118.7 | j | q1:MSTRG.5|MSTRG.5.2|6|0.000000|0.000000|0.000000|3060 |
| 2 | TCONS_00000003 | XLOC_000001 | ENSG00000228794.11|ENST00000445118.7 | j | q1:MSTRG.5|MSTRG.5.3|6|0.000000|0.000000|0.000000|3060 |
| 3 | TCONS_00000004 | XLOC_000002 | ENSG00000187608.10|ENST00000649529.1 | = | q1:MSTRG.6|MSTRG.6.1|2|0.000000|0.000000|0.000000|700 |
| 4 | TCONS_00000005 | XLOC_000003 | ENSG00000217801.11|ENST00000456409.6 | = | q1:MSTRG.8|MSTRG.8.1|4|0.000000|0.000000|0.000000|937 |
#col0 gene_ids qJ:gene_id
tracking_df["gene_id"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[0].split(':')[1])
#col1 transcript_id
tracking_df["transcript_id"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[1])
#col2 num_exons
tracking_df["num_exons"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[2])
#col3 FPKM
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[3]).unique()
#col4 TPM
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[4]).unique()
#col5 cov
tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[5]).unique()
#col6 len
tracking_df["len"]=tracking_df['qJ:gene_id|transcript_id|num_exons|FPKM|TPM|cov|len'].apply(lambda x: x.split('|')[6])
tracking_df=tracking_df[['query_transfrag_id','query_locus_id','reference_gene_id','class_code','gene_id','transcript_id','num_exons','len']]
tracking_df.transcript_id.unique().shape
(25855,)
tracking_df
#check - gffcompare matching at gene vs transcript level
| query_transfrag_id | query_locus_id | reference_gene_id | class_code | gene_id | transcript_id | num_exons | len | |
|---|---|---|---|---|---|---|---|---|
| 0 | TCONS_00000001 | XLOC_000001 | ENSG00000228794.11|ENST00000670780.1 | j | MSTRG.5 | MSTRG.5.1 | 7 | 3219 |
| 1 | TCONS_00000002 | XLOC_000001 | ENSG00000228794.11|ENST00000445118.7 | j | MSTRG.5 | MSTRG.5.2 | 6 | 3060 |
| 2 | TCONS_00000003 | XLOC_000001 | ENSG00000228794.11|ENST00000445118.7 | j | MSTRG.5 | MSTRG.5.3 | 6 | 3060 |
| 3 | TCONS_00000004 | XLOC_000002 | ENSG00000187608.10|ENST00000649529.1 | = | MSTRG.6 | MSTRG.6.1 | 2 | 700 |
| 4 | TCONS_00000005 | XLOC_000003 | ENSG00000217801.11|ENST00000456409.6 | = | MSTRG.8 | MSTRG.8.1 | 4 | 937 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 25850 | TCONS_00025851 | XLOC_011477 | ENSG00000131002.13|ENST00000693214.1 | = | MSTRG.12167 | MSTRG.12167.2 | 3 | 2533 |
| 25851 | TCONS_00025852 | XLOC_011477 | ENSG00000131002.13|ENST00000253320.8 | c | MSTRG.12166 | MSTRG.12166.1 | 2 | 1820 |
| 25852 | TCONS_00025853 | XLOC_011478 | ENSG00000198692.10|ENST00000361365.7 | = | MSTRG.12169 | MSTRG.12169.1 | 7 | 1392 |
| 25853 | TCONS_00025854 | XLOC_011479 | ENSG00000012817.16|ENST00000541639.5 | j | MSTRG.12168 | MSTRG.12168.1 | 8 | 1527 |
| 25854 | TCONS_00025855 | XLOC_011480 | ENSG00000187191.16|ENST00000382365.7 | j | MSTRG.12170 | MSTRG.12170.1 | 12 | 2761 |
25855 rows × 8 columns
Looking at the distribution of class codes (similar to SQANTI)#
| Code | Relationship to reference transcript |
|---|---|
| = | complete, exact intron chain match |
| c | contained in reference transcript (intron compatible) |
| k | contains reference transcript (reverse containment) |
| m | retained intron(s) compared to reference, full intron chain match
everywhere else |
| n | completely overlaps intron from reference transcript, partial or no intron
chain match everywhere else |
| j | multi-exon with at least one junction match |
| e | single exon that partially covers an intron from reference |
| o | other same strand overlap with reference exons |
| s | intron match on the opposite strand (likely a mapping error) |
| x | exonic overlap on the opposite strand |
| i | fully contained within a reference intron |
| y | contains a reference within its intron(s) |
| p | possible polymerase run-on (close to reference but no overlap) |
| r | repeat (at least 50% bases are soft-masked) |
| u | none of the above (unknown, intergenic) |
sns.histplot(tracking_df.class_code).set(title='Distribution of class codes - Gffcompare \n Ref Guided track Bulk Kinnex ')
[Text(0.5, 1.0, 'Distribution of class codes - Gffcompare \n Ref Guided track Bulk Kinnex ')]
We have combined count matrix, generated with Isoquant in Quantification mode, quantifying against a reconstrutced gtf - stringtie_merged.gtf
isoquant_counts_df = pd.read_csv(wd+'/refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv', sep='\t')
isoquant_counts_df.head(5)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[13], line 1
----> 1 isoquant_counts_df = pd.read_csv(wd+'/refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv', sep='\t')
2 isoquant_counts_df.head(5)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1026, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)
1013 kwds_defaults = _refine_defaults_read(
1014 dialect,
1015 delimiter,
(...) 1022 dtype_backend=dtype_backend,
1023 )
1024 kwds.update(kwds_defaults)
-> 1026 return _read(filepath_or_buffer, kwds)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:620, in _read(filepath_or_buffer, kwds)
617 _validate_names(kwds.get("names", None))
619 # Create the parser.
--> 620 parser = TextFileReader(filepath_or_buffer, **kwds)
622 if chunksize or iterator:
623 return parser
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1620, in TextFileReader.__init__(self, f, engine, **kwds)
1617 self.options["has_index_names"] = kwds["has_index_names"]
1619 self.handles: IOHandles | None = None
-> 1620 self._engine = self._make_engine(f, self.engine)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1880, in TextFileReader._make_engine(self, f, engine)
1878 if "b" not in mode:
1879 mode += "b"
-> 1880 self.handles = get_handle(
1881 f,
1882 mode,
1883 encoding=self.options.get("encoding", None),
1884 compression=self.options.get("compression", None),
1885 memory_map=self.options.get("memory_map", False),
1886 is_text=is_text,
1887 errors=self.options.get("encoding_errors", "strict"),
1888 storage_options=self.options.get("storage_options", None),
1889 )
1890 assert self.handles is not None
1891 f = self.handles.handle
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pandas/io/common.py:873, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
868 elif isinstance(handle, str):
869 # Check whether the filename is to be opened in binary mode.
870 # Binary mode does not support 'encoding' and 'newline'.
871 if ioargs.encoding and "b" not in ioargs.mode:
872 # Encoding
--> 873 handle = open(
874 handle,
875 ioargs.mode,
876 encoding=ioargs.encoding,
877 errors=errors,
878 newline="",
879 )
880 else:
881 # Binary mode
882 handle = open(handle, ioargs.mode)
FileNotFoundError: [Errno 2] No such file or directory: '../test_data//refGuided_isoquants_quantify_combined_trancripts_counts_matrix.tsv'
isoquant_counts_df.shape
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 isoquant_counts_df.shape
NameError: name 'isoquant_counts_df' is not defined
isoquant_counts_df.rename(columns={'#feature_id': 'isoform_id'}, inplace=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 isoquant_counts_df.rename(columns={'#feature_id': 'isoform_id'}, inplace=True)
NameError: name 'isoquant_counts_df' is not defined
tracking_df[['reference_gene_id','class_code','transcript_id']]
| reference_gene_id | class_code | transcript_id | |
|---|---|---|---|
| 0 | ENSG00000228794.11|ENST00000670780.1 | j | MSTRG.5.1 |
| 1 | ENSG00000228794.11|ENST00000445118.7 | j | MSTRG.5.2 |
| 2 | ENSG00000228794.11|ENST00000445118.7 | j | MSTRG.5.3 |
| 3 | ENSG00000187608.10|ENST00000649529.1 | = | MSTRG.6.1 |
| 4 | ENSG00000217801.11|ENST00000456409.6 | = | MSTRG.8.1 |
| ... | ... | ... | ... |
| 25850 | ENSG00000131002.13|ENST00000693214.1 | = | MSTRG.12167.2 |
| 25851 | ENSG00000131002.13|ENST00000253320.8 | c | MSTRG.12166.1 |
| 25852 | ENSG00000198692.10|ENST00000361365.7 | = | MSTRG.12169.1 |
| 25853 | ENSG00000012817.16|ENST00000541639.5 | j | MSTRG.12168.1 |
| 25854 | ENSG00000187191.16|ENST00000382365.7 | j | MSTRG.12170.1 |
25855 rows × 3 columns
print(isoquant_counts_df.shape)
# 59734
isoquant_counts_df[isoquant_counts_df.isoform_id.str.contains('_')]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 print(isoquant_counts_df.shape)
2 # 59734
3 isoquant_counts_df[isoquant_counts_df.isoform_id.str.contains('_')]
NameError: name 'isoquant_counts_df' is not defined
isoquant_wRefids=isoquant_counts_df[~isoquant_counts_df.isoform_id.str.contains('_')].merge(tracking_df[['reference_gene_id','class_code','transcript_id']], \
how='left', left_on='isoform_id', right_on='transcript_id')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 isoquant_wRefids=isoquant_counts_df[~isoquant_counts_df.isoform_id.str.contains('_')].merge(tracking_df[['reference_gene_id','class_code','transcript_id']], \
2 how='left', left_on='isoform_id', right_on='transcript_id')
NameError: name 'isoquant_counts_df' is not defined
isoquant_wRefids.loc[isoquant_wRefids.isoform_id.isin(['MSTRG.5415.3','MSTRG.12284.4',
'MSTRG.6514.2','MSTRG.6514.1',
'MSTRG.18675.1','MSTRG.19411.3'])]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.isoform_id.isin(['MSTRG.5415.3','MSTRG.12284.4',
2 'MSTRG.6514.2','MSTRG.6514.1',
3 'MSTRG.18675.1','MSTRG.19411.3'])]
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.shape)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 print(isoquant_wRefids.shape)
NameError: name 'isoquant_wRefids' is not defined
#print(isoquant_wRefids.reference_gene_id.head(10))
#print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[-1]).head(10))
#print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[0]).head(10))
isoquant_wRefids[['ref_gene_id','ref_transcript_id']]=\
isoquant_wRefids['reference_gene_id'].str.split('|',expand=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 6
1 #print(isoquant_wRefids.reference_gene_id.head(10))
2 #print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[-1]).head(10))
3 #print(isoquant_wRefids.reference_gene_id.apply(lambda x: x.split('|')[0]).head(10))
5 isoquant_wRefids[['ref_gene_id','ref_transcript_id']]=\
----> 6 isoquant_wRefids['reference_gene_id'].str.split('|',expand=True)
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',].shape)
isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 1
----> 1 print(isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',].shape)
2 isoquant_wRefids.loc[isoquant_wRefids['reference_gene_id']!='-',]
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids=isoquant_wRefids[['transcript_id','ref_gene_id','ref_transcript_id','class_code',\
'HG002_E0','HG002_E0','HG002_E2','K562_E0','K562_E1','K562_E2']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 isoquant_wRefids=isoquant_wRefids[['transcript_id','ref_gene_id','ref_transcript_id','class_code',\
2 'HG002_E0','HG002_E0','HG002_E2','K562_E0','K562_E1','K562_E2']]
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.to_csv(wd+'isoquant_counts_mapped.tsv', sep='\t')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 isoquant_wRefids.to_csv(wd+'isoquant_counts_mapped.tsv', sep='\t')
NameError: name 'isoquant_wRefids' is not defined
print(isoquant_wRefids.shape)
isoquant_wRefids.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 1
----> 1 print(isoquant_wRefids.shape)
2 isoquant_wRefids.head()
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['strgt_gene_ids']=\
isoquant_wRefids.transcript_id.apply(lambda x: ('.').join(x.split('.')[0:2]))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 2
1 isoquant_wRefids['strgt_gene_ids']=\
----> 2 isoquant_wRefids.transcript_id.apply(lambda x: ('.').join(x.split('.')[0:2]))
NameError: name 'isoquant_wRefids' is not defined
(isoquant_wRefids[['ref_gene_id','strgt_gene_ids']]\
.groupby('strgt_gene_ids')['ref_gene_id'].nunique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 1
----> 1 (isoquant_wRefids[['ref_gene_id','strgt_gene_ids']]\
2 .groupby('strgt_gene_ids')['ref_gene_id'].nunique())
NameError: name 'isoquant_wRefids' is not defined
Gffcompare may assign distinct Reference Genes and Transcripts to the same Stringtie Gene Ids. We would want to preserve that assignemnt.
Here in this sub-section we explore these cases.
stringtie_gene_id_counts=isoquant_wRefids.groupby('strgt_gene_ids').agg({'ref_gene_id':lambda x: len(pd.unique(x))})
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 1
----> 1 stringtie_gene_id_counts=isoquant_wRefids.groupby('strgt_gene_ids').agg({'ref_gene_id':lambda x: len(pd.unique(x))})
NameError: name 'isoquant_wRefids' is not defined
stringtie_gene_id_counts=stringtie_gene_id_counts.add_suffix('_Count').reset_index()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 1
----> 1 stringtie_gene_id_counts=stringtie_gene_id_counts.add_suffix('_Count').reset_index()
NameError: name 'stringtie_gene_id_counts' is not defined
stringtie_gene_id_counts
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 1
----> 1 stringtie_gene_id_counts
NameError: name 'stringtie_gene_id_counts' is not defined
len(isoquant_wRefids.loc\
[isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids)]\
.strgt_gene_ids.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[31], line 1
----> 1 len(isoquant_wRefids.loc\
2 [isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids)]\
3 .strgt_gene_ids.unique())
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc\
[isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids),\
['strgt_gene_ids','transcript_id','ref_gene_id','ref_transcript_id','class_code']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[32], line 1
----> 1 isoquant_wRefids.loc\
2 [isoquant_wRefids.strgt_gene_ids.isin(stringtie_gene_id_counts[stringtie_gene_id_counts.ref_gene_id_Count>1].strgt_gene_ids),\
3 ['strgt_gene_ids','transcript_id','ref_gene_id','ref_transcript_id','class_code']]
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[33], line 1
----> 1 len(isoquant_wRefids.transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.ref_transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 1
----> 1 len(isoquant_wRefids.ref_transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
len(isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].ref_transcript_id.unique())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[35], line 1
----> 1 len(isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].ref_transcript_id.unique())
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].sort_values(by=['ref_transcript_id'])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[36], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].sort_values(by=['ref_transcript_id'])
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids.\
duplicated(subset=['ref_transcript_id'])].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[37], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids.\
2 duplicated(subset=['ref_transcript_id'])].groupby(['class_code']).size()
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids=isoquant_wRefids.fillna(value="")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[38], line 1
----> 1 isoquant_wRefids=isoquant_wRefids.fillna(value="")
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['isoform_id']=isoquant_wRefids[['transcript_id', 'ref_gene_id', 'ref_transcript_id','class_code']]\
.agg('^'.join, axis=1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[39], line 1
----> 1 isoquant_wRefids['isoform_id']=isoquant_wRefids[['transcript_id', 'ref_gene_id', 'ref_transcript_id','class_code']]\
2 .agg('^'.join, axis=1)
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.head(5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[40], line 1
----> 1 isoquant_wRefids.head(5)
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids['isoform_id']=isoquant_wRefids.isoform_id.str.replace("=", "eq")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[41], line 1
----> 1 isoquant_wRefids['isoform_id']=isoquant_wRefids.isoform_id.str.replace("=", "eq")
NameError: name 'isoquant_wRefids' is not defined
isoquant_wRefids.loc[isoquant_wRefids['class_code']=="=",].head(3)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[42], line 1
----> 1 isoquant_wRefids.loc[isoquant_wRefids['class_code']=="=",].head(3)
NameError: name 'isoquant_wRefids' is not defined
isoquant_counts_combined_wRefids=isoquant_wRefids[['isoform_id',\
'HG002_E0','HG002_E1','HG002_E2','K562_E0','K562_E1','K562_E2']]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[43], line 1
----> 1 isoquant_counts_combined_wRefids=isoquant_wRefids[['isoform_id',\
2 'HG002_E0','HG002_E1','HG002_E2','K562_E0','K562_E1','K562_E2']]
NameError: name 'isoquant_wRefids' is not defined
isoquant_counts_combined_wRefids.head(5)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[44], line 1
----> 1 isoquant_counts_combined_wRefids.head(5)
NameError: name 'isoquant_counts_combined_wRefids' is not defined
isoquant_counts_combined_wRefids.to_csv(\
wd+'/isoquant_counts_combined_wRefids.tsv', \
sep='\t')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[45], line 1
----> 1 isoquant_counts_combined_wRefids.to_csv(\
2 wd+'/isoquant_counts_combined_wRefids.tsv', \
3 sep='\t')
NameError: name 'isoquant_counts_combined_wRefids' is not defined
Updating the gtf with Ref Ensemble ids
gtf_f=wd+'stringtie_merged.gtf'
# as PyRanges-object
gr_obj = pr.read_gtf(gtf_f)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:230, in skiprows(f)
228 import gzip
--> 230 fh = gzip.open(f)
231 for i, l in enumerate(fh):
File ~/.asdf/installs/python/3.12.10/lib/python3.12/gzip.py:65, in open(filename, mode, compresslevel, encoding, errors, newline)
64 if isinstance(filename, (str, bytes, os.PathLike)):
---> 65 binary_file = GzipFile(filename, gz_mode, compresslevel)
66 elif hasattr(filename, "read") or hasattr(filename, "write"):
File ~/.asdf/installs/python/3.12.10/lib/python3.12/gzip.py:201, in GzipFile.__init__(self, filename, mode, compresslevel, fileobj, mtime)
200 if fileobj is None:
--> 201 fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
202 if filename is None:
FileNotFoundError: [Errno 2] No such file or directory: '../test_data/stringtie_merged.gtf'
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[46], line 3
1 gtf_f=wd+'stringtie_merged.gtf'
2 # as PyRanges-object
----> 3 gr_obj = pr.read_gtf(gtf_f)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:318, in read_gtf(f, full, as_df, nrows, duplicate_attr, ignore_bad)
245 def read_gtf(
246 f,
247 full=True,
(...) 251 ignore_bad: bool = False,
252 ):
253 """Read files in the Gene Transfer Format.
254
255 Parameters
(...) 315 >>> # 18 hidden columns: gene_name, gene_source, gene_biotype, transcript_id, transcript_version, transcript_name, transcript_source, transcript_biotype, tag, transcript_support_level, ... (+ 8 more.)
316 """
--> 318 _skiprows = skiprows(f)
320 if full:
321 gr = read_gtf_full(f, as_df, nrows, _skiprows, duplicate_attr, ignore_bad=ignore_bad)
File ~/checkouts/readthedocs.org/user_builds/kinnex-documentation-external/envs/latest/lib/python3.12/site-packages/pyranges/readers.py:235, in skiprows(f)
233 break
234 except (OSError, TypeError): # not a gzipped file, or StringIO
--> 235 fh = open(f)
236 for i, l in enumerate(fh):
237 if l[0] != "#":
FileNotFoundError: [Errno 2] No such file or directory: '../test_data/stringtie_merged.gtf'
gr_df=gr_obj.df
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[47], line 1
----> 1 gr_df=gr_obj.df
NameError: name 'gr_obj' is not defined
print(gr_df.shape)
gr_df.head()
# 244000 number of unique transcripts from gffcompare tracking file
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[48], line 1
----> 1 print(gr_df.shape)
2 gr_df.head()
3 # 244000 number of unique transcripts from gffcompare tracking file
NameError: name 'gr_df' is not defined
gr_df[gr_df.Feature=='transcript'].transcript_id.unique().shape
# 25855 transcripts in gtf
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[49], line 1
----> 1 gr_df[gr_df.Feature=='transcript'].transcript_id.unique().shape
2 # 25855 transcripts in gtf
NameError: name 'gr_df' is not defined
#MSTRG.5
gr_df.loc[gr_df['gene_id']=='MSTRG.5',].transcript_id.unique()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[50], line 2
1 #MSTRG.5
----> 2 gr_df.loc[gr_df['gene_id']=='MSTRG.5',].transcript_id.unique()
NameError: name 'gr_df' is not defined
print(isoquant_wRefids.shape)
isoquant_wRefids.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[51], line 1
----> 1 print(isoquant_wRefids.shape)
2 isoquant_wRefids.head()
NameError: name 'isoquant_wRefids' is not defined
gr_wRefids=gr_df.merge(\
isoquant_wRefids[['transcript_id','strgt_gene_ids','ref_gene_id','ref_transcript_id','isoform_id','class_code']],\
how='left', on='transcript_id')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[52], line 1
----> 1 gr_wRefids=gr_df.merge(\
2 isoquant_wRefids[['transcript_id','strgt_gene_ids','ref_gene_id','ref_transcript_id','isoform_id','class_code']],\
3 how='left', on='transcript_id')
NameError: name 'gr_df' is not defined
print(gr_wRefids.shape)
gr_wRefids.head(5)
# 7618 transcript ids in the gtf do not exist in isoquant results.
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[53], line 1
----> 1 print(gr_wRefids.shape)
2 gr_wRefids.head(5)
3 # 7618 transcript ids in the gtf do not exist in isoquant results.
NameError: name 'gr_wRefids' is not defined
gr_wRefids[['ref_transcript_id','class_code']].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[54], line 1
----> 1 gr_wRefids[['ref_transcript_id','class_code']].groupby(['class_code']).size()
NameError: name 'gr_wRefids' is not defined
gr_wRefids[gr_wRefids['ref_transcript_id']==""].groupby(['class_code']).size()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[55], line 1
----> 1 gr_wRefids[gr_wRefids['ref_transcript_id']==""].groupby(['class_code']).size()
NameError: name 'gr_wRefids' is not defined
gr_wRefids=gr_wRefids.rename(columns={
'strgt_gene_ids':'mstrg_gene',
'transcript_id':'mstrg_transcript',
'isoform_id':'transcript_id',
'ref_gene_id': 'ensmb_gene',
'ref_transcript_id':'ensmb_transcript'})
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[56], line 1
----> 1 gr_wRefids=gr_wRefids.rename(columns={
2 'strgt_gene_ids':'mstrg_gene',
3 'transcript_id':'mstrg_transcript',
4 'isoform_id':'transcript_id',
5 'ref_gene_id': 'ensmb_gene',
6 'ref_transcript_id':'ensmb_transcript'})
NameError: name 'gr_wRefids' is not defined
len(gr_wRefids[gr_wRefids['transcript_id'].isnull()])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[57], line 1
----> 1 len(gr_wRefids[gr_wRefids['transcript_id'].isnull()])
NameError: name 'gr_wRefids' is not defined
gr_wRefids['transcript_id']=gr_wRefids['transcript_id'].fillna("")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[58], line 1
----> 1 gr_wRefids['transcript_id']=gr_wRefids['transcript_id'].fillna("")
NameError: name 'gr_wRefids' is not defined
len(gr_wRefids.loc[(gr_wRefids['class_code']=='u') & (gr_wRefids['Feature']=='transcript'),])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[59], line 1
----> 1 len(gr_wRefids.loc[(gr_wRefids['class_code']=='u') & (gr_wRefids['Feature']=='transcript'),])
NameError: name 'gr_wRefids' is not defined
gr_wRefids=gr_wRefids.drop(['mstrg_gene','mstrg_transcript','ensmb_gene',\
'ensmb_transcript','class_code'], axis=1)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[60], line 1
----> 1 gr_wRefids=gr_wRefids.drop(['mstrg_gene','mstrg_transcript','ensmb_gene',\
2 'ensmb_transcript','class_code'], axis=1)
NameError: name 'gr_wRefids' is not defined
gr_wRefids.head(10)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[61], line 1
----> 1 gr_wRefids.head(10)
NameError: name 'gr_wRefids' is not defined
Writing gtf - replacing transcript Ids in strtingtie GTF with Long Ids#
gr_wRefids_pr= pr.PyRanges(gr_wRefids)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[62], line 1
----> 1 gr_wRefids_pr= pr.PyRanges(gr_wRefids)
NameError: name 'gr_wRefids' is not defined
gr_wRefids_pr.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[63], line 1
----> 1 gr_wRefids_pr.head()
NameError: name 'gr_wRefids_pr' is not defined
gr_wRefids_pr.to_gtf(\
wd+"stringtie_merged.wRefIds.gtf")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[64], line 1
----> 1 gr_wRefids_pr.to_gtf(\
2 wd+"stringtie_merged.wRefIds.gtf")
NameError: name 'gr_wRefids_pr' is not defined
#! cat stringtie_merged.wRefIds.gtf | head -5