OSCR

An in vivo parallelized reporter assay to uncover tissue-specific splicing regulatory sequences in a multicellular animal.

Code ↔ Paper

1 match between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 1 match
  1. [1] § Materials and methods › Phylogenetically averaged motif score analysis on C. elegans neuronal/muscle switch-like splicing events ↔ gff_exon-intron_annotations.py, lines 567–615 · score 0.68 · retained introns, skipped exons, ss exons, GFF, metadata, splicing event

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 736 lines · 44 KB · Apache-2.0 · 1 match

  1. import os
  2. import sys
  3. import numpy as np
  4. import matplotlib.pyplot as plt
  6. gff_file = sys.argv[1]
  7. keywords_file = sys.argv[2]
  8. intron_annot_file = sys.argv[3]
  9. exon_intron_annot_file = sys.argv[4]
  10. intron_cords_file = sys.argv[5]
  11. coordinate_type = sys.argv[6] # relative coordinates or absolute coordinates.
  12. median = int(sys.argv[7])
  14. # Depending on the length you want your sequences to be resized to, you can either choose to enter a length here
  15. # alternatively if you want to use the median length for resizing you can comment this out and use the code below on
  16. # median = round_to_nearest_multiple_of_10(median_size)
  17. # To keep it simple, I have used the variable name median and it depends on user to use the actual median size or
  18. # specify a length they'd like the introns to be resized to.
  20. # CAN BE CONVERTED TO BED FILE EASILY WITH EXCEL AND R - crdnts_to_bed.r ###
  23. #Simple logic here is if the exon spans an intron, it likely is a retained intron event whereas if an intron spans an exon, that's not an actual intron and gets filtered out by intron_filter().
  24. def retained_introns(model):
  25. model['splicing_event'] = np.nan
  26. exons = model[model['group'] == 'exon']
  27. introns = model[model['group'] == 'intron']
  28. exoncoords = set(map(tuple, exons[['start', 'end']].values))
  30. for index, row in introns.iterrows():
  31. for exon_start, exon_end in exoncoords:
  32. if exon_start < row['start'] and exon_end > row['end']:
  33. model.loc[index, 'splicing_event'] = 'retained intron'
  34. exon_to_delete = exons[(exons['start'] == exon_start) & (exons['end'] == exon_end)]
  35. model.loc[exon_to_delete.index, 'splicing_event'] = 'drop'
  37. # I am dropping those exons which aren't real exons but two exons flanking a retained intron.
  38. model = model[model['splicing_event'] != 'drop']
  39. return model
  42. def exon_cat(gene_model):
  43. group=gene_model[gene_model['group']=='exon']
  44. group.sort_values(by=['transcript_id'],inplace=True)
  45. group['transcript_id'] = group['transcript_id'].astype(str).str.replace('transcript_id=','')
  46. smallest_exon_start_tr = group.groupby(['transcript_id'])[['start','end']].min().reset_index()
  47. greatest_exon_end_tr = group.groupby(['transcript_id'])[['start','end']].max().reset_index()
  49. group['splicing_event'] = np.nan
  50. smallest_exon_rows = group[group[['transcript_id','start', 'end']].apply(tuple, axis=1).isin(smallest_exon_start_tr[['transcript_id','start', 'end']].apply(tuple, axis=1))]
  51. greatest_exon_rows = group[group[['transcript_id','start', 'end']].apply(tuple, axis=1).isin(greatest_exon_end_tr[['transcript_id','start', 'end']].apply(tuple, axis=1))]
  53. if '+' in group['strand'].unique():
  54. group.loc[smallest_exon_rows.index, 'splicing_event'] = 'first_exon'
  55. group.loc[greatest_exon_rows.index, 'splicing_event'] = 'last_exon'
  57. elif '-' in group['strand'].unique():
  58. group.loc[smallest_exon_rows.index, 'splicing_event'] = 'last_exon'
  59. group.loc[greatest_exon_rows.index, 'splicing_event'] = 'first_exon'
  61. excluded = group[
  62. (group['splicing_event'] == 'first_exon') |
  63. (group['splicing_event'] == 'last_exon')
  64. ][['start','end', 'transcript_id']]
  66. excluded_set = set(map(tuple, excluded.values))
  68. excluded_tuples = group[['start', 'end', 'transcript_id']].apply(tuple, axis=1).isin(excluded_set)
  69. filtered_group = group[~(excluded_tuples)]
  70. constitutive_alltrs_dic={}
  71. tr_numfilter = filtered_group['transcript_id'].nunique()
  72. alltrexons = filtered_group.groupby(['start','end'])['transcript_id'].count().reset_index(name='transcript_id_count')
  73. Const_overall_trs_condition = alltrexons[alltrexons['transcript_id_count'] == tr_numfilter]
  74. if not Const_overall_trs_condition.empty:
  75. constitutive_alltrs_tuples = set(Const_overall_trs_condition[['start', 'end']].itertuples(index=False, name=None))
  76. else:
  77. constitutive_alltrs_tuples=set()
  79. #Processing skipped and constitutive separately.
  80. first_exon_rows = group[group['splicing_event'] == 'first_exon']
  81. first_exon_tuples = first_exon_rows[['start','end','transcript_id']].apply(tuple, axis=1)
  84. # Grouping transcripts by first exons.
  85. tr_subgroup_df = first_exon_rows.groupby(['start','end'])['transcript_id'].apply(list).reset_index()
  86. tr_subgroup = tr_subgroup_df['transcript_id'].to_dict()
  88. skipped_dic = {}
  89. constitutive_dic = {}
  90. # group transcripts into sub-groups based on first exon and then locate the exon under question in the transcripts that have first exons coordinates smaller than exon in quesstion !
  91. # I am creating transcript subgroups with all exons except the first ones. Hence, I use filtered_group for subgroup.
  92. for key in tr_subgroup:
  93. subgroup = group[group['transcript_id'].isin(tr_subgroup[key])]
  94. tr_num = subgroup['transcript_id'].nunique()
  95. subgroup_filter =subgroup[~(subgroup['splicing_event']=='first_exon') & ~(subgroup['splicing_event']=='last_exon')]
  96. ex_counts = subgroup_filter.groupby(['start', 'end'])['transcript_id'].count().reset_index(name='transcript_id_count')
  97. skipped_condition = ex_counts[ex_counts['transcript_id_count'] < tr_num]
  98. skipped = list(skipped_condition[['start', 'end']].itertuples(index=False, name=None))
  99. skipped_dic[key] = skipped
  100. constitutive_condition = ex_counts[ex_counts['transcript_id_count'] == tr_num]
  101. constitutive = list((constitutive_condition[['start', 'end']].itertuples(index=False, name=None)))
  102. constitutive_dic[key]=constitutive
  104. # Next approach will be parsing out the true constitutive and skipped based on first exon coordinates.
  105. constitutive_true_dic={}
  106. constitutive_tuples_pre = set(tuple_ for tuples_list in constitutive_dic.values() for tuple_ in tuples_list)
  108. for coords in constitutive_tuples_pre:
  109. start, end = coords
  110. if '+' in first_exon_rows['strand'].unique():
  111. # With this line of code I am checking if the exon in question is an internal exon in any of the transcripts.
  112. # Only when I locate transcript block by its first exon and then if the exon in question is internal to it ( start > tup[1] (or first_exon end)) then, I extract the transcript_ids
  113. subset_condition = group.loc[
  114. group.apply(
  115. lambda row: any(start > tup[1]
  116. for tup in first_exon_tuples
  117. if ((row['end'] == tup[1]) & (row['splicing_event'] == 'first_exon'))
  118. ),axis=1)]
  119. trs = subset_condition['transcript_id'].unique()
  120. subset_trs = filtered_group[filtered_group['transcript_id'].isin(trs)]
  122. elif '-' in first_exon_rows['strand'].unique():
  123. subset_condition = group.loc[
  124. group.apply(
  125. lambda row: any(start < tup[1]
  126. for tup in first_exon_tuples
  127. if ((row['end'] == tup[1]) & (row['splicing_event'] == 'first_exon'))
  128. ),axis=1)]
  129. trs = subset_condition['transcript_id'].unique()
  130. subset_trs = filtered_group[filtered_group['transcript_id'].isin(trs)]
  132. second_ex_counts = subset_trs.groupby(['start', 'end'])['transcript_id'].count().reset_index(name='transcript_id_count')
  133. tr_num_subset = subset_trs['transcript_id'].nunique()
  135. skipped_sec_condition = second_ex_counts[second_ex_counts['transcript_id_count'] < tr_num_subset][['start','end']]
  136. if start in list(skipped_sec_condition['start']):
  137. skipped_dic[start]=[]
  138. skipped_dic[start].append(coords)
  140. const_real = second_ex_counts[second_ex_counts['transcript_id_count'] == tr_num_subset][['start','end']]
  141. if start in list(const_real['start']):
  142. constitutive_true_dic[start]=[]
  143. constitutive_true_dic[start].append(coords)
  145. skipped_tuples = set(tuple_ for tuples_list in skipped_dic.values() for tuple_ in tuples_list)
  146. constitutive_tuples = set(tuple_ for tuples_list in constitutive_true_dic.values() for tuple_ in tuples_list)
  148. #I am accounting for first and last exons here, they need to be excluded when I am trying to classify exons of other transcripts in which the first exon of trs1 might be a second or third exon.
  149. tr_num = filtered_group['transcript_id'].nunique()
  150. ex_1start_more_ends = filtered_group.groupby(['start']).filter(lambda g: g['end'].nunique() > 1)
  151. ex_1start_more_ends_count = ex_1start_more_ends.groupby(['start'])['transcript_id'].transform('count')
  152. composite_condition1 = ex_1start_more_ends[(ex_1start_more_ends_count < tr_num)]
  154. ex_1end_more_starts = filtered_group.groupby(['end']).filter(lambda g: g['start'].nunique() > 1)
  155. ex_1end_more_starts_count = ex_1end_more_starts.groupby(['end'])['transcript_id'].transform('count')
  156. composite_condition2 = ex_1end_more_starts[(ex_1end_more_starts_count < tr_num)]
  158. if '+' in group['strand'].unique():
  160. group.loc[group['start'].isin(ex_1start_more_ends['start']), 'splicing_event'] ='alt_5_prime'
  161. group.loc[group['end'].isin(ex_1end_more_starts['end']), 'splicing_event'] ='alt_3_prime'
  162. if len(composite_condition1)>0:
  163. group.loc[group['start'].isin(composite_condition1['start']), 'splicing_event'] ='composite_alt_5_skipped'
  164. if len(composite_condition2)>0:
  165. group.loc[group['end'].isin(composite_condition2['end']), 'splicing_event'] ='composite_alt_3_skipped'
  166. group.loc[group.apply(lambda row: (row['start'], row['end']) in skipped_tuples and pd.isna(row['splicing_event']), axis=1),
  167. 'splicing_event'] = 'skipped_exon'
  168. group.loc[group.apply(lambda row: ((row['start'], row['end']) in constitutive_tuples) and ((row['start'], row['end']) in constitutive_alltrs_tuples), axis=1),
  169. 'splicing_event'] = 'constitutive'
  170. group.loc[group.apply(lambda row: ((row['start'], row['end']) in constitutive_tuples) and not ((row['start'], row['end']) in constitutive_alltrs_tuples), axis=1),
  171. 'splicing_event'] = 'constitutive_totranscript'
  173. elif '-' in group['strand'].unique():
  174. group.loc[group['start'].isin(ex_1start_more_ends['start']), 'splicing_event'] ='alt_3_prime(-)'
  175. if len(composite_condition1)>0:
  176. group.loc[group['start'].isin(composite_condition1['start']), 'splicing_event'] ='composite_alt_3_skipped(-)'
  177. group.loc[group['end'].isin(ex_1end_more_starts['end']), 'splicing_event'] ='alt_5_prime(-)'
  178. if len(composite_condition2)>0:
  179. group.loc[group['end'].isin(composite_condition2['end']), 'splicing_event'] ='composite_alt_5_skipped(-)'
  180. group.loc[group.apply(lambda row: (row['start'], row['end']) in skipped_tuples and pd.isna(row['splicing_event']), axis=1),
  181. 'splicing_event'] = 'skipped_exon(-)'
  182. group.loc[group.apply(lambda row: ((row['start'], row['end']) in constitutive_tuples) and ((row['start'], row['end']) in constitutive_alltrs_tuples), axis=1),
  183. 'splicing_event'] = 'constitutive(-)'
  184. group.loc[group.apply(lambda row: ((row['start'], row['end']) in constitutive_tuples) and not ((row['start'], row['end']) in constitutive_alltrs_tuples), axis=1),
  185. 'splicing_event'] = 'constitutive_totranscript(-)'
  186. # This snippet is repeated so that some exons, which get labelled wrongly as any other events, even though they show those characteristics, to keep it simple, this will override and will always label such exons as first exons over any other thing.
  187. if '+' in group['strand'].unique():
  188. group.loc[smallest_exon_rows.index, 'splicing_event'] = 'first_exon'
  189. group.loc[greatest_exon_rows.index, 'splicing_event'] = 'last_exon'
  191. elif '-' in group['strand'].unique():
  192. group.loc[smallest_exon_rows.index, 'splicing_event'] = 'last_exon'
  193. group.loc[greatest_exon_rows.index, 'splicing_event'] = 'first_exon'
  195. gene_model = gene_model.merge(
  196. group[['start', 'end', 'splicing_event']],
  197. on=['start', 'end'],
  198. how='left'
  199. )
  200. gene_model['splicing_event'] = gene_model['splicing_event_y'].combine_first(gene_model['splicing_event_x'])
  201. gene_model.drop(['splicing_event_x', 'splicing_event_y'], axis=1, inplace=True)
  203. return gene_model
  205. def intron_filter(group):
  206. exon_condition = group[
  207. group['splicing_event'].astype(str).str.contains('constitutive') |
  208. group['splicing_event'].astype(str).str.contains('skipped_exon') |
  209. group['splicing_event'].astype(str).str.contains('first') |
  210. group['splicing_event'].astype(str).str.contains('last')]
  211. exon_range = set(exon_condition[['start','end']].itertuples(index=False, name=None))
  212. first_exons = group[group['splicing_event'].astype(str).str.contains('first')]
  213. last_exons = group[group['splicing_event'].astype(str).str.contains('last')]
  214. alt_3_exons = group[group['splicing_event'].astype(str).str.contains('alt_3')]
  215. alt_5_exons = group[group['splicing_event'].astype(str).str.contains('alt_5')]
  216. skipped_exons = group[group['splicing_event'].astype(str).str.contains('skipped_exon')]
  217. skipped_range = set(skipped_exons[['start', 'end']].itertuples(index=False, name=None))
  218. constitutive = group[group['splicing_event'].astype(str).str.contains('constitutive')]
  220. if '+' in group['strand'].unique():
  221. limitcoord_alt3 = alt_3_exons.groupby('end').apply(
  222. lambda group: group.loc[group['start'].idxmin()]
  223. ).reset_index(drop=True)
  225. limitcoord_alt5 = alt_5_exons.groupby('start').apply(
  226. lambda group: group.loc[group['end'].idxmax()]
  227. ).reset_index(drop=True)
  229. elif '-' in group['strand'].unique():
  230. limitcoord_alt5 = alt_5_exons.groupby('end').apply(
  231. lambda group: group.loc[group['start'].idxmin()]
  232. ).reset_index(drop=True)
  234. limitcoord_alt3 = alt_3_exons.groupby('start').apply(
  235. lambda group: group.loc[group['end'].idxmax()]
  236. ).reset_index(drop=True)
  238. group['correct_coords'] = np.nan
  239. introns = group[group['group']=='intron'].copy()
  241. for index, row in introns.iterrows():
  242. start, end, strand = row['start'], row['end'], row['strand']
  243. tr_id = row['transcript_id']
  244. is_in_exon_range = any(
  245. (exon_start <= end and exon_end >= start)
  246. for exon_start, exon_end in exon_range)
  248. if strand == '+':
  249. matches_first = ((start in (first_exons['end'] + 1).values) or (end in (first_exons['start'] - 1).values)) & (tr_id in first_exons['transcript_id'])
  251. matches_last = ((start in (last_exons['end'] + 1).values) or (end in (last_exons['start'] - 1).values)) & (tr_id in last_exons['transcript_id'])
  253. matches_alt3 = not limitcoord_alt3.empty and (
  254. end in (limitcoord_alt3['start'] - 1).values)
  256. matches_alt5 = not limitcoord_alt5.empty and (
  257. start in (limitcoord_alt5['end'] + 1).values)
  259. matches_constitutive = not constitutive.empty and (
  260. start in (constitutive['end'] + 1).values or end in (constitutive['start'] - 1).values)
  262. matches_skipped = not skipped_exons.empty and (
  263. start in (skipped_exons['end'] + 1).values or end in (skipped_exons['start'] - 1).values)
  265. elif (strand == '-'):
  266. matches_first = ((start in (first_exons['end'] + 1).values) or (end in (first_exons['start'] - 1).values)) & (tr_id in first_exons['transcript_id'])
  268. matches_last = ((start in (last_exons['end'] + 1).values) or (end in (last_exons['start'] - 1).values)) & (tr_id in last_exons['transcript_id'])
  270. matches_alt3 = not limitcoord_alt3.empty and (
  271. start in (limitcoord_alt3['end'] + 1).values)
  273. matches_alt5 = not limitcoord_alt5.empty and (
  274. end in (limitcoord_alt5['start'] - 1).values)
  276. matches_constitutive = not constitutive.empty and (
  277. start in (constitutive['end'] + 1).values or end in (constitutive['start'] - 1).values)
  279. matches_skipped = not skipped_exons.empty and (
  280. start in (skipped_exons['end'] + 1).values or end in (skipped_exons['start'] - 1).values)
  283. if (matches_alt3 or matches_alt5 or matches_constitutive or matches_skipped or matches_first or matches_last) and not is_in_exon_range:
  284. group.at[index, 'correct_coords'] = 'True'
  285. else:
  286. group.at[index, 'correct_coords'] = None
  288. for index, row in group.iterrows():
  289. type = row['group']
  290. if type == 'intron' and pd.isna(row['correct_coords']):
  291. group.drop(index, inplace=True)
  293. return group
  296. def intron_cat(group):
  297. introns = group[group['group']=='intron']
  298. first_int = group[group['splicing_event'].str.contains('first', na = False)][['start','end','strand','transcript_id']]
  299. last_int = group[group['splicing_event'].str.contains('last', na = False)][['start','end','strand','transcript_id']]
  301. for index, row in introns.iterrows():
  302. if (row['strand']=='+'):
  303. if (row['start'] in first_int['end'].values + 1) and (row['transcript_id'] in first_int['transcript_id'].values):
  304. current_val = group.at[index, 'splicing_event']
  305. group.at[index, 'splicing_event'] = f"{current_val}, first_intron" if pd.notna(current_val) else "first_intron"
  306. elif (row['end'] in last_int['start'].values - 1) and (row['transcript_id'] in last_int['transcript_id'].values):
  307. current_val = group.at[index, 'splicing_event']
  308. group.at[index, 'splicing_event'] = f"{current_val}, last_intron" if pd.notna(current_val) else "last_intron"
  309. elif (row['strand']=='-'):
  310. if (row['end'] in first_int['start'].values - 1) and (row['transcript_id'] in first_int['transcript_id'].values):
  311. current_val = group.at[index, 'splicing_event']
  312. group.at[index, 'splicing_event'] = f"{current_val}, first_intron(-)" if pd.notna(current_val) else "first_intron(-)"
  313. elif (row['start'] in last_int['end'].values + 1) and (row['transcript_id'] in last_int['transcript_id'].values):
  314. current_val = group.at[index, 'splicing_event']
  315. group.at[index, 'splicing_event'] = f"{current_val}, last_intron(-)" if pd.notna(current_val) else "last_intron(-)"
  317. constex_1 = group[group['splicing_event'].str.contains('constitutive',na=False)][['start','end','strand']]
  319. # Labels introns flanking constitutive exons on both strands
  320. for index, row in introns.iterrows():
  321. if (row['start'] in constex_1['end'].values + 1):
  322. if (row['strand']=='+'):
  323. group.at[index, 'splicing_event'] = 'constitutive_downstream'
  324. elif (row['strand']=='-'):
  325. group.at[index, 'splicing_event'] = 'constitutive_upstream(-)'
  327. elif (row['end'] in constex_1['start'].values - 1):
  328. if (row['strand']=='+'):
  329. group.at[index, 'splicing_event'] = 'constitutive_upstream'
  330. elif (row['strand']=='-'):
  331. group.at[index, 'splicing_event'] = 'constitutive_downstream(-)'
  333. constex = group[group['splicing_event'].str.contains('constitutive_totranscript',na=False)][['start','end','strand']]
  335. # Labels introns flanking constitutive exons grouped by first exon coords, on both strands
  336. for index, row in introns.iterrows():
  337. if (row['start'] in constex['end'].values + 1):
  338. if (row['strand']=='+'):
  339. group.at[index, 'splicing_event'] = 'constitutive_totranscript_downstream'
  340. elif (row['strand']=='-'):
  341. group.at[index, 'splicing_event'] = 'constitutive_totranscript_upstream(-)'
  343. elif (row['end'] in constex['start'].values - 1):
  344. if (row['strand']=='+'):
  345. group.at[index, 'splicing_event'] = 'constitutive_totranscript_upstream'
  346. elif (row['strand']=='-'):
  347. group.at[index, 'splicing_event'] = 'constitutive_totranscript_downstream(-)'
  349. # Labels introns flanking skipped exons on both strands
  350. skipex = group[group['splicing_event'].str.contains('skipped_exon',na=False)][['start','end','strand']]
  351. for index, row in introns.iterrows():
  352. if (row['start'] in skipex['end'].values + 1):
  353. if (row['strand']=='+'):
  354. current_val = group.at[index, 'splicing_event']
  355. if pd.notna(current_val):
  356. group.at[index, 'splicing_event'] = current_val + ', skipped_downstream'
  357. else:
  358. group.at[index, 'splicing_event'] = 'skipped_downstream'
  359. elif (row['strand']=='-'):
  360. current_val = group.at[index, 'splicing_event']
  361. if pd.notna(current_val):
  362. group.at[index, 'splicing_event'] = current_val + ', skipped_upstream(-)'
  363. else:
  364. group.at[index, 'splicing_event'] = 'skipped_upstream(-)'
  366. elif (row['end'] in skipex['start'].values - 1):
  367. if (row['strand']=='+'):
  368. current_val = group.at[index, 'splicing_event']
  369. if pd.notna(current_val):
  370. group.at[index, 'splicing_event'] = current_val + ', skipped_upstream'
  371. else:
  372. group.at[index, 'splicing_event'] = 'skipped_upstream'
  373. elif (row['strand']=='-'):
  374. current_val = group.at[index, 'splicing_event']
  375. if pd.notna(current_val):
  376. group.at[index, 'splicing_event'] = current_val + ', skipped_downstream(-)'
  377. else:
  378. group.at[index, 'splicing_event'] = 'skipped_downstream(-)'
  380. #Labels alternative 5'ss exon flanking introns on both strands:
  381. alt_5 = group[group['splicing_event'].str.contains('alt_5_prime',na=False)][['gene_id','group','start','end','strand','transcript_id']]
  382. alt_5_ex = alt_5[alt_5['group']=='exon']
  383. if '+' in alt_5['strand'].unique():
  384. alt5_ends = alt_5_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  385. max_end_indices = alt5_ends.groupby(['gene_id', 'start'])['end'].idxmax() #Now I am only selecting the start and end coordinates which are #containing a greater end position i.e. the second 5'ss Because this is where the intron truly begins.
  386. alt5_grt_end = alt5_ends.loc[max_end_indices]
  387. elif '-' in alt_5['strand'].unique():
  388. alt_5_starts = alt_5_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  389. min_start_indices = alt_5_starts.groupby(['gene_id', 'end'])['start'].idxmin()
  390. neg_alt_small_starts = alt_5_starts.loc[min_start_indices]
  392. for index, row in introns.iterrows():
  393. if (row['strand']=='+') and not alt_5.empty:
  394. if (row['start'] in alt5_grt_end['end'].values + 1):
  395. current_val = group.at[index, 'splicing_event']
  396. if pd.notna(current_val):
  397. group.at[index, 'splicing_event'] = current_val + ', alt_5_downstream'
  398. else:
  399. group.at[index, 'splicing_event'] = 'alt_5_downstream'
  400. elif (row['end'] in alt5_grt_end['start'].values - 1):
  401. current_val = group.at[index, 'splicing_event']
  402. if pd.notna(current_val):
  403. group.at[index, 'splicing_event'] = current_val + ', alt_5_upstream'
  404. else:
  405. group.at[index, 'splicing_event'] = 'alt_5_upstream'
  407. elif (row['strand']=='-') and not alt_5.empty:
  408. if (row['start'] in neg_alt_small_starts['end'].values + 1):
  409. current_val = group.at[index, 'splicing_event']
  410. if pd.notna(current_val):
  411. group.at[index, 'splicing_event'] = current_val + ', alt_5_upstream(-)'
  412. else:
  413. group.at[index, 'splicing_event'] = 'alt_5_upstream(-)'
  414. elif (row['end'] in neg_alt_small_starts['start'].values - 1):
  415. current_val = group.at[index, 'splicing_event']
  416. if pd.notna(current_val):
  417. group.at[index, 'splicing_event'] = current_val + ', alt_5_downstream(-)'
  418. else:
  419. group.at[index, 'splicing_event'] = 'alt_5_downstream(-)'
  421. #Labelling composite alternative 5'ss skipped exon flanking introns
  422. composite_alt_5 = group[group['splicing_event'].str.contains('composite_alt_5_skipped',na=False)][['gene_id','group','start','end','strand','transcript_id']]
  423. composite_alt_5_ex = composite_alt_5[composite_alt_5['group']=='exon']
  424. if '+' in composite_alt_5['strand'].unique():
  425. comp_alt5_ends = composite_alt_5_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  426. comp_max_end_indices = comp_alt5_ends.groupby(['gene_id', 'start'])['end'].idxmax() #Now I am only selecting the start and end coordinates which are #containing a greater end position i.e. the second 5'ss Because this is where the intron truly begins.
  427. comp_alt5_grt_end = comp_alt5_ends.loc[comp_max_end_indices]
  428. elif '-' in composite_alt_5['strand'].unique():
  429. composite_alt_5_starts = composite_alt_5_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  430. comp_min_start_indices = composite_alt_5_starts.groupby(['gene_id', 'end'])['start'].idxmin()
  431. comp_neg_alt_small_starts = composite_alt_5_starts.loc[comp_min_start_indices]
  433. for index, row in introns.iterrows():
  434. if (row['strand']=='+') and not composite_alt_5.empty:
  435. if (row['start'] in comp_alt5_grt_end['end'].values + 1):
  436. current_val = group.at[index, 'splicing_event']
  437. if pd.notna(current_val):
  438. group.at[index, 'splicing_event'] = current_val + ', composite_alt_5_downstream'
  439. else:
  440. group.at[index, 'splicing_event'] = 'composite_alt_5_downstream'
  441. elif (row['end'] in comp_alt5_grt_end['start'].values - 1):
  442. current_val = group.at[index, 'splicing_event']
  443. if pd.notna(current_val):
  444. group.at[index, 'splicing_event'] = current_val + ', composite_alt_5_upstream'
  445. else:
  446. group.at[index, 'splicing_event'] = 'composite_alt_5_upstream'
  448. elif (row['strand']=='-') and not composite_alt_5.empty:
  449. if (row['start'] in comp_neg_alt_small_starts['end'].values + 1):
  450. current_val = group.at[index, 'splicing_event']
  451. if pd.notna(current_val):
  452. group.at[index, 'splicing_event'] = current_val + ', composite_alt_5_upstream(-)'
  453. else:
  454. group.at[index, 'splicing_event'] = 'composite_alt_5_upstream(-)'
  455. elif (row['end'] in comp_neg_alt_small_starts['start'].values - 1):
  456. current_val = group.at[index, 'splicing_event']
  457. if pd.notna(current_val):
  458. group.at[index, 'splicing_event'] = current_val + ', composite_alt_5_downstream(-)'
  459. else:
  460. group.at[index, 'splicing_event'] = 'composite_alt_5_downstream(-)'
  462. #Labels alternative 5'ss exon flanking introns on both strands:
  463. alt_3 = group[group['splicing_event'].str.contains('alt_3_prime',na=False)][['gene_id','group','start','end','strand','transcript_id']]
  464. alt_3_ex = alt_3[alt_3['group']=='exon']
  465. if '+' in alt_3['strand'].unique():
  466. alt3_starts = alt_3_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  467. min_start_indices = alt3_starts.groupby(['gene_id', 'end'])['start'].idxmin()
  468. alt3_small_starts = alt3_starts.loc[min_start_indices]
  469. elif '-' in alt_3['strand'].unique():
  470. N_alt3_ends = alt_3_ex.groupby(['gene_id', 'start', 'end'])['transcript_id'].nunique().reset_index()
  471. max_end_indices = N_alt3_ends.groupby(['gene_id', 'start'])['end'].idxmax()
  472. N_alt3_grt_end = N_alt3_ends.loc[max_end_indices]
  475. for index, row in introns.iterrows():
  476. if (row['strand']=='+') and not alt_3.empty:
  477. if (row['start'] in alt3_small_starts['end'].values + 1):
  478. current_val = group.at[index, 'splicing_event']
  479. if pd.notna(current_val):
  480. group.at[index, 'splicing_event'] = current_val + ', alt_3_downstream'
  481. else:
  482. group.at[index, 'splicing_event'] = 'alt_3_downstream'
  483. elif (row['end'] in alt3_small_starts['start'].values - 1):
  484. current_val = group.at[index, 'splicing_event']
  485. if pd.notna(current_val):
  486. group.at[index, 'splicing_event'] = current_val + ', alt_3_upstream'
  487. else:
  488. group.at[index, 'splicing_event'] = 'alt_3_upstream'
  490. elif (row['strand']=='-') and not alt_3.empty:
  491. if (row['start'] in N_alt3_grt_end['end'].values + 1):
  492. current_val = group.at[index, 'splicing_event']
  493. if pd.notna(current_val):
  494. group.at[index, 'splicing_event'] = current_val + ', alt_3_upstream(-)'
  495. else:
  496. group.at[index, 'splicing_event'] = 'alt_3_upstream(-)'
  497. elif (row['end'] in N_alt3_grt_end['start'].values - 1):
  498. current_val = group.at[index, 'splicing_event']
  499. if pd.notna(current_val):
  500. group.at[index, 'splicing_event'] = current_val + ', alt_3_downstream(-)'
  501. else:
  502. group.at[index, 'splicing_event'] = 'alt_3_downstream(-)'
  504. #Labelling composite alternative 3'ss skipped exon flanking introns
  505. composite_alt_3 = group[group['splicing_event'].str.contains('composite_alt_3_skipped',na=False)][['gene_id','group','start','end','strand','transcript_id']]
  506. composite_alt_3_ex = composite_alt_3[composite_alt_3['group']=='exon']
  507. if '+' in composite_alt_3['strand'].unique():
  508. comp_alt3_starts = composite_alt_3_ex.groupby(['gene_id', 'end', 'start'])['transcript_id'].nunique().reset_index()
  509. comp_min_start_indices = comp_alt3_starts.groupby(['gene_id', 'end'])['start'].idxmin()
  510. comp_alt3_small_starts = comp_alt3_starts.loc[comp_min_start_indices]
  511. elif '-' in composite_alt_3['strand'].unique():
  512. comp_N_alt3_ends = composite_alt_3_ex.groupby(['gene_id', 'start', 'end'])['transcript_id'].nunique().reset_index()
  513. comp_max_end_indices = comp_N_alt3_ends.groupby(['gene_id', 'start'])['end'].idxmax()
  514. comp_N_alt3_grt_end = comp_N_alt3_ends.loc[comp_max_end_indices]
  517. for index, row in introns.iterrows():
  518. if (row['strand']=='+') and not composite_alt_3.empty:
  519. if (row['start'] in comp_alt3_small_starts['end'].values + 1):
  520. current_val = group.at[index, 'splicing_event']
  521. if pd.notna(current_val):
  522. group.at[index, 'splicing_event'] = current_val + ', composite_alt_3_downstream'
  523. else:
  524. group.at[index, 'splicing_event'] = 'composite_alt_3_downstream'
  525. elif (row['end'] in comp_alt3_small_starts['start'].values - 1):
  526. current_val = group.at[index, 'splicing_event']
  527. if pd.notna(current_val):
  528. group.at[index, 'splicing_event'] = current_val + ', composite_alt_3_upstream'
  529. else:
  530. group.at[index, 'splicing_event'] = 'composite_alt_3_upstream'
  532. elif (row['strand']=='-') and not composite_alt_3.empty:
  533. if (row['start'] in comp_N_alt3_grt_end['end'].values + 1):
  534. current_val = group.at[index, 'splicing_event']
  535. if pd.notna(current_val):
  536. group.at[index, 'splicing_event'] = current_val + ', composite_alt_3_upstream(-)'
  537. else:
  538. group.at[index, 'splicing_event'] = 'composite_alt_3_upstream(-)'
  539. elif (row['end'] in comp_N_alt3_grt_end['start'].values - 1):
  540. current_val = group.at[index, 'splicing_event']
  541. if pd.notna(current_val):
  542. group.at[index, 'splicing_event'] = current_val + ', composite_alt_3_downstream(-)'
  543. else:
  544. group.at[index, 'splicing_event'] = 'composite_alt_3_downstream(-)'
  546. return group
  549. keywords = []
  551. import pandas as pd
  553. r2 = pd.read_csv(gff_file, sep='\t', names=['Chr', 'source', 'group', 'start', 'end', '_', 'strand', '-', 'info'])
  554. r2.dropna(axis=0, inplace=True)
  555. r2_Ex_Int = r2[(r2['group'] == 'exon') | (r2['group'] == 'intron')].copy()
  557. # r2_Ex_Int[['ID','Parent','gene_id', 'gene_symbol','transcript_id','transcript_symbol']] = r2_Ex_Int[
  558. # 'info'].str.split(';', expand=True) #- Drosophila gff has these entries in metadata column
  560. print('\n \\ The original gff dataframe. Parsing out the metadata column - info. \\\n')
  561. print(r2_Ex_Int.head(10))
  563. with open(keywords_file, 'r') as file:
  564. for w in file:
  565. keywords.append(str(w.rstrip()))
  567. # Converting metadata column identifiers (keywords) into columns and extracting relevant information.
  568. pattern = '|'.join(r'(?P<{}>{}=[^;]+)'.format(k, k) for k in keywords)
  569. extracted_info = r2_Ex_Int['info'].str.extractall(pattern)
  570. grouped_info = extracted_info.fillna('').astype(str).groupby(level=0).agg(lambda x: ';'.join(filter(None, x)))
  571. grouped_info.columns = [f'{col}' for col in grouped_info.columns]
  572. r2_Ex_Int = r2_Ex_Int.join(grouped_info, how='left')
  574. print('\n \\ Here\'s a glimpse of the modified gff dataframe. \\ \n')
  575. print(r2_Ex_Int.head(10))
  576. r2_Ex_Int.to_csv('gff_parsed_meta.csv', sep='\t')
  579. r2_Ex_Int.sort_values(by=['gene_id', 'start'], ascending=True, inplace=True)
  580. r2_Ex_Int.drop(labels='info', axis=1, inplace=True)
  581. r2_Ex_Int.loc[:, 'sort_order'] = r2_Ex_Int.apply(lambda x: x['start'] if x['strand'] == '+' else -x['start'], axis=1)
  582. r2_Ex_Int = r2_Ex_Int.sort_values(by=['sort_order'])
  583. r2_Ex_Int = r2_Ex_Int.drop(columns=['sort_order'])
  585. r2_Ex_Int['order_of_appearance'] = r2_Ex_Int.groupby(['group', 'transcript_id']).cumcount() + 1
  587. r2_Ex_Int_1 = r2_Ex_Int.groupby(['gene_id']).apply(lambda x: retained_introns(x))
  588. r2_Ex_Int_1.reset_index(drop=True, inplace=True)
  589. r2_Exons_annot_composite = r2_Ex_Int_1.groupby(['gene_id']).apply(lambda x: exon_cat(x))
  590. r2_Exons_annot_composite.reset_index(drop=True, inplace=True)
  592. retained_introns = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'retained intron'])
  593. constitutive = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'constitutive']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'constitutive(-)'])
  594. skipped = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'skipped_exon']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'skipped_exon(-)'])
  595. alt_5 = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'alt_5_prime']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'alt_5_prime(-)'])
  596. alt_3 = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'alt_3_prime']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'alt_3_prime(-)'])
  597. composite = len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'composite_alt_3_skipped']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'composite_alt_3_skipped(-)']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'composite_alt_5_skipped']) + len(r2_Exons_annot_composite[r2_Exons_annot_composite['splicing_event'] == 'composite_alt_5_skipped(-)'])
  599. print(f'Total number of constitutive exons: {constitutive}')
  600. print(f'Total number of skipped exons: {skipped}')
  601. print(f'Total number of alternative 5\'ss exons: {alt_5}')
  602. print(f'Total number of alternative 3\'ss exons: {alt_3}')
  603. print(f'Total number of composite alternative 5\'ss or alternative 3\'ss containing and skipped exons: {composite}')
  604. print(f'Total number of retained introns: {retained_introns}')
  605. df_final = r2_Exons_annot_composite.groupby(['gene_id']).apply(lambda x: intron_filter(x))
  606. df_final.reset_index(drop=True, inplace=True)
  607. all_annot_df = df_final.groupby(['gene_id']).apply(lambda x: intron_cat(x))
  609. df_final.to_csv('./Exonslabelled_intronsnot.csv', sep='\t')
  611. print(all_annot_df.head(10))
  613. all_annot_df.reset_index(drop=True, inplace=True)
  615. all_annot_df.sort_values(by=['gene_id', 'transcript_id', 'order_of_appearance'], ascending=True, inplace=True)
  617. ############################################################# The functions below, scaled_coordinates () and scale_group () were written for cases where we only have sequence fasta files of orthologues and reference sps
  618. # For example, the analysis I did on C.elegans switch like splicing events. I had relative coordinates that were inferred from the fasta files. I didn't have absolute coordinates ######## If this isn't the case, the lines
  619. # where these functions are being called can be commented out. ##### The lines are only to represent everything in terms of sclaed or relative coordinates.#################################################################
  620. if coordinate_type == 'relative':
  621. def scale_coordinates(row, start, end):# removed strand from args, since same function works for both strands.
  622. max_length = start - end
  623. qstart = 0
  624. scaled_rstart = qstart + (row['start'] - start)
  625. scaled_rend = qstart + (row['end'] - start)
  626. return pd.Series([scaled_rstart, scaled_rend], index=['scaled_rstart', 'scaled_rend'])
  628. # Group by 'gene' and apply scaling function
  629. def scale_group(group):
  630. # For each group (gene), calculate min and max start and end, works for both strands since start and ends are always
  631. # smaller start and greater end, only the order in negative strand is reversed, i.e., exon with smallest start is the
  632. # last exon. Draw it out to avoid confusion.
  633. strand = group['strand'].iloc[0]
  634. min_rstart = group['start'].min()
  635. max_rend = group['end'].max()
  636. group[['scaled_rstart', 'scaled_rend']] = group.apply(scale_coordinates, axis=1, args=(min_rstart, max_rend))
  638. group = group.rename(columns={'start': 'absolute_start', 'end': 'absolute_end'})
  639. group = group.rename(columns={'scaled_rstart': 'start', 'scaled_rend': 'end'})
  640. return group
  642. # Assuming 'Parent' identifies different genes
  643. scaled_df = all_annot_df.groupby('Parent').apply(scale_group).reset_index(drop=True)
  644. # Assuming 'Parent' identifies different genes
  645. all_annot_Introns = scaled_df[scaled_df['group'] == 'intron']
  647. elif coordinate_type =='absolute':
  648. all_annot_Introns = all_annot_df[all_annot_df['group'] == 'intron']
  650. all_annot_Introns = all_annot_Introns.copy()
  651. #Unknown_Intron_indices = all_annot_Introns[all_annot_Introns['splicing_event_y_y'] == 'intron'].index
  652. #all_annot_Introns.drop(Unknown_Intron_indices, inplace=True)
  653. all_annot_Introns['ID'] = all_annot_Introns['Chr'].astype(str) + '_' + all_annot_Introns['start'].astype(str) + '_' + \
  654. all_annot_Introns['end'].astype(str)
  655. all_annot_Introns.groupby(['ID', 'gene_id', 'splicing_event'])['splicing_event'].nunique()
  656. Intron_df = all_annot_Introns
  658. # 1. Calculating size of introns
  659. all_annot_Introns['size'] = abs(all_annot_Introns['start'] - all_annot_Introns['end'])
  660. print('The median size of introns is : ')
  661. print(all_annot_Introns['size'].median())
  663. median_size = all_annot_Introns['size'].median()
  666. def round_to_nearest_multiple_of_10(number):
  667. return round(number / 10) * 10
  669. if median == 0:
  670. median = round_to_nearest_multiple_of_10(median_size)
  672. print('The introns bigger than the median size are being resized to - %d bp. \n' % median)
  675. # 2. Eliminating extremely small introns
  676. small_introns = all_annot_Introns[all_annot_Introns['size'] < 10].index
  677. all_annot_Introns = all_annot_Introns.drop(small_introns)
  679. # 3. Creating a separate dataframe for introns smaller than median bp in length
  680. Intron_size_median = all_annot_Introns[all_annot_Introns['size'] < median]
  682. # 4. Processing introns greater than median bp in length
  683. Intron_size_great_median = all_annot_Introns[all_annot_Introns['size'] >= median]
  684. Intron_size_great_median['End_medianbp_apart'] = Intron_size_great_median['start'] + median
  685. Intron_size_great_median['Start_medianbp_apart'] = Intron_size_great_median['end'] - median
  687. # Creating new ID columns for starting median bp and ending median bp of introns
  688. Intron_size_great_median['ID_Starting_median'] = Intron_size_great_median['Chr'].astype(str) + '_' + \
  689. Intron_size_great_median[
  690. 'strand'] + '_' + Intron_size_great_median['start'].astype(
  691. str) + '_' + Intron_size_great_median[
  692. 'End_medianbp_apart'].astype(str)
  694. Intron_size_great_median['ID_Ending_median'] = Intron_size_great_median['Chr'].astype(str) + '_' + \
  695. Intron_size_great_median[
  696. 'strand'] + '_' + Intron_size_great_median[
  697. 'Start_medianbp_apart'].astype(str) + '_' + Intron_size_great_median[
  698. 'end'].astype(str)
  700. merged_intron_df = pd.concat([Intron_size_great_median, Intron_size_median], ignore_index=True)
  701. merged_intron_df.drop_duplicates(inplace=True)
  702. merged_intron_df.to_csv(intron_annot_file, sep='\t')
  704. # Extracting the coordinates and arranging them based on their strands. ##
  705. # NOTE: THESE NEED TO BE CONVERTED TO BED FILE ###
  706. # CAN BE CONVERTED TO BED FILE EASILY WITH EXCEL AND R - GTF-to-bed.r###
  707. intron_beginningmedian = Intron_size_great_median[
  708. ['Parent','Chr','group', 'start', 'End_medianbp_apart', 'strand', 'ID_Starting_median']]
  709. intron_endingmedian = Intron_size_great_median[['Parent','Chr','group', 'Start_medianbp_apart', 'end', 'strand', 'ID_Ending_median']]
  710. small_intron_coordinates = Intron_size_median[['Parent','Chr','group', 'start', 'end', 'strand', 'ID']]
  711. intron_beginningmedian.rename(columns={'End_medianbp_apart': 'end', 'ID_Starting_median': 'ID'}, inplace=True)
  712. intron_endingmedian.rename(columns={'Start_medianbp_apart': 'start', 'ID_Ending_median': 'ID'}, inplace=True)
  713. big_introns = pd.concat([intron_beginningmedian, intron_endingmedian], ignore_index=True)
  714. all_intron_coordinates = pd.concat([small_intron_coordinates, big_introns], ignore_index=True)
  715. all_intron_coordinates.drop_duplicates(inplace=True)
  716. all_intron_coordinates.to_csv(intron_cords_file, sep='\t')
  719. if coordinate_type == 'relative':
  720. all_exon_coordinates = scaled_df[scaled_df['group']=='exon'][['gene_id','Parent','Chr','group','start', 'end','strand', 'ID']]
  721. all_coordinates = pd.concat([all_intron_coordinates, all_exon_coordinates], ignore_index=True)
  722. all_coordinates.sort_values(by=['gene_id', 'start'], ascending=[True, True], inplace=True)
  724. all_coordinates.to_csv('./CelsWS15_Exon_IntronFragment_coordinates_April2025.tsv', sep='\t', index=False)
  725. scaled_df.to_csv(exon_intron_annot_file, sep='\t')
  727. elif coordinate_type == 'absolute':
  728. all_exon_coordinates = all_annot_df[all_annot_df['group']=='exon'][['gene_id','Parent','Chr','group','start', 'end','strand', 'ID']]
  729. all_coordinates = pd.concat([all_intron_coordinates, all_exon_coordinates], ignore_index=True)
  730. all_coordinates.sort_values(by=['gene_id', 'start'], ascending=[True, True], inplace=True)
  732. all_coordinates.to_csv('./Exon_IntronFragment_April2025_coordinates.tsv', sep='\t', index=False)
  733. all_annot_df.to_csv(exon_intron_annot_file, sep='\t')
  734. # The end - this will give an annotated intron coordinates file which can be further analyzed in excel.###
  736. #%%

