Генетические визуализатор последовательность генерирования изображения большого размера


Вдохновленный этой реддите пост, я написал программу, которая принимает последовательность ДНК (в формате fasta формат файла) и генерирует графическое представление. HIV virus Ebola virus

Вирус ВИЧ и Эболы

Моя программа работает так:

  1. Принять и обработать входной файл
  2. Разбить входной файл в массив (с помощью пользовательского класса array largearray.Array())
  3. Перебрать и разобрать, что построить путь, который хранится в другом largearray.Array что содержит координаты и цвет каждого пикселя на пути
  4. Создать образ и нарисовать путь на него
  5. Сохранить Изображение

Я разместил здесь ранее, прошу помочь с оптимизацией шаги 1-4. Graipher помог мне с его оптимизированной версии моего класса Array. Я теперь в состоянии обрабатывать очень большие входные файлы, но я ударил еще одним узким местом в отношении того, насколько большой изображение пильный можете создать. (Шаг 4)

Входной файл, который больше, чем примерно 450 МБ, создает образы, которые, видимо, слишком большой для обработки. Пильная работает в MemoryError в PIL.Image.new() функция. Что не удивительно, учитывая, что программа пытается создать изображения размером более 100 000 пикселей с большими файлами.

Traceback (most recent call last):
  File ".\genome-imager.py", line 113, in <module>
    with Image.new("RGBA", dim, None) as img:
  File "C:\Program Files\Python3\lib\site-packages\PIL\Image.py", line 2291, in new
    return Image()._new(core.new(mode, size))
MemoryError

Я потратил много времени на обдумывание решения, и я придумал это разделение изображения на блоки/квадраты, собирая по пути на каждый кусок по отдельности и затем сшивая их вместе. Я провел некоторое время, пытаясь сделать это, но я не могу найти то, что работает и мой недостаток опыта мешает.

Если у кого есть какие-либо советы или решения о том, как я могу преодолеть это узкое место я высоко ценю его. Моя цель-преодолеть 400 МБ барьер и перейти к геномы больше, чем ГБ. Возможно, геном человека в будущем.

Использования (работать с Python х64!):

python3 genome-imager.py path/to/file

Несколько примеров входных файлов и их результаты можно найти здесь. Читать contents.txt для получения дополнительной информации и который геномов нарваться на ошибку.

Важно:

Эта программа создает очень большие временные файлы (~100 МБ каждый) и может сбрасывать до 15 ГБ (и более) за раз. Файлы изображений также оказаться довольно большой (~100-200 Мб). Убедитесь, что у вас достаточно места. Эти временные файлы очищаются в конце программы, если он прогрессирует к завершению, но некоторые из них не удаляются, если программа останавливается на полпути. Помните, чтобы удалить их!

Эта программа также достаточно интенсивно

Код:

genome-imager.py

#!/usr/bin/env python3
# Python3
#
#   Run from command line
#

import logging
from argparse import ArgumentParser
from copy import deepcopy, copy
from datetime import timedelta
from math import ceil
from os import remove, makedirs
from os.path import exists
from re import sub
from time import time

from PIL import Image, ImageDraw

from largearray import Array

uuid = id(time())

parser = ArgumentParser()
parser.add_argument("file", help="Location of input file. path/to/file (FASTA file formats supported)")
parser.add_argument("-i", "--image-name",
                    help="Where to save finished image file. path/to/file (Default: Name_of_input_file.png)")
parser.add_argument("-s", "--dump-size", help="The size of temp files to dump to disk. (Default & Max: 5)", type=int)
parser.add_argument("-t", "--temp", help="Where to dump temp files. path/to/directory/ (Default: .cache/)", type=str)
parser.add_argument("-d", "--debug-file", help="Where to store debug file. path/to/file (Default ./cache/debug.log")
args = parser.parse_args()
filepath = args.file
ipath = ".".join(filepath.split(".")[:-1]) + ".png"
if args.image_name:
    ipath = args.image_name
print(ipath)
dsize = 5
if args.dump_size:
    dsize = args.dump_size
cachedir = ".cache/"
if args.temp:
    cachedir = args.temp
debugpath = '.cache/debug%d.log' % uuid
if args.debug_file:
    debugpath = args.debug_file
