Скрипт №1 | Учебный сайт Екатерины Швецовой

Скрипт №1

Данный скрипт читает gbk-файл с описаниями и последовательностями генов, анализирует его и создает текстовый файл с таблицей, содержащей координаты и напревление генов. Скрипт также разбирает гены с неточными координатами и гены, содержащие несколько фрагментов.

Загрузить скрипт

Разбор скрипта

Для удобства вводим класс ProteinCodingGene. Его инициализирующая функция:


        def __init__ (self):

            self.start = ""

            self.end = ""

            self.direction = "1"

В данном скрипте я использую следующие методы класса ProteinCodingGene.
Метод, возвращающий строку с информацией о гене - объекте данного класса:


        def get_data (self):
      
            string = (self.start + "\t" + self.end + "\t"  + self.direction + "\n")

            return string

Метод, определяющий направление гена. Эта функция анализирует строку, содержащую координаты гена и, если та начинается с "complement", то присваивает гену направление -1 (по умолчанию все гены имеют направление 1, что видно из инициализирующей функции класса). Функция возвращает строку с координатами гена, но избавленную от "complement("


        def complement (self, string):

            x = string.split("(", 1)

            if x[0] == "complement":

                self.direction = "-1"

            x = string.replace("complement(","")

            return x

Данная функция анализирует строку с координатами гена и присваивает объекту класса (гену) координаты начала и конца. Она подходит и для генов, содержащих несколько фрагментов, т. к. за координату конца гена берется последний элемент списка, получившегося после разбиения строки по ".."


        def str_analysis (self, string):

            x = string.split("..")

            self.start = x[0]

            self.end = x[-1]

Метод, анализирующий гены с неточными координатами. Данная функция формирует список из символов строки, содержащей координаты гена и ищет среди них знаки больше или равно. Если таковые имееются, то они удаляются и к строке направления гена прибавляется 2 единицы. (предполагается, что мы уже определили изначальное направление гена, работая с ним как со стандартным геном, т. е. он уже имел direction 1 или -1. 2 единицы прибавляются, чтобы в будущем как-то отличить обычные гены от генов с неточными координатами.


        def inexact_coord (self, s):

            n = list(s)

            flag = 0

            for i in n:

                if i == ">" or i == "<":

                    self.start = self.start.replace (">","")

                    self.end = self.end.replace (">","")

                    self.start = self.start.replace("<","")

                    self.end = self.end.replace ("<","")

                    flag = 1

            if flag == 1:

                self.direction = (self.direction + "11")

Следующая часть скрипта получает данные об именах входного и выходного файлов (входной файл содержит информацию о генах, в выходной файл будет записана таблица с координатами и направлением генов) из коммандной строки и сохраняет в переменных filename1 и filename2.


        import sys

        if len(sys.argv) != 3:

            print "Script should receive 2 arguments:"

            print "First: name of the input file"

            print "Second: name of the output file"

            sys.exit()

        filename1 = sys.argv[1]
        filename2 = sys.argv[2]

Затем мы открываем первый файл для чтения, второй - для записи. Записываем шапку таблицы в выходной файл и создаем пустой список генов.


        myfile = open(filename1)

        newfile = open(filename2, "w")

        newfile.write ("start\tend\tdirection\n")

        gene_list = []

Начинаем читать файл построчно, пропуская пустые строки. Затем разбиваем строку по первому пробелу. Результатом, записанным в переменной step_str, будет список из двух строк.


        for step_str in myfile:
         
            if len(step_str.strip()) == 0:

                continue

            step_str = step_str.split (None,1)

Если нулевой элемент списка - "CDS" (или, что то же самое, если строка начинается с "CDS"), то мы создаем новую переменную gene, объект класса ProteinCodingGene. Далее работаем со строкой, содержащей координаты гена (она является первым элементом списка step_str): избавляемся в ней от незначащих символов, закрывающих скобок, а затем применяем к ней и к нашему гену метод complement. Т. е. на данном этапе определено направление гена, а в рабочей переменной w_str содержатся координаты гена, но уже без "complement(".


            if step_str[0] == "CDS":

                gene = ProteinCodingGene()

                w_str = step_str[1].strip()

                w_str = w_str.replace(")","")

                w_str = gene.complement(w_str)

Теперь возможно два варианта: либо ген стандартный, и тогда в w_str записано что-то вида "890786..896543", либо ген состоит из нескольких фрагментов и тогда строчка будет начинаться с "join(". Чтобы разделить эти случаи и сначала разбила строчку по знаку открывающейся скобки (если ген обычный, то скобок в строке вообще нет и список будет содержать один элемент). Если нулевой элемент списка - "join", то анализируем его первый элемент (который, что естественно, уже не будет содержать "join") с помощью метода str_analysis и приписываем к направлению гена единицу, чтобы отличать "join" гены от стандартных. Если же нулевой элемент списка - не "join", то ген стандартный, соответственно, анализируем нулевой элемент списка. На данном этапе gene уже имеет координаты и ориентацию.


                w_str = w_str.split("(")

                if w_str[0] == "join":

                    gene.str_analysis (w_str[1])

                    gene.direction = (gene.direction + "1")

                else:

                    gene.str_analysis (w_str[0])

Анализируем гены с неточными координатами. Для этого формируем строку, содержащую уже полученные нами координаты и направление гена (для этого используем метод get_data). Применяем к гену и данной строке метод inexact_coord. Наконец, добавляем полученный ген в ранее созданный список.


                gene_str = gene.get_data()

                gene.inexact_coord(gene_str)

                gene_list.append(gene)

Проходимся по полученному списку генов и записываем информацию о его координатах и направлении в выходной файл. Закрываем файл.


        for i in gene_list:

            newfile.write (i.get_data())

        newfile.close()