gff_exon-intron_annotations.py, under Apache-2.0 · at the source

Overview

Authors: Sanjana Bhatnagar1, Jade Ho1, Nour H Sadek1, Isha Singh1, Michael Zoberman1, Bina Koterniak1, Yufang Liu1, Daniel D Fusca2, Asher D Cutter2, Alan M Moses1, John A Calarco1
ORCID iDs: Sanjana Bhatnagar, Jade Ho, Daniel D Fusca, Asher D Cutter, Alan M Moses, John A Calarco
  1. Department of Cell and Systems Biology, 25 Harbord Street, Toronto, Ontario, M5S 3G5, Canada
  2. Department of Ecology and Evolutionary Biology, 25 Willcocks Street, Toronto, Ontario, M5S 3B2, Canada
Institutions: University of Toronto (Canada); Mayo Clinic (United States); Hospital for Sick Children (Canada)
Journal: Nucleic acids research, volume 54, issue 10, article gkag451
Dates: received 27 June 2025; accepted 22 April 2026; published online 25 May 2026; in print May 2026
Type: Research article · Language: English
License: CC BY-NC
Identifiers: DOI 10.1093/nar/gkag451 · PMID 42179040 · PMCID PMC13199693 · OpenAlex W4413385803
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.08.16.670696
Status: code verified
Categories: C. elegans (organism), cellular / molecular (subfield)
Methods: Statistics, Machine learning
MeSH: Alternative Splicing*, Caenorhabditis elegans*, Genes, Reporter*, Introns*, Animals, Enhancer Elements, Genetic, Muscle Cells, Neurons, Nucleotide Motifs, RNA-Binding Proteins, Silencer Elements, Transcriptional (* major topic)
Topic: RNA Research and Splicing (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: CIHR (180365, 156300); Canadian Institutes of Health Research (180365, 156300); Natural Sciences and Engineering Research Council of Canada (RGPIN-2017-06573)
Citations: not cited yet (Europe PMC); 137 references in the paper

Abstract

Introns play a critical role in regulating alternative splicing. However, identifying functional intronic motifs is challenging due to their short and degenerate sequence composition. Massively parallel reporter assays have provided insights into cis-regulatory logic governing alternative splicing, but these approaches are generally performed in cell culture, limiting their ability to capture tissue-specific contexts. Here, we implemented in vivo parallelized reporter assays (PRA) in Caenorhabditis elegans neurons and muscle cells to screen for intronic enhancer and silencer motifs among thousands of randomized sequences. We identified nearly 200 sequences regulating splicing in these tissues. We uncovered core sub-sequences with tissue-biased enhancing and silencing activity, including motifs recognized by well-characterized RNA-binding proteins, and unmapped motifs with no obvious cognate binding protein. Mapping our PRA-derived motifs to native introns flanking tissue-biased alternative exons revealed their conservation across nematodes, supporting their functional relevance. Additionally, individual intronic regions frequently contained diverse combinations of these motifs, indicative of complex engagement of these sequences by RNA-binding proteins. Finally, we performed targeted mutagenesis of PRA-derived intronic enhancers flanking a neuronal microexon, identifying key cis-regulatory determinants of microexon splicing. Together, our study provides a framework to explore the role of intronic regions in tissue-specific splicing regulation within a multicellular organism.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 1 match between paragraphs and lines of code.

Zenodo 15757714

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
4 files
At the source: doi.org/10.5281/zenodo.15757714

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • 1 match between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

The raw sequence data (fastq files) are available at the NCBI SRA:

https://www.ncbi.nlm.nih.gov/sra/PRJNA1359015

The code and additional data are available through Zenodo:

https://doi.org/10.5281/zenodo.17593183

Inferring Exon and Intron Metadata from .gff file (v1.0.1).

https://doi.org/10.5281/zenodo.15757714

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 11 MeSH terms, 3 funders, 134 references.

Cite

This paper

Bhatnagar, S., Ho, J., Sadek, N. H., Singh, I., Zoberman, M., Koterniak, B., Liu, Y., Fusca, D. D., Cutter, A. D., Moses, A. M., & Calarco, J. A. (2026). An in vivo parallelized reporter assay to uncover tissue-specific splicing regulatory sequences in a multicellular animal. Nucleic acids research, 54(10), gkag451. https://doi.org/10.1093/nar/gkag451

BibTeX

@article{bhatnagar2026vivo,
author = {Bhatnagar, Sanjana and Ho, Jade and Sadek, Nour H and Singh, Isha and Zoberman, Michael and Koterniak, Bina and Liu, Yufang and Fusca, Daniel D and Cutter, Asher D and Moses, Alan M and Calarco, John A},
title = {{An in vivo parallelized reporter assay to uncover tissue-specific splicing regulatory sequences in a multicellular animal}},
journal = {Nucleic acids research},
year = {2026},
month = may,
volume = {54},
number = {10},
pages = {gkag451},
publisher = {Oxford University Press},
issn = {0305-1048},
doi = {10.1093/nar/gkag451},
url = {https://doi.org/10.1093/nar/gkag451},
pmid = {42179040},
pmcid = {PMC13199693}
}

RIS

TY - JOUR
AU - Bhatnagar, Sanjana
AU - Ho, Jade
AU - Sadek, Nour H
AU - Singh, Isha
AU - Zoberman, Michael
AU - Koterniak, Bina
AU - Liu, Yufang
AU - Fusca, Daniel D
AU - Cutter, Asher D
AU - Moses, Alan M
AU - Calarco, John A
TI - An in vivo parallelized reporter assay to uncover tissue-specific splicing regulatory sequences in a multicellular animal
T2 - Nucleic acids research
J2 - Nucleic Acids Res
PY - 2026
DA - 2026/05/01
VL - 54
IS - 10
SP - gkag451
SN - 0305-1048
PB - Oxford University Press
DO - 10.1093/nar/gkag451
UR - https://doi.org/10.1093/nar/gkag451
LA - en
ER -

CSL-JSON

{
"id": "10.1093/nar/gkag451",
"type": "article-journal",
"title": "An in vivo parallelized reporter assay to uncover tissue-specific splicing regulatory sequences in a multicellular animal",
"container-title": "Nucleic acids research",
"author": [
{
"family": "Bhatnagar",
"given": "Sanjana"
},
{
"family": "Ho",
"given": "Jade"
},
{
"family": "Sadek",
"given": "Nour H"
},
{
"family": "Singh",
"given": "Isha"
},
{
"family": "Zoberman",
"given": "Michael"
},
{
"family": "Koterniak",
"given": "Bina"
},
{
"family": "Liu",
"given": "Yufang"
},
{
"family": "Fusca",
"given": "Daniel D"
},
{
"family": "Cutter",
"given": "Asher D"
},
{
"family": "Moses",
"given": "Alan M"
},
{
"family": "Calarco",
"given": "John A"
}
],
"container-title-short": "Nucleic Acids Res",
"volume": "54",
"issue": "10",
"page": "gkag451",
"DOI": "10.1093/nar/gkag451",
"PMID": "42179040",
"PMCID": "PMC13199693",
"ISSN": "0305-1048",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/nar/gkag451",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: pandas, Matplotlib, NumPy, 2 references, author John A Calarco
[2] doi:10.1038/s41467-026-76320-x [code]
Intrinsic and non-cell-autonomous roles for the neurodevelopmental syndrome-linked transcription factor UNC-3/EBF.
Journal: Nature communications
In common: C. elegans, 4 references
[3] doi:10.1186/s13024-026-00944-2
TDP-43: [GU]-ardian of the transcriptome.
Journal: Molecular neurodegeneration
In common: cellular / molecular, 4 references
[4] doi:10.1038/s41467-026-75352-7 [code]
Mechanosensory encoding of surface mechanics optimizes locomotion.
Journal: Nature communications
In common: pandas, Matplotlib, NumPy, C. elegans, cellular / molecular, 2 references
[5] doi:10.1126/sciadv.aec9329
Cohesin and NuRD antagonistically drive alternative neuronal fates via PLZF transcription factors.
Journal: Science advances
In common: C. elegans, cellular / molecular, 4 references
[6] doi:10.1038/s41586-026-10501-y [code]
Long-term editing of brain circuits using an engineered electrical synapse.
Journal: Nature
In common: pandas, Matplotlib, NumPy, C. elegans, cellular / molecular, 2 references
[7] doi:10.1016/j.celrep.2026.117288
Interlocked transcription factor feedback loops maintain and restore adult touch sensation.
Journal: Cell reports
In common: cellular / molecular, 4 references
[8] doi:10.1371/journal.pcbi.1013441 [code]
Large vision model framework for automated C. elegans analysis: From static morphometry to dynamic neural activity.
Journal: PLoS computational biology
In common: pandas, Matplotlib, NumPy, C. elegans, 2 references
[9] doi:10.1016/j.xhgg.2026.100629 [code]
Positive selection on brain cis-regulatory elements in the human lineage drives changes in gene expression and susceptibility to neuropsychiatric disorders.
Journal: HGG advances
In common: pandas, Matplotlib, NumPy, cellular / molecular, 2 references
[10] doi:10.1038/s41586-026-10348-3 [code]
An enteric neuron ionotropic receptor regulates salt stress resistance.
Journal: Nature
In common: pandas, Matplotlib, NumPy, C. elegans, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.