if not exists(filepath):
    raise Exception("Path of input file does not exist")
print("Debug at %s" % debugpath)
if exists(debugpath):
    remove(debugpath)
if not exists(cachedir):
    makedirs(cachedir)
logging.basicConfig(filename=debugpath, level=logging.DEBUG)
logging.info("Init: %d" % uuid)

del parser, ArgumentParser, remove, exists,

print("Generating vizualization of %s" % filepath)
starttime = time()
file = open(filepath, 'r')
logging.info("File opened")
logging.info("Serializing %s ..." % filepath)
raw = ''.join([n for n in file.readlines() if not n.startswith('>')]).replace('\n', "").lower()
logging.info("Replaced FASTA info")
file.close()
del file
raw = sub(r'[rykmswbdhv-]', "n", raw)  # Handles miscellaneous FASTA characters
raw = sub(r'[^atgcn]', "", raw)  # Handles 4 bases and not-known

sequence = Array(name="sequence", cachedirectory=cachedir, a=list(raw), maxitems=(dsize * 10))
sequence.trim()
logging.info("Parsed characters (%d items)" % len(sequence))
del sub, raw
endtime = [time()]
print("The input file has been serialized. %s (%d items) Calculating path..." % (
    str(timedelta(seconds=(endtime[0] - starttime))), len(sequence)))

action = {  # The different bases and their respective colours and movements
    "a": ((0, 255, 0), 0, -1),  # green - Moves up
    "t": ((255, 0, 0), 0, 1),  # red - Moves Down
    "g": ((255, 0, 255), -1, 0),  # hot pink - Moves Left
    "c": ((0, 0, 255), 1, 0),  # blue - Moves Right
    "n": ((0, 0, 0), 0, 0),  # black - Stays on spot
}

maxp = [[0, 0], [0, 0]]  # Top left and bottom right corners of completed path
curr = [0, 0]

pendingactions = Array(name="pendingactions", cachedirectory=cachedir, maxitems=dsize)
logging.info("%d temp files will be created [pendingactions]" % ceil(len(sequence) / pendingactions.maxitems))

for i in sequence:
    # get the actions associated from dict
    curr[0] += action[i][1]
    curr[1] += action[i][2]
    if curr[0] > maxp[0][0]:
        maxp[0][0] = curr[0]
    elif curr[0] < maxp[1][0]:
        maxp[1][0] = curr[0]
    if curr[1] > maxp[0][1]:
        maxp[0][1] = curr[1]
    elif curr[1] < maxp[1][1]:
        maxp[1][1] = curr[1]
    pendingactions.append((copy(curr), action[i][0]))
pendingactions.trim()
del sequence.a
del sequence, copy, deepcopy

# Final dimensions of image + 10px border
dim = (abs(maxp[0][0] - maxp[1][0]) + 20, abs(maxp[0][1] - maxp[1][1]) + 20)
endtime.append(time())
print("The path has been calculated. %s Rendering image... %s" % (
    str(timedelta(seconds=(endtime[1] - starttime))), "(" + str(dim[0]) + "x" + str(dim[1]) + ")"))

with Image.new("RGBA", dim, None) as img:
    logging.info("Initial image created. (%d x %d)" % (dim[0], dim[1]))
    draw = ImageDraw.Draw(img)
    logging.info("Draw object created")

    for i in pendingactions:
        draw.point((i[0][0] + abs(maxp[1][0]) + 10, i[0][1] + abs(maxp[1][1]) + 10), fill=i[1])
    logging.info("Path Drawn")


    def mean(n):  # I couldn't find an average function in base python
        s = float(n[0] + n[1]) / 2
        return s


    # Start and End points are dynamically sized to the dimensions of the final image
    draw.ellipse([((abs(maxp[1][0]) + 10) - ceil(mean(dim) / 500), (abs(maxp[1][1]) + 10) - ceil(mean(dim) / 500)),
                  ((abs(maxp[1][0]) + 10) + ceil(mean(dim) / 500), (abs(maxp[1][1]) + 10) + ceil(mean(dim) / 500))],
                 fill=(255, 255, 0), outline=(255, 255, 0))  # yellow
    draw.ellipse([((curr[0] + abs(maxp[1][0]) + 10) - ceil(mean(dim) / 500),
                   (curr[1] + abs(maxp[1][1]) + 10) - ceil(mean(dim) / 500)), (
                      (curr[0] + abs(maxp[1][0]) + 10) + ceil(mean(dim) / 500),
                      (curr[1] + abs(maxp[1][1]) + 10) + ceil(mean(dim) / 500))], fill=(51, 255, 255),
                 outline=(51, 255, 255))  # neon blue
    logging.info("Start and End points drawn")

    del pendingactions.a
    del maxp, curr, mean, dim, draw, ImageDraw, pendingactions

    endtime.append(time())
    print("The image has been rendered. %s Saving..." % str(timedelta(seconds=(endtime[2] - endtime[1]))))
    img.save(ipath, "PNG", optimize=True)
    logging.info("Image saved at %s" % ipath)

