scmkl.get_region_groupings

  1import pandas as pd
  2import numpy as np
  3import re
  4
  5
  6def _find_overlap(start1 : int, end1 : int, start2 : int, end2 : int) -> bool:
  7    """
  8    Function to determine whether two regions on the same chromosome 
  9    overlap.
 10    
 11    Parameters
 12    ----------
 13    start1 : int
 14        The start position for region 1
 15    
 16    end1 : int
 17        The end position for region 1
 18    
 19    start2 : int
 20        The start position for region 2
 21    
 22    end2: int
 23        The end postion for region 2
 24    
 25    Returns
 26    -------
 27    is_overlapping : bool
 28        `True` if the regions overlap by 1bp. `False` if the regions 
 29        do not overlap
 30    """
 31    return max(start1, start2) <= min(end1, end2)
 32
 33
 34def get_tss(row : pd.DataFrame) -> int:
 35    """
 36    Takes a row from a DataFrame as a DataFrame with columns 
 37    ['start', 'end', 'strand'] and returns the transcription start site
 38    depending on gene strandedness.
 39    
 40    Parameters
 41    ----------
 42    row : pd.DataFrame 
 43        A gtf dataframe containing only one row and columns 
 44        `['start', 'end', 'strand']`.
 45    
 46    Returns
 47    -------
 48    tss : int
 49        The transcription start site for row's respective annotation.
 50    """
 51    if row.iloc[2] == '+':
 52        return row.iloc[0]
 53    
 54    elif row.iloc[2] == '-':
 55        return row.iloc[1]
 56    
 57    else:
 58        raise ValueError("Strand symbol must be '+' or '-'")
 59    
 60
 61def calc_range(row : pd.DataFrame, len_up : int, len_down : int) -> list:
 62    """
 63    Returns an infered promotor region for given annotation range 
 64    depending on transcription start site and user-defined upstream 
 65    and downstream adjustments.
 66
 67    Parameters
 68    ----------
 69    row : pd.DataFrame 
 70        A gtf as a dataframe containing only one row and columns 
 71        `['tss', 'strand']` where `'tss'` is the transcription start 
 72        site.
 73    
 74    len_up : int 
 75        Number of base pairs upstream from the transcription 
 76        start site the promotor range should be adjusted for.
 77
 78    len_down : int
 79        Number of base pairs downstream from the transcription start 
 80        site the promotor range should be adjusted for.
 81    
 82    Returns
 83    -------
 84    start_end : list
 85        A 2 element list where the first element is the adjusted start
 86        position and the second the the adjusted end position.
 87    """
 88    if row.iloc[1] == '+':
 89        start = row.iloc[0] - len_up
 90        end = row.iloc[0] + len_down
 91
 92    elif row.iloc[1] == '-':
 93        start = row.iloc[0] - len_down
 94        end = row.iloc[0] + len_up
 95
 96    else:
 97        raise ValueError('Strand symbol must be "+" or "-"')
 98    
 99    return [start, end]
