Всестороннее введение в ваш геном с помощью стека SciPy

Биоинформатика — это междисциплинарная область, в которой разрабатываются методы и программные средства для анализа и понимания биологических данных.

Проще говоря, вы можете просто думать об этом как о науке о данных для биологии.

Среди многих типов биологических данных данные геномики являются одними из наиболее широко анализируемых. Особенно в связи с быстрым развитием технологий секвенирования ДНК следующего поколения (NGS) объем геномных данных растет в геометрической прогрессии. По словам Стивенса, Закари Д. и др., сбор геномных данных составляет эксабайт в год.

В этом посте я демонстрирую пример анализа файла GFF3 для генома человека с помощью стека SciPy. Generic Feature Format Version 3 (GFF3) — это текущий стандартный формат текстовых файлов для хранения геномных признаков. В частности, в этом посте вы узнаете, как использовать стек SciPy, чтобы ответить на следующие вопросы о геноме человека:

1. Какая часть генома неполная?
2. Сколько генов в геноме?
3. Какова длина типичного гена?
4. Как выглядит распределение генов по хромосомам?

Последний файл GFF3 для генома человека можно скачать здесь. Файл README, который находится в этом каталоге, содержит краткое описание этого формата данных, а более подробную спецификацию можно найти здесь.

Мы будем использовать Pandas, основной компонент стека SciPy, предоставляющий быстрые, гибкие и выразительные структуры данных, для обработки и понимания файла GFF3.

Настраивать

Прежде всего, давайте настроим виртуальную среду с установленным стеком SciPy. Этот процесс может занять много времени, если собирать его из исходного кода вручную, поскольку стек включает множество пакетов, некоторые из которых зависят от внешнего кода FORTRAN или C. Здесь я рекомендую использовать Miniconda, что делает настройку очень простой.

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

Флаг -b в строке bash указывает на выполнение в пакетном режиме. После того, как приведенные выше команды будут использованы для успешной установки Miniconda, запустите новую виртуальную среду для геномики, а затем установите стек SciPy.

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

Обратите внимание, что мы указали только 3 пакета, которые будем использовать в этом посте. Если вам нужны все пакеты, перечисленные в стеке SciPy, просто добавьте их в конец команды conda create .

Если вы не уверены в точном названии пакета, попробуйте conda search . Давайте активируем виртуальную среду и запустим IPython.

 source activate venv/ ipython

IPython — значительно более мощная замена интерфейса интерпретатора Python по умолчанию, поэтому все, что вы делали в интерпретаторе Python по умолчанию, также можно делать в IPython. Я настоятельно рекомендую каждому программисту Python, который еще не использовал IPython, попробовать его.

Загрузите файл аннотации

Теперь, когда наша настройка завершена, давайте загрузим файл аннотации генома человека в формате GFF3.

Это около 37 МБ, очень маленький файл по сравнению с информационным содержанием генома человека, который в обычном тексте составляет около 3 ГБ. Это связано с тем, что файл GFF3 содержит только аннотацию последовательностей, в то время как данные последовательности обычно хранятся в другом формате файла, называемом FASTA. Если вам интересно, вы можете скачать FASTA здесь, но мы не будем использовать данные последовательности в этом руководстве.

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

Префикс ! сообщает IPython, что это команда оболочки, а не команда Python. Однако IPython также может обрабатывать некоторые часто используемые команды оболочки, такие как ls , pwd , rm , mkdir , rmdir даже без префикса ! .

Взглянув на заголовок файла GFF, вы увидите множество строк метаданных/прагм/директив, начинающихся с ## или #! .

Согласно README, ## означает, что метаданные стабильны, а #! значит экспериментальный.

Позже вы также увидите ### , еще одну директиву с еще более тонким значением, основанным на спецификации.

Предполагается, что читаемые человеком комментарии должны быть после одного # . Для простоты мы будем рассматривать все строки, начинающиеся с # , как комментарии и просто игнорировать их при анализе.

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

Первая строка указывает, что в этом файле используется формат GFF версии 3.

Далее следуют сводки всех областей последовательности. Как мы увидим позже, такая информация также может быть найдена в основной части файла.

Строки, начинающиеся с #! показывает информацию о конкретной сборке генома, GRCh38.p7, к которой относится этот файл аннотации.

Genome Reference Consortium (GCR) — это международный консорциум, который наблюдает за обновлениями и улучшениями нескольких эталонных сборок генома, в том числе для человека, мыши, рыбки данио и курицы.

Просматривая этот файл, вот первые несколько строк аннотации.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

Столбцы seqid , источник , тип , начало , конец , счет , цепь , фаза , атрибуты . Некоторые из них очень легко понять. В качестве примера возьмем первую строку:

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

Это аннотация первой хромосомы с seqid 1, которая начинается с первого основания до 24 895 622-го основания.