endtime.append(time())
del img, Image
print("Done! %s Image is saved as: %s" % (str(timedelta(seconds=(endtime[3] - endtime[2]))), ipath))
print("Program took %s to run" % str(timedelta(seconds=(endtime[3] - starttime))))
logging.info("%s | %s | %s | %s # Parsing File | Computing Path | Rendering | Saving" % (
    str(timedelta(seconds=(endtime[0] - starttime))), str(timedelta(seconds=(endtime[1] - starttime))),
    str(timedelta(seconds=(endtime[2] - starttime))), str(timedelta(seconds=(endtime[3] - starttime)))))

largearray.py

#!/usr/bin/env python3
# Python3
#
#   Simple array class that dynamically saves temp files to disk to conserve memory
#

import logging
import pickle
from datetime import timedelta
from itertools import islice
from os import makedirs, remove
from os.path import exists
from shutil import rmtree
from time import time

startime = time()
logging.getLogger(__name__).addHandler(logging.NullHandler())
class Array():
    """1D Array class
    Dynamically saves temp files to disk to conserve memory"""

    def __init__(self, name="Array", cachedirectory=".cache/", a=None, maxitems=1):
        # How much data to keep in memory before dumping to disk
        self.maxitems = int(maxitems*1e6)
        self.fc = 0  # file counter
        self.uuid = id(self)
        self.name = name
        logging.debug("[largearray.Array] Instance %d %s created | %s" % (self.uuid, self.name, str(timedelta(seconds=time()-startime))))
        self.dir = cachedirectory + str(self.uuid) # make a unique subfolder (unique as long as the array exists)
        if exists(self.dir):
            rmtree(self.dir)
        makedirs(self.dir)
        logging.debug("[largearray.Array] Instance %d caches in %s with %d items per file" % (self.uuid, self.dir, self.maxitems))
        self.path = self.dir + "/temp%d.dat"  # Name of temp files
        self.hastrim = False
        self.a = []
        if a is not None:
            self.extend(a)

    def append(self, n):
        """Append n to the array.
        If size exceeds self.maxitems, dump to disk
        """
        if self.hastrim:
            raise Exception("ERROR: Class [array] methods append() and extend() cannot be called after method trim()")
        else:
            self.a.append(n)
            if len(self.a) >= self.maxitems:
                logging.debug("[largearray.Array] Instance %d dumps temp %d | %s" % (self.uuid, self.fc, str(timedelta(seconds=time()-startime))))
                with open(self.path % self.fc, 'wb') as pfile:
                    pickle.dump(self.a, pfile)  # Dump the data
                self.a = []
                self.fc += 1

    def trim(self):
        """If there are remaining values in the array stored in memory, dump them to disk (even if there is less than maxitems.
        NOTE: Only run this after all possible appends and extends have been done
        WARNING: This cannot be called more than once, and if this has been called, append() and extend() cannot be called again"""
        if len(self.a) > 0:
            if self.hastrim:
                raise Exception("ERROR: Class [array] method trim() can only be called once")
            else:
                self.hastrim = True
                self.trimlen = len(self.a)
                logging.debug("[largearray.Array] Instance %d trims temp %d | %s" % (self.uuid, self.fc, str(timedelta(seconds=time()-startime))))
                with open(self.path % self.fc, 'wb') as pfile:
                    pickle.dump(self.a, pfile)  # Dump the data
                self.a = []
                self.fc += 1

    def extend(self, values):
        """Convenience method to append multiple values"""
        for n in values:
            self.append(n)

    def __iter__(self):
        """Allows iterating over the values in the array.
        Loads the values from disk as necessary."""
        for fc in range(self.fc):
            logging.debug("[largearray.Array] Instance %d iterates temp %d | %s" % (self.uuid, fc, str(timedelta(seconds=time()-startime))))
            with open(self.path % fc, 'rb') as pfile:
                yield from pickle.load(pfile)
        yield from self.a

    def __repr__(self):
        """The values currently in memory"""
        s = "[..., " if self.fc else "["
        return s + ", ".join(map(str, self.a)) + "]"

    def __getitem__(self, index):
        """Get the item at index or the items in slice.
        Loads all dumps from disk until start of slice for the latter."""
        if isinstance(index, slice):
            return list(islice(self, index.start, index.stop, index.step))
        else:
            fc, i = divmod(index, self.maxitems)
            with open(self.path % fc, 'rb') as pfile:
                return pickle.load(pfile)[i]

    def __len__(self):
        """Length of the array (including values on disk)"""
        if self.hastrim:
            return (self.fc-1) * self.maxitems + self.trimlen
        return self.fc * self.maxitems + len(self.a)

    def __delattr__(self, item):
        """Calling del <object name>.a
        will delete entire array"""
        if item == 'a':
            super().__delattr__('a')
            rmtree(self.dir)
            logging.debug("[largearray.Array] Instance %d deletes all array data | %s" % (self.uuid, str(timedelta(seconds=time()-startime))))
        else:
            super(Array, self).__delattr__(item)

    def __setitem__(self, key, value):
        if isinstance(key, slice):
            l = list(islice(self, key.start, key.stop, key.step))
            for i in l:
                l[i].__setitem__(value)
                set()
        else:
            fc, i = divmod(key, self.maxitems)
            with open(self.path % fc, 'rb') as pfile:
                l = pickle.load(pfile)
                l[i] = value
            remove(self.path % fc)
            with open(self.path % fc, 'wb') as pfile:
                pickle.dump(l, pfile)

    def __delitem__(self, key):
        fc, i = divmod(key, self.maxitems)
        with open(self.path % fc, 'rb') as pfile:
            l = pickle.load(pfile)
            del l[i]
        remove(self.path % fc)
        with open(self.path % fc, 'wb') as pfile:
            pickle.dump(l, pfile)