100
101
102def adjust_regions(gene_anno : pd.DataFrame, len_up : int, len_down : int
103                    ) -> pd.DataFrame:
104    """
105    Takes a GTF file as a pd.DataFrame and adjusts start and end 
106    positions to represent promotor regions given user-defined 
107    upstream and downstream adjustments relative to the transcription
108    start site.
109
110    Parameters
111    ----------
112    gene_anno : pd.DataFrame 
113        A dataframe with columns `['chr', 'start', 'end', 'strand', 
114        'gene_name']` created from a GTF file.
115    
116    len_up : int
117        Number of base pairs upstream of the transcription start site 
118        the promotor region should be adjusted for.
119    
120    len_down : int
121        Number base pairs downstream of the transcription start site 
122        the promotor region should be adjusted for.
123    
124    Returns
125    -------
126    gene_anno : pd.DataFrame
127        A dataframe where `['start', 'end']` columns represent the 
128        start and end positions of inferred promotor regions for each 
129        annotation.
130    """
131    # Subsetting DataFrame to only required data
132    required_cols = ['chr', 'start', 'end', 'strand', 'gene_name']
133    gene_anno = gene_anno[required_cols]
134
135    # Adding a column that annotates transcription start site based on strand
136    tss = gene_anno[['start', 'end', 'strand']].apply(get_tss, axis = 1)
137    gene_anno.insert(1, column = 'tss', value = tss)
138
139    # Calculating start and end positions based on the tss
140    adj_regions = gene_anno[['tss', 'strand']].apply(lambda row:
141                                                     calc_range(
142                                                        row = row,
143                                                        len_up = len_up,
144                                                        len_down = len_down), 
145                                                   axis = 1, 
146                                                   result_type = 'expand'
147                                                   )
148    
149    gene_anno.loc[:, ['start', 'end']] = adj_regions.values
150    
151    return gene_anno
152
153
154def create_region_dicts(gene_anno : pd.DataFrame) -> dict:
155    """
156    Takes a GTF as a pd.DataFrame and returns dictionaries required for
157    region comparisons between gene_annotations and assay features.
158
159    Parameters
160    ----------
161    gene_anno : pd.DataFrame 
162        A dataframe with columns `['chr', 'start', 'end', 'strand', 
163        'gene_name']` created from a GTF file.
164    
165    Returns
166    -------
167    peak_gene_dict : dict
168        Keys are regions and values are genes.
169        
170    ga_regions : dict
171        Chromosomes are keys and regions are values.
172    """
173    peak_gene_dict = {}
174    ga_regions = {}
175
176    for _, anno in gene_anno.iterrows():
177
178        cur_chr = anno['chr']
179        cur_start = int(anno['start'])
180        cur_end = int(anno['end'])
181        cur_name = anno['gene_name']
182
183        peak_gene_dict[(cur_chr, cur_start, cur_end)] = cur_name
184
185        cur_region = np.array([[cur_start, cur_end]], dtype = int)
186        if cur_chr in ga_regions.keys():
187            ga_regions[cur_chr] = np.concatenate((ga_regions[cur_chr], 
188                                                  cur_region), axis = 0)
189        
190        else:
191            ga_regions[cur_chr] = cur_region
192
193    return peak_gene_dict, ga_regions
194
195
196def create_feature_dict(feature_names : list | set | np.ndarray | pd.Series
197                         ) -> dict:
198    """
199    Takes an array of feature names and returns data formatted for
200    feature comparisons.
201    
202    Returns
203    -------
204    feature_names : array_like
205        Array of region names from single-cell epigenetic data matrix.
206    
207    Returns
208    -------
209    feature_dict : dict
210        Keys are chromosomes and values are regions.
211    """
212    feature_dict = {}
213    feature_names = [re.split(":|-", peak) for peak in feature_names]
214    for peak_set in feature_names:
215        if peak_set[0] in feature_dict.keys():
216            feature_dict[peak_set[0]] = np.concatenate(
217                                                (feature_dict[peak_set[0]],
218                                                    np.array([[peak_set[1], 
219                                                                peak_set[2]]], 
220                                                                dtype = int)), 
221                                                                axis = 0)
222        else:
223            feature_dict[peak_set[0]] = np.array([[peak_set[1], 
224                                                   peak_set[2]]], 
225                                                   dtype = int)
226
227    return feature_dict
228
229
230def compare_regions(feature_dict : dict, ga_regions : dict,
231                     peak_gene_dict : dict, gene_sets : dict, chr_sep : str
232                     ) -> dict:
233    """
234    Takes features from a single-cell data matrix and regions from
235    a gene annotation file to return an region grouping where regions 
236    from feature_dict and regions from gene annotations overlap.
237
238    Parameters
239    ----------
240    feature_dict : dict
241        Keys are chromosomes and values are regions. This data should 
242        come from a single-cell experiment.
243
244    ga_regions : dict
245        Keys are chromosomes and values are regions. This data should 
246        come from a gene annotations (gtf) file.
247
248    peak_gene_dict : dict
249        Keys are peaks from gene annotation file and values are the 
250        gene they are associated with.
251
252    gene_sets : dict
253        Keys are gene set names and values are an iterable object of 
254        gene names.
255
256    chr_sep : str
257        The character that separates the chromosome from the rest of 
258        the region in the original feature array.
259    
260    Returns
261    -------
262    epi_grouping : dict
263        Keys are the names from gene_sets and values are a list of 
264        regions from `feature_names` that overlap with promotor regions 
265        respective to genes in `gene_sets` (i.e., if region in 
266        `feature_names` overlaps with promotor region from a gene in a 
267        gene set from `gene_sets`, that region will be added to the new 
268        dictionary under the respective gene set name).
269    """
270    epi_grouping = {group : [] for group in gene_sets.keys()}
271
272    for chrom in feature_dict.keys():
273        if chrom not in ga_regions.keys():
274            continue
275        for region in feature_dict[chrom]:
276            for anno in ga_regions[chrom]:
277                if _find_overlap(region[0], region[1], anno[0], anno[1]):
278                    gene = peak_gene_dict[(chrom, anno[0], anno[1])]
279                    for group in gene_sets.keys():
280                        if gene in gene_sets[group]:
281                            feat_region = ''.join((chrom, chr_sep, 
282                                                  str(region[0]), "-", 
283                                                  str(region[1])))
284                            epi_grouping[group].append(feat_region)
285
286    return epi_grouping
287
288
289def get_region_groupings(gene_anno : pd.DataFrame, gene_sets : dict, 
290                       feature_names : np.ndarray | pd.Series | list | set,
291                       len_up : int = 5000, len_down : int = 5000) -> dict:
292    """
293    Creates a peak set where keys are gene set names from `gene_sets` 
294    and values are arrays of features pulled from `feature_names`. 
295    Features are added to each peak set given overlap between regions 
296    in single-cell data matrix and inferred gene promoter regions in 
297    `gene_anno`.
298
299    Parameters
300    ----------
301    gene_anno : pd.DataFrame
302        Gene annotations in GTF format as a pd.DataFrame with columns
303        `['chr', 'start', 'end', 'strand', 'gene_name']`.
304
305    gene_sets : dict
306        Gene set names as keys and an iterable object of gene names
307        as values.
308
309    feature_names : array_like | set
310        Feature names corresponding to a single_cell epigenetic data 
311        matrix.
312
313    Returns
314    -------
315    epi_grouping : dict
316        Keys are the names from `gene_sets` and values
317        are a list of regions from `feature_names` that overlap with 
318        promotor regions respective to genes in `gene_sets` (i.e., if 
319        region in `feature_names` overlaps with promotor region from a 
320        gene in a gene set from `gene_sets`, that region will be added 
321        to the new dictionary under the respective gene set name).
322
323    Examples
324    --------
325    >>> # Reading in a gene set and the peak names from dataset
326    >>> gene_sets = np.load("data/RNA_hallmark_groupings.pkl", 
327    ...                     allow_pickle = True)
328    >>> peaks = np.load("data/MCF7_region_names.npy", 
329    ...                 allow_pickle = True)
330    >>> 
331    >>> # Reading in GTF file
332    >>> gtf_path = "data/hg38_subset_protein_coding.annotation.gtf"
333    >>> gtf = scmkl.read_gtf(gtf_path, filter_to_coding=True)
334    >>>
335    >>> region_grouping = scmkl.get_region_groupings(gene_anno = gtf,
336    ...                                         gene_sets = gene_sets,
337    ...                                         feature_names = peaks)
338    >>>
339    >>> region_grouping.keys()
340    dict_keys(['HALLMARK_TNFA_SIGNALING_VIA_NFKB', ...])
341    """
342    # Getting a unique set of gene names from gene_sets
343    all_genes = {gene for group in gene_sets.keys()
344                 for gene in gene_sets[group]}
345    
346    # Filtering out NaN values
347    all_genes = [gene for gene in all_genes if type(gene) != float]
348
349    # Filtering out annotations for genes not present in gene_sets
350    gene_anno = gene_anno[np.isin(gene_anno['gene_name'], all_genes)]
351
352    # Adjusting start and end columns to represent promotor regions
353    gene_anno = adjust_regions(gene_anno = gene_anno, 
354                                len_up = len_up, len_down = len_down)
355
356    # Creating a dictionary from assay features where [chr] : (start, end)
357    feature_dict = create_feature_dict(feature_names)
358
359    # Creating data structures from gene_anno for comparing regions
360    peak_gene_dict, ga_regions = create_region_dicts(gene_anno)
361
362    # Capturing the separator type used in assay
363    chr_sep = ':' if ':' in feature_names[0] else '-'
364
365    epi_groupings = compare_regions(feature_dict = feature_dict,
366                                     ga_regions = ga_regions,
367                                     peak_gene_dict = peak_gene_dict,
368                                     gene_sets = gene_sets,
369                                     chr_sep = chr_sep)
370    
371    return epi_groupings