Update 2023-Spring/2023-02-14-session-03/Analysis.ipynb

%% Cell type:markdown id:61dfe56d-4ab8-4107-97dc-6c3c4dad1c07 tags:

## Analysis or GT repeats

Example with repeats found on chromosome "NC_000024.10" for its 185 genes defined within a `json` file.

The repeats are located in plain text files `Repeats_Y_+.txt` and `Repeats_Y_-.txt` named after their `+` or `-` orientation.

The following files should be within the same directory for ease of use:

- Analysis.ipynb (this noteboo

- NC_000024.10_genes.json

- Repeats_Y_+.txt

- Repeats_Y_-.txt

The purpose of the code is to evaluate if a repeat is a *subsequence* of a gene, and then report the *begining* and *end* ranges for both the genes and the repeats.

This Python code therefore illustrates the use of `class` to define Python objects and they `methods` contained within that are the computations we want to run.

---

Thanks to Cristian for volunteering and demonstrating this real life example.

%% Cell type:code id:8523646c-39ed-40a3-93c7-92eaf03073a5 tags:

``` python

# -*- coding: utf-8 -*-

import json

def open_annotation(annotation_file):

    """Open anotation file.

    Parameters

    ----------

    annotation_file : string

        Annotation file name.

    Returns

    -------

    genes : List

        List containing genes.

    """

    # open file

    file = open(annotation_file, 'r')

    data = json.load(file)

    file.close()

    # List of genes

    genes = []

    for geneID in data:

        # get gene data

        gene_dictionary = data[geneID]

        seq_range = gene_dictionary['seq_range']

        seq_type = 'gene'

        orientation = gene_dictionary['orientation']

        gene = Gene(seq_type, seq_range, orientation, geneID)

        genes.append(gene)

    return genes

def open_repeats(repeat_file, orientation):

    """Open repeat data file.

    Parameters

    ----------

    repeat_file : String

        Repeats data file.

    orientation : String

        Repeat sequence orientation, it can be '+' or '-'.

    Returns

    -------

    repeats : List

        List of repeats sequences.

    """

    # Open data file

    file = open(repeat_file, 'r')

    lines = file.readlines()

    file.close()

    # repeats list

    repeats = []

    for index in range(1, len(lines)):

        line = lines[index].split()

        start = int(line[0])

        end = int(line[1])

        repeat = Sequence('repeat', (start, end), orientation)

        repeats.append(repeat)

    return repeats

```

%% Cell type:code id:ca39b0c3-cbcd-4c86-a277-840dc97e62d0 tags:

``` python

class Sequence():

    """ DNA sequence Class"""

    def __init__(self, seq_type, seq_range, orientation):

        self.seq_type = seq_type

        self.seq_range = seq_range

        self.orientation = orientation

    def __str__(self):

        info = '{0} range {1}'.format(self.seq_type, self.seq_range)

        return info

    # Define subsequence

    def is_subsequence(self, other_sequence):

        other_start, other_end = other_sequence.seq_range

        start, end = self.seq_range

        same_orientation = self.orientation == other_sequence.orientation

        return start >= other_start and end <= other_end and same_orientation

```

%% Cell type:code id:f290bdea-5a35-4f94-8787-7ceb354ec7d7 tags:

``` python

s = Sequence('Repeat', (125, 136), '+')

```

%% Cell type:code id:c8a0628e-10a6-41a4-a13e-50d3bb7d9fa1 tags:

``` python

s.seq_type

```

%% Output

    'Repeat'

%% Cell type:code id:08b538ea-8285-4f35-b559-81d8696ded3f tags:

``` python

print(s.seq_type)

```

%% Output

    Repeat

%% Cell type:code id:691b813e-c255-4324-98d4-3b53665e21f4 tags:

``` python

g = Sequence('Gene', (0,1000), '+')

s = Sequence('Repeat', (125, 1360), '+')

```

%% Cell type:code id:c5eea42f-be57-4dba-9f81-7a5a57fa00f1 tags:

``` python

print(s.is_subsequence(g))

```

%% Output

    False

%% Cell type:code id:0f545bcb-ccfb-45e0-95bb-2eee770d1835 tags:

``` python

# New class

class Gene(Sequence):

    def __init__(self, seq_type, seq_range, orientation, geneID):

        self.seq_type = seq_type

        self.seq_range =  seq_range

        self.orientation = orientation

        self.geneID = geneID

        self.repeats = []

    def add_repeat(self, repeat):

        self.repeats.append(repeat)

```

%% Cell type:code id:6c245bd3-e808-42a4-bc26-f895e4c65e1e tags:

``` python

# Function to compare each gene to all other repeats

def analysis(genes, repeats):

    for gene in genes:

        for repeat in repeats:

            result = repeat.is_subsequence(gene)

            if result:

                print(gene, repeat)

```

%% Cell type:code id:a026dd34-2417-4f87-ad3e-a3449ad4c9de tags:

``` python

# Make a list of repeats and a list of genes

r = open_repeats('Repeats_Y_+.txt', '+')

a = open_annotation('NC_000024.10_genes.json')

```

%% Cell type:code id:ec9562e6-9fd6-4087-bd0f-adbb7c5ed43f tags:

``` python

analysis(a,r)

```

%% Output

    gene range [2935381, 2982508] repeat range (2941936, 2941963)

    gene range [5000044, 5742228] repeat range (5064668, 5064701)

    gene range [5000044, 5742228] repeat range (5068544, 5068579)

    gene range [5000044, 5742228] repeat range (5085781, 5085812)

    gene range [5000044, 5742228] repeat range (5096439, 5096476)

    gene range [5000044, 5742228] repeat range (5113798, 5113825)

    gene range [5000044, 5742228] repeat range (5154627, 5154650)

    gene range [5000044, 5742228] repeat range (5156163, 5156189)

    gene range [5000044, 5742228] repeat range (5167519, 5167562)

    gene range [5000044, 5742228] repeat range (5318715, 5318746)

    gene range [5000044, 5742228] repeat range (5351296, 5351333)

    gene range [5000044, 5742228] repeat range (5360331, 5360362)

    gene range [5000044, 5742228] repeat range (5410039, 5410066)

    gene range [5000044, 5742228] repeat range (5525945, 5525968)

    gene range [5000044, 5742228] repeat range (5634788, 5634820)

    gene range [5000044, 5742228] repeat range (5732569, 5732592)

    gene range [6908594, 7107154] repeat range (6975903, 6975929)

    gene range [7458527, 7520566] repeat range (7468429, 7468452)

    gene range [9801153, 9813245] repeat range (9801487, 9801523)

    gene range [9801153, 9813245] repeat range (9803282, 9803325)

    gene range [12701231, 12860839] repeat range (12754207, 12754232)

    gene range [12701231, 12860839] repeat range (12767608, 12767639)

    gene range [12701231, 12860839] repeat range (12769778, 12769801)

    gene range [14522616, 14845654] repeat range (14526438, 14526461)

    gene range [14522616, 14845654] repeat range (14534769, 14534806)

    gene range [14522616, 14845654] repeat range (14540100, 14540124)

    gene range [14522616, 14845654] repeat range (14553568, 14553601)

    gene range [14522616, 14845654] repeat range (14637862, 14637884)

    gene range [14522616, 14845654] repeat range (14730989, 14731032)

    gene range [14522616, 14845654] repeat range (14792852, 14792891)

    gene range [20575776, 20593154] repeat range (20590112, 20590137)

    gene range [21511338, 21527212] repeat range (21513475, 21513520)

    gene range [21534879, 21559683] repeat range (21537016, 21537061)

    gene range [21534879, 21559683] repeat range (21551864, 21551901)

    gene range [22403410, 22419317] repeat range (22405572, 22405611)

    gene range [22438940, 22484714] repeat range (22441783, 22441819)

    gene range [22438940, 22484714] repeat range (22465311, 22465339)

%% Cell type:code id:4cd7a420-40b6-419d-af47-9cc08ca7d398 tags:

``` python

```