Комментарии
2 ответа

Парсер

Ваш парсер есть ошибка в строке 62:

raw = ''.join([n for n in file.readlines() if not n.startswith('>')]).replace('\n', "").lower()

будет объединять все последовательности, несмотря на их потенциально принадлежащих к различным последовательностям ДНК. Вот Википедия (лучший источник для цитирования хD):


Множественный формат fasta последовательности может быть получена путем объединения нескольких отдельных файлов в формате fasta последовательности. Это не означает противоречие [...] следовательно, заставляя всех последующих последовательностей начать с ">" [...].

Вы же не уважают Заголовок (;) либо. Далее


Последовательности могут быть последовательностями белка или нуклеиновой кислоты последовательности

Теперь моя биология не хватает, в лучшем случае. Я предполагаю, что вы не уважаете, потому что вы можете сделать предположение только когда-либо иметь дело с нуклеиновыми кислотами в ДНК? Поэтому вы игнорируете проверку на что?

Вы, кажется, дальше, похоже, удалось сократить это только A C G N T что я думаю, также имеет очень хороший биологическая причина того, что я не знаю.

Это, вероятно, взорвется, если вы не осторожны с тем, что файлы, которые вы помещаете в.


Большой Массив

Есть несколько оптимизаций можно сделать, но только 1, что действительно имеет смысл здесь: Бен она. Вам не нужно это.

Вы предварительно вычислить много информации о том, что вы повторно использовать ровно один раз и которое можно вычислить на лету. Нет никакой выгоды в этом. На самом деле, вы делаете это хуже, прочитав .fna раз и потом читать такие версии, что на 40+ раз больше. Не круто :D