Другими словами, длина первой хромосомы составляет около 25 миллионов оснований.

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

В последнем столбце атрибутов указано, что хромосома 1 также имеет три псевдонима, а именно CM000663.2, chr1 и NC_000001.11. Примерно так выглядит файл GFF3, но мы не будем проверять их построчно, поэтому пришло время загрузить весь файл в Pandas.

Pandas хорошо подходит для работы с форматом GFF3, потому что это файл с разделителями табуляцией, а Pandas очень хорошо поддерживает чтение CSV-подобных файлов.

Обратите внимание, что одним исключением из формата с разделителями табуляцией является случай, когда GFF3 содержит ##FASTA .

Согласно спецификации, ##FASTA указывает на конец части аннотации, за которой будет следовать одна или несколько последовательностей в формате FASTA (без разделителей табуляцией). Но это не относится к файлу GFF3, который мы собираемся анализировать.

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

Последняя строка выше загружает весь файл GFF3 с помощью метода pandas.read_csv .

Поскольку это не стандартный файл CSV, нам нужно немного настроить вызов.

Сначала мы сообщаем Pandas о недоступности информации о заголовке в GFF3 с помощью header=None , а затем указываем точное имя для каждого столбца с names=col_names .

Если аргумент name не указан, Pandas будет использовать возрастающие числа в качестве names для каждого столбца.

sep='\t' сообщает Pandas, что столбцы разделены табуляцией, а не запятой. Значение sep= на самом деле может быть регулярным выражением (регулярным выражением). Это удобно, если в файле используются разные разделители для каждого столбца (о да, такое бывает). comment='#' означает, что строки, начинающиеся с # , считаются комментариями и будут игнорироваться.

compression='gzip' сообщает Pandas, что входной файл сжат с помощью gzip.

Кроме того, pandas.read_csv имеет богатый набор параметров, которые позволяют читать различные форматы файлов, подобные CSV.

Тип возвращаемого значения — DataFrame , самая важная структура данных в Pandas, используемая для представления 2D-данных.

Pandas также имеет структуру данных Series и Panel для 1D и 3D данных соответственно. Пожалуйста, обратитесь к документации для ознакомления со структурами данных Pandas.

Давайте посмотрим на первые несколько записей с помощью метода .head .

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

Вывод хорошо отформатирован в табличном формате с более длинными строками в столбце атрибутов, частично замененными на ... .

Вы можете заставить Pandas не пропускать длинные строки с помощью pd.set_option('display.max_colwidth', -1) . Кроме того, в Pandas есть множество опций, которые можно настроить.

Далее давайте получим базовую информацию об этом кадре данных с помощью метода .info .

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

Это показывает, что GFF3 имеет 2 601 848 аннотированных строк, и каждая строка имеет девять столбцов.

Для каждого столбца он также показывает их типы данных.

Эти начало и конец имеют тип int64, целые числа, представляющие позиции в геноме.

Все остальные столбцы имеют тип object , что, вероятно, означает, что их значения состоят из смеси целых чисел, чисел с плавающей запятой и строк.

Размер всей информации составляет около 178,7+ МБ, хранящихся в памяти. Это получается более компактно, чем несжатый файл, который будет весить около 402 МБ. Быстрая проверка показана ниже.

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

С точки зрения высокого уровня мы загрузили весь файл GFF3 в объект DataFrame в Python, и весь наш последующий анализ будет основан на этом единственном объекте.

Теперь давайте посмотрим, что представляет собой первый столбец seqid.

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid — это один из способов доступа к данным столбца из фрейма данных. Другой способ — df['seqid'] , который является более общим синтаксисом, поскольку, если имя столбца является зарезервированным ключевым словом Python (например class ) или содержит расширение . или символ пробела, первый способ ( df.seqid ) не будет работать.

Выходные данные показывают, что существует 194 уникальных последовательности, которые включают хромосомы с 1 по 22, X, Y и митохондриальную (МТ) ДНК, а также 169 других последовательностей.

Секиды, начинающиеся с KI и GL, представляют собой последовательности ДНК — или каркасы — в геноме, которые не были успешно собраны в хромосомы.

Для тех, кто не знаком с геномикой, это важно.

Хотя первый проект генома человека был опубликован более 15 лет назад, текущий геном человека все еще неполный. Сложность сборки этих последовательностей во многом связана со сложными повторяющимися участками генома.

Далее, давайте взглянем на исходный столбец.

В README говорится, что источником является свободный текстовый квалификатор, предназначенный для описания алгоритма или рабочей процедуры, создавшей эту функцию.

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

Это пример использования метода value_counts , чрезвычайно полезного для быстрого подсчета категориальных переменных.

Из результата мы видим, что для этого столбца есть семь возможных значений, и большинство записей в файле GFF3 происходят из havana, ensembl и ensembl_havana.

Вы можете узнать больше о том, что означают эти источники и о взаимоотношениях между ними в этом посте.

Для простоты мы сосредоточимся на записях из источников GRCh38, havana, ensembl и ensembl_havan.a.

Какая часть генома неполная?

Информация о каждой целой хромосоме находится в записях из источника GRCh38, поэтому давайте сначала отфильтруем все остальное и присвоим отфильтрованный результат новой переменной gdf .

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

Фильтровать в Pandas легко.

Если вы проверите значение, полученное из выражения df.source == 'GRCh38' , это будет серия значений True и False для каждой записи с тем же индексом, что и df . Передача его в df[] вернет только те записи, для которых соответствующие значения равны True.

В df[] есть 194 ключа, для которых df.source == 'GRCh38' .

Как мы видели ранее, в столбце seqid также есть 194 уникальных значения, что означает, что каждая запись в gdf соответствует определенному seqid.

Затем мы случайным образом выбираем 10 записей с помощью sample метода, чтобы рассмотреть их поближе.

Вы можете видеть, что несобранные последовательности относятся к типу суперконтигов, а другие — к хромосомам. Чтобы вычислить долю неполного генома, нам сначала нужно знать длину всего генома, которая представляет собой сумму длин всех секидов.

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

В приведенном выше фрагменте мы сначала сделали копию gdf с помощью .copy() . В противном случае исходный gdf является просто фрагментом df , и его прямое изменение приведет к SettingWithCopyWarning (подробности см. здесь).

Затем мы вычисляем длину каждой записи и добавляем ее обратно в gdf в виде нового столбца с именем «длина». Общая длина оказывается около 3,1 миллиарда, а доля несобранных последовательностей — около 0,37%.

Вот как работает нарезка в последних двух командах.

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

Знак минус ( - ) добавляется к началу индекса, чтобы отменить выборку, потому что мы на самом деле хотим все, чего нет в списке (т.е. мы хотим несобранные, начинающиеся с KI и GL)…

Примечание. Поскольку собранные и несобранные последовательности различаются столбцом типа, последнюю строку можно альтернативно переписать следующим образом, чтобы получить одинаковые результаты.

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

Сколько существует генов?

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

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

Метод isin снова используется для фильтрации.

Затем быстрый подсчет значений показывает, что большинство записей представляют собой экзон, кодирующую последовательность (CDS) и нетранслируемую область (UTR).

Это субгенные элементы, но мы в основном ищем количество генов. Как показано, их 42 470, но мы хотим знать больше.

В частности, как их зовут и чем они занимаются? Чтобы ответить на эти вопросы, нам нужно внимательно изучить информацию в столбце атрибутов.

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

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

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

Во-первых, мы извлекаем имена генов.

В регулярном выражении Name=(?P<gene_name>.+?); , +? используется вместо + , потому что мы хотим, чтобы он не был жадным и чтобы поиск останавливался на первой точке с запятой; в противном случае результат будет соответствовать последней точке с запятой.

Кроме того, регулярное выражение сначала компилируется с помощью re.compile , а не используется напрямую, как в re.search для повышения производительности, потому что мы будем применять его к тысячам строк атрибутов.

extract_gene_name служит вспомогательной функцией для использования в pd.apply , который является методом, используемым, когда функция должна применяться к каждой записи фрейма данных или серии.

В этом конкретном случае мы хотим извлечь имя гена для каждой записи в ndf.attributes и добавить имена обратно в ndf в новый столбец с именем gene_name .

Идентификаторы генов и описание извлекаются аналогичным образом.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

Регулярное выражение для RE_GENE_ID немного более конкретное, так как мы знаем, что каждый gene_id должен начинаться с ENSG , где ENS означает ансамбль, а G означает ген.

Для записей, у которых нет описания, мы вернем пустую строку. После того, как все будет извлечено, мы больше не будем использовать столбец атрибутов, поэтому давайте удалим его, чтобы все было красиво и чисто с помощью метода .drop :

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

В приведенном выше вызове attributes указывают конкретный столбец, который мы хотим удалить.

axis=1 означает, что мы удаляем столбец вместо строки (по умолчанию axis=0 ).

inplace=True означает, что сброс выполняется на самом DataFrame вместо того, чтобы возвращать новую копию с указанным удаленным столбцом.

Быстрый взгляд на .head показывает, что столбец атрибутов действительно исчез, и были добавлены три новых столбца: gene_name , gene_id и desc .

Из любопытства посмотрим, уникальны ли все gene_id и gene_name :

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

Удивительно, но количество имен генов меньше, чем количество идентификаторов генов, что указывает на то, что некоторое имя_гена должно соответствовать нескольким идентификаторам генов. Давайте узнаем, что они из себя представляют.

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