Короткая версия заключается в том, что largearray это упаковочная питон list()С магазина 5-tupels (x,y,r,g,b) стоит 40 байт на запись. strawberry.fna берет ~8.5 GB для представления в этом формате, в то время как вы можете, вместо того чтобы придерживаться одного персонажа стоит 1 байт и посмотри вверх (r,g,b) значение, а также вычислить (x,y) используя predecesor (которых вы знаете). Это позволит сэкономить 40(!) раз диск/памяти (до200 MB после оптимизации).

Вот версия, которая использует largearrayно спасает только персонаж (стоит 1 байт) вместо tupel (40 байт): https://gist.github.com/FirefoxMetzger/9cfdf91f9b838f755d8ae68ab2a52ac8

Я также добавил numpy как зависимость и рефакторинг кода (много), потому что я не мог сделать головы и хвосты от него.

Вот еще более короткий вариант, что никто больше не нужен largearray:
https://gist.github.com/FirefoxMetzger/95b4f30fee91a21ab376cd4db9ccce42

Я также снял ГРМ замечания, потому что это может быть сделано более эффективно cProfile или другой профайлер на ваш выбор. В этом случае, я чувствовал, что отпечатки мешали много кода, теперь легче читать. Если мигающий курсор без выходных вас раздражает, не стесняйтесь добавить больше отпечатков (в ущерб удобочитаемости и чуть-чуть во время выполнения).


Выход

Я побежал strawberry.fna и он пытался создать 94980x100286 изображения. Это ~9.5 Gpx (гигапикселя)! 29.5 GB как несжатый образ не весело. По крайней мере, на моей машине он не вписывается в несжатом памяти. Аналогично, если я смотрю на свой strawberry.png не похоже, верно; я только красноватые тени.

На вашей машине (если у вас есть больше одного, так как ты, видимо может поместиться 30+ГБ в памяти без проблем) вы могли бы получить больше от этого.

К счастью, большинство из этих пикселей пустой / прозрачный. На самом деле для strawberry.fna только 2,2% пикселей не прозрачный. Возможно, даже меньше, потому что я заметил, что они иногда пересекаются, например, в "ГК" последовательность будет виден как "C" торчал в сторону.

Моей попыткой будет использование библиотек matplotlib, однако это все еще довольно медленно и, как известно, не слишком хорошо справляется с большими данными. Она работает немного быстрее, чем ваш первоначальный вариант и требования к памяти лучше; все равно его не хватает strawberries.fna на моей машине. Однако, вы получаете интерактивный рисунок, что вы можете увеличить в и "скриншоты". Мне нравится этот способ лучше изучить последовательность.

Еще одним преимуществом является то, что можно указывать формат файла и просмотра, легенды, оси и так далее.

Вот код для этого:

#!/usr/bin/env python3
# Python3

import logging
from argparse import ArgumentParser
from os import remove
from os.path import exists
from re import sub
from time import time

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

uuid = id(time())

parser = ArgumentParser()
parser.add_argument("file", help="Location of input file. path/to/file (FASTA file formats supported)")
parser.add_argument("-d", "--debug-file", help="Where to store debug file. path/to/file (Default ./cache/debug.log")
args = parser.parse_args()
filepath = args.file
if not exists(filepath):
raise Exception("Path of input file does not exist")

debugpath = '.cache/debug%d.log' % uuid
if args.debug_file:
debugpath = args.debug_file
if exists(debugpath):
remove(debugpath)

logging.basicConfig(filename=debugpath, level=logging.DEBUG)
logging.info("Init: %d" % uuid)

with open(filepath, 'r') as file:
logging.info("Serializing %s ..." % filepath)
raw = ''.join([n for n in file.readlines() if not n.startswith('>')]).replace('\n', "").lower()

raw = sub(r'[rykmswbdhv-]', "n", raw) # Handles miscellaneous FASTA characters
raw = sub(r'[^atgcn]', "", raw) # Remove everything except 4 bases and not-known
test_length = min(5e6, len(raw))
print("Processing %d elements. Have fun waiting! I'll go grab a coffee"
% test_length)

action_lookup = {
"a": np.array((0, -1)),
"t": np.array((0, 1)),
"g": np.array((-1, 0)),
"c": np.array((1, 0)),
"n": np.array((0, 0))
}
color_lookup = {
"a": (0, 1, 0),
"t": (1, 0, 0),
"g": (1, 0, 1),
"c": (0, 0, 1),
"n": (0, 0, 0)
}

fig, ax = plt.subplots()

lines = np.zeros((int(test_length), 2, 2), dtype=np.int16)
colors = list()
cursor = np.array((0, 0))
old_cursor = cursor.copy()
for idx, base in enumerate(raw):
if idx >= test_length:
break
# get the actions associated from dict
cursor += action_lookup[base]
lines[idx, 0, :] = old_cursor
lines[idx, 1, :] = cursor
colors.append(color_lookup[base])
old_cursor = cursor.copy()

print("Done gathering data - plotting now")

linesegments = LineCollection(
lines,
colors=colors
)
ax.add_collection(linesegments)
plt.plot()
plt.show()

И вот некоторые сексуальные образы. (Они всегда начинаются в (0,0)).
ВИЧ.Фна:
hiv.fna
random2.Фна:
random2.fna
Демонстрируя урожай/зум-способность на random2:
Crop of Random2.fna

Я тоже могу ограничить количество обрабатываемых элементов. Таким образом, я могу визуализировать первые 10 миллионов элементов strawberry.fna. Я нахожу эту визуализацию, чтобы быть ограниченное применение, но думаю, что это имеет больше смысла для вас?

10M Strawberry.fna
Зум в о (0,0):
Strawberry zoom
И еще более полезный / бесполезный зум на (0,0):
Even more Zoom

6
ответ дан 14 апреля 2018 в 06:04 Источник Поделиться

Более красивые петли

В for i in sequence петли, вы могли бы избавиться от многих обращается с []:


  • color, dx, dy = action[i] получить доступ к action[i] только один раз и использовать кортеж распаковка для предотвращения доступа к элементам 1 и 2

  • использовать cur_x и cur_y вместо curr (это также устраняет необходимость для copy)

  • использовать max/min строение

  • определение 4 переменные вместо maxp

Вы получите что-то вроде:

max_x, max_y, min_x, min_y = 0, 0, 0, 0
cur_x, cur_y = [0, 0]
pendingactions = []

for i in sequence:
color, dx, dy = action[i]
cur_x += dx
cur_y += dy
max_x = max(cur_x, max_x)
min_x = min(cur_x, min_x)
max_y = max(cur_x, max_y)
min_y = min(cur_x, min_y)
pendingactions.append(([cur_x, cur_y], color))

Кроме того, вы можете получить мин/макс значений позвонив по телефону min/max на заключительном массива.

Такая же логика применяется во многих местах в коде.

Раз меры

Вместо того:

starttime = time()
...
endtime = [time()]
...str(timedelta(seconds=(endtime[0] - starttime)))
...
endtime.append(time())
...str(timedelta(seconds=(endtime[3] - endtime[2]))
...str(timedelta(seconds=(endtime[3] - starttime)))

Который должен индексы должны быть обновлены, если вы добавите частичное время, вы можете:


  • сохранить все меры времени, включая первый в массиве

  • полагаться на [-1] чтобы получить доступ к последнему элементу.

Фальшивый портной выше становится (/!\непроверено)

times = [time()]
...
times.append(time())
...str(timedelta(seconds=(times[-1] - times[0])))
...
times.append(time())
...str(timedelta(seconds=(times[-1] - times[-2])))
...str(timedelta(seconds=(times[-1] - times[0])))

Используйте инструменты для перевода из Python вместо регулярного выражения

Вместо регулярных выражений, может быть, вы могли бы использовать [maketrans] для обработки характер преобразований. Это может быть быстрее, может быть медленнее, я не выполнял каких-либо ориентира...

Магия чисел

У вас есть много магических чисел в коде: через "константы" для это поможет сделать код более удобным для чтения, имея значимое имя, а также легче поддерживать путем изменения значения в одном месте, если это необходимо.

3
ответ дан 14 апреля 2018 в 04:04 Источник Поделиться