Мы группируем все записи по значению gene_name , затем подсчитываем количество элементов в каждой группе с помощью .count() .

Если вы проверите вывод ndf.groupby('gene_name').count() , все столбцы подсчитываются для каждой группы, но большинство из них имеют одинаковые значения.

Обратите внимание, что значения NA не будут учитываться при подсчете, поэтому учитывайте только первый столбец, seqid (мы используем .ix[:, 0] , чтобы гарантировать отсутствие значений NA).

Затем отсортируйте значения счетчика с помощью .sort_values ​​и измените порядок с помощью .ix[::-1] .

В результате имя гена может быть общим для семи идентификаторов генов.

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

Более пристальный взгляд на все гены SCARNA20 показывает, что они действительно все разные.

Хотя они имеют одно и то же имя, они расположены в разных позициях генома.

Их описания, однако, не кажутся очень полезными для их различения.

Суть здесь в том, чтобы знать, что имена генов не уникальны для всех идентификаторов генов, и около 0,15% имен являются общими для нескольких генов.

Какова длина типичного гена?

Подобно тому, что мы делали, когда пытались понять неполноту генома, мы можем легко добавить столбец length в ndf :

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() вычисляет простую статистику на основе значений длины:

• Средняя длина гена составляет около 36 000 оснований.

• Средняя длина гена составляет около 5200 оснований.

• Минимальная и максимальная длина гена составляет около восьми и 2,3 миллиона оснований соответственно.

Поскольку среднее значение намного больше медианы, это означает, что распределение длин смещено вправо. Чтобы иметь более конкретный вид, давайте построим распределение.

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

Pandas предоставляет простой интерфейс для matplotlib, чтобы сделать графику очень удобной с DataFrames или сериями.

В этом случае он говорит, что нам нужен график гистограммы ( kind='hist' ) с 50 ячейками, и пусть ось Y будет в логарифмическом масштабе ( logy=True ).

Из гистограммы видно, что большинство генов находятся в первой ячейке. Однако длина некоторых генов может превышать два миллиона оснований. Давайте узнаем, что они из себя представляют:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

Как видите, самый длинный ген называется CNTNAP2, что является сокращением от контактин-ассоциированного белка-подобного 2. Согласно его странице в Википедии,

Этот ген охватывает почти 1,6% хромосомы 7 и является одним из крупнейших генов в геноме человека.

Конечно! Мы только что убедились в этом сами. Напротив, как насчет самых маленьких генов? Оказывается, они могут быть короче восьми оснований.

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

Длины двух крайних случаев различаются на пять порядков (2,3 миллиона против 8), что огромно и может указывать на уровень разнообразия жизни.

Один ген может быть транслирован во множество различных белков с помощью процесса, называемого альтернативным сплайсингом, который мы еще не исследовали. Такая информация также находится внутри файла GFF3, но выходит за рамки этого поста.

Распределение генов среди хромосом

Последнее, что я хотел бы обсудить, — это распределение генов по хромосомам, которое также служит примером для введения метода .merge для объединения двух DataFrames. Интуитивно понятно, что более длинные хромосомы, вероятно, содержат больше генов. Давайте посмотрим, правда ли это.

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

Мы позаимствовали переменную chrs из предыдущего раздела и использовали ее для фильтрации несобранных последовательностей. Судя по результату, самая большая хромосома 1 действительно содержит больше всего генов. Хотя хромосома Y имеет наименьшее количество генов, это не самая маленькая хромосома.

Обратите внимание, что в митохондриях (МТ) вроде бы нет генов, что не соответствует действительности.

Еще немного фильтрации первого DataFrame df , возвращенного pd.read_csv , показывает, что все гены MT взяты из источника insdc (которые были отфильтрованы ранее при создании edf , где мы рассматривали только источники havana, ensembl или ensembl_havana).

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

1. Около 0,37% генома человека до сих пор не завершены, несмотря на то, что первый набросок был опубликован более 15 лет назад.
2. Согласно этому конкретному файлу GFF3, который мы использовали, в геноме человека содержится около 42 000 генов.
3. Длина гена может варьироваться от нескольких десятков до более чем двух миллионов оснований.
4. Гены неравномерно распределены по хромосомам. В целом, чем больше хромосома, тем больше в ней генов, но для подмножества хромосом корреляция может быть отрицательной.

Файл GFF3 очень богат аннотационной информацией, и мы только поверхностно коснулись его. Если вы заинтересованы в дальнейшем исследовании, вот несколько вопросов, с которыми вы можете поиграть:

1. Сколько транскриптов обычно имеет ген? Какой процент генов имеет более 1 транскрипта?
2. Сколько изоформ обычно имеет транскрипт?
3. Сколько экзонов, CDS и UTR обычно имеет транскрипт? Каких размеров они?
4. Можно ли классифицировать гены на основе их функции, как описано в столбце описания?