суббота, 19 сентября 2015 г.

Загрузка файла формата GFF в рабочее пространство R среды.

Использование базовой R функции read.table.
Существует несколько путей загрузки данных, хранящихся в файлах формата GFF, в рабочее пространство R среды. Первый путь, самый простой, но не самый подходящий – с помощью стандартной функции read.table. Для корректной загрузки содержимого файлов формата GFF следует обратить внимание лишь на аргумент sep этой функции и с его помощью задать табуляцию как способ разделения столбцов таблицы. Значения остальных аргументов могут быть выбраны по умолчанию (в том числе аргумент comment.char, так как по умолчанию он имеет значение "#"). Ниже представлен код для загрузки, а также фрагмент итоговой таблицы и ее структура. В качестве примера взят файл формата GFF3, содержащий данные 81-го релиза геномного обозревателя Ensembl (GRCh38.p3 сборка генома человека, июнь 2015 года).

rm(list = ls())
setwd("D:/Transcriptome")
DataGFF = read.table(file = "Ensembl genes. GRCh38.p3, release 81, June 2015.gff3", sep = "\t")
print(DataGFF[1:5, ])
V1
V2
V3
V4
V5
V6
V7
V8

1
1
GRCh38
chromosome
1
248956422
.
.
.

2
1
havana
gene
11869
14409
.
+
.

3
1
havana
processed_transcript
11869
14409
.
+
.

4
1
havana
exon
11869
12227
.
+
.

5
1
havana
exon
12613
12721
.
+
.

V9
1
ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11
2
ID=gene:ENSG00000223972;Name=DDX11L1;biotype=transcribed_unprocessed_pseudogene;
description=DEAD/H (Asp-Glu-Ala-Asp/His) box helicase 11 like 1 [Source:HGNC Symbol%3BAcc:
HGNC:37102];gene_id=ENSG00000223972;havana_gene=OTTHUMG00000000961;havana_version=2;
logic_name=havana;version=5
3
ID=transcript:ENST00000456328;Parent=gene:ENSG00000223972;Name=DDX11L1-002;
biotype=processed_transcript;havana_transcript=OTTHUMT00000362751;havana_version=1;tag=basic;
transcript_id=ENST00000456328;transcript_support_level=1;version=2
4
Parent=transcript:ENST00000456328;Name=ENSE00002234944;constitutive=0;ensembl_end_phase=-1;
ensembl_phase=-1;exon_id=ENSE00002234944;rank=1;version=1
5
Parent=transcript:ENST00000456328;Name=ENSE00003582793;constitutive=0;ensembl_end_phase=-1;
ensembl_phase=-1;exon_id=ENSE00003582793;rank=2;version=1
str(DataGFF)
'data.frame': 1659566 obs. of 9 variables:
 $ V1: Factor w/ 492 levels "1","10","11",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ V2: Factor w/ 5 levels "ensembl","ensembl_havana",..: 3 4 4 4 4 4 4 4 4 4 ...
 $ V3: Factor w/ 31 levels "aberrant_processed_transcript",..: 4 7 17 5 5 5 19 5 5 5 ...
 $ V4: int  1 11869 11869 11869 12613 13221 12010 12010 12179 12613 ...
 $ V5: int  248956422 14409 14409 12227 12721 14409 13670 12057 12227 12697 ...
 $ V6: Factor w/ 1 level ".": 1 1 1 1 1 1 1 1 1 1 ...
 $ V7: Factor w/ 3 levels "-",".","+": 2 3 3 3 3 3 3 3 3 3 ...
 $ V8: Factor w/ 4 levels ".","0","1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ V9: Factor w/ 1117867 levels
 "ID=CDS:ENSP00000000442;Parent=transcript:ENST00000000442;protein_id=ENSP00000000442",..: 64816
 81414 153621 617878 617880 617879 150945 602003 601998 602000 ...

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

Использование функции read.gff.
Функция read.gff специализирована на загрузке файлов формата GFF с ftp-сервера NCBI. Эта функция доступна в виде исходного кода со страницы https://github.com/cstubben/genomes/blob/master/R/read.gff.R  сайта GitHub, либо как часть пакета genomes. Пакет genomes может быть установлен с официального сайта проекта Bioconductor с помощью стандартной процедуры:

source("https://bioconductor.org/biocLite.R")
biocLite("genomes")

Пример использования функции read.gff для загрузки файла формата GFF с ftp-сервера NCBI, а также фрагмент итогового объекта показаны ниже.

rm(list = ls())
setwd("D:/Transcriptome")
library(genomes)
sourceFile = "ftp://ftp.ncbi.nih.gov/genomes/Bacteria/Yersinia_pestis_CO92_uid57621/NC_003132.gff"
DataGFF = read.gff(file = sourceFile, locus.tags = TRUE, nrows = -1)
print(DataGFF)
GRanges object with 9 ranges and 5 metadata columns:
seqnames
ranges
strand
|
locus
pid
feature
gene
<Rle>
<IRanges>
<Rle>
|
<character>
<character>
<character>
<character>
[1]
NC_003132
[    87,
1109]
+
|
YPPCP1.01
NP_395227
CDS
[2]
NC_003132
[1106,
1888]
+
|
YPPCP1.02
NP_395228
CDS
[3]
NC_003132
[2925,
3119]
+
|
YPPCP1.03
NP_395229
CDS
rop
[4]
NC_003132
[4355,
4780]
+
|
YPPCP1.04
NP_395230
CDS
pim
[5]
NC_003132
[4815,
5888]
-
|
YPPCP1.05c
NP_395231
CDS
pst
[6]
NC_003132
[6006,
6422]
+
|
YPPCP1.06
NP_395232
CDS
[7]
NC_003132
[6665,
7603]
+
|
YPPCP1.07
NP_395233
CDS
pla
[8]
NC_003132
[7790,
8089]
-
|
YPPCP1.08c
NP_395234
CDS
[9]
NC_003132
[8089,
8436]
-
|
YPPCP1.09c
NP_395235
CDS
-------
seqinfo: 1 sequence from an unspecified genome

Полученный объект DataGFF принадлежит к объектам класса GRanges, структура и пути создания которых будут рассмотрены в разделе «Организация и создание объектов класса GRanges», а возможности по их дальнейшему использованию в разделе «Использование объектов класса GRanges». Поскольку объекты такого класса являются специализированными контейнерами для хранения координат интервалов генома и описательной информации по этим интервалам, то это предоставляет широкие возможности для дальнейшего анализа импортированных данных. Однако применимость функции read.gff только к файлам формата GFF, находящихся на ftp-сервере NCBI, а также отсутствие дополнительных возможностей по выбору загружаемых интервалов генома, полей и атрибутов данных существенно ограничивают ее использование.
Первое из ограничений функции read.gff можно снять, использовав вместо нее схожую функцию readGff3 из пакета genomeIntervals (http://www.bioconductor.org/packages/release/bioc/html/genomeIntervals.html), так как она работает с любыми локально расположенными файлами формата GFF3. Однако и эта функция не предоставляет пользователю каких-то дополнительных возможностей по выбору и загрузке только необходимых частей исходного файла.

Использование функций импорта из пакета rtracklayer.
Третий путь основан на применении либо базовой функции import, либо одной из специализированных функций семейства import.gff из пакета rtracklayer. Пакет rtracklayer может быть установлен с официального сайта проекта Bioconductor с помощью стандартной процедуры:

source("https://bioconductor.org/biocLite.R")
biocLite("rtracklayer")

Для корректной работы базовой функции import обязательным условием является правильное задание значений нескольких аргументов:
q con – путь к файлу или URL-адрес файла на удаленном сервере. Путь или URL-адрес задается вместе с названием файла. Архивы файлов не обязательно предварительно распаковывать, так как функция способна самостоятельно справиться с расширениями gz, bz2 и xz.
q format – формат файла. Функция работает с форматами gff, gff1, gff2, gff3, gvf и gtf. Если формат был указан в названии файла, то этот аргумент можно проигнорировать.
q version – версия формата GFF. Аргумент может принимать значения "" (версия формата будет взята функций автоматически из директивы ##gff-version * загружаемого файла, или же будет присвоено значение “1”), "1", "2" или "3". Этот аргумент следует прописывать только в том случае, если аргумент format был задан как “gff” без указаний на версию формата.
q genome – референсный геном. Задается в виде UCSC-идентификатора (например, “hg38”), либо сохраняется по умолчанию как NA (если геном не известен или же нет необходимости в указании генома). Указывать конкретный геном есть смыл только в том случае, если предполагается получение сиквенс-информации о геноме из соответствующего пакета BSgenome или непосредственно от геномного браузера UCSC через on-line соединение.
q asRangedData – аргумент управляет типом выводимого объекта. Если asRangedData = FALSE (по умолчанию), то объектом будет GRanges, если же asRangedData = TRUE, то RangedData.
q colnames – названия полей, которые необходимо подвергнуть парсингу. Парсингу могут подвергаться как постоянные поля файлов формата GFF (такие, например, как “source” или “strand”), так и поле “attributes”, содержащееся в файлах с форматами GFF2 или GFF3. По умолчанию этот аргумент имеет значение NULL, при котором парсингу будут подвержены все поля файла.
q which – диапазон данных (интервалы в геноме), которые необходимо загрузить. Диапазон задается с помощью объекта RangesList или GRanges. По умолчанию этот аргумент имеет значение NULL, при котором будут загружены все данные из файла.
q feature.type – тип признаков, которые необходимо загрузить. По умолчанию этот аргумент имеет значение NULL, при котором загружаться будут все типы признаков.
Функции import.gff, import.gff1, import.gff2 и import.gff3 из того же пакета rtracklayer предназначены для загрузки файлов в соответствующих версиях формата GFF. Эти функции имеют те же аргументы, что и функция import, но не нуждаются в аргументах format и version.
Ниже представлен код загрузки файла формата GFF с помощью базовой функции import, а также фрагмент итогового объекта.

rm(list = ls())
setwd("D:/Transcriptome")
library(rtracklayer)
DataGFF = import(con = "D:/Transcriptome/Ensembl genes. GRCh38.p3, release 81, June 2015.gff3",
   format = "gff", version = "3", genome = NA, asRangedData = FALSE, colnames = NULL,
   which = NULL, feature.type = NULL)
print(DataGFF)
GRanges object with 2605185 ranges and 26 metadata columns:
seqnames
ranges
strand
|
source
type
score
<Rle>
<IRanges>
<Rle>
|
<factor>
<factor>
<numeric>
[1]
1
[        1,
248956422]
*
|
GRCh38
chromosome
<NA>
[2]
1
[11869,
14409]
+
|
havana
gene
<NA>
[3]
1
[11869,
14409]
+
|
havana
processed_transcript
<NA>
[4]
1
[11869,
12227]
+
|
havana
exon
<NA>
[5]
1
[12613,
12721]
+
|
havana
exon
<NA>
...
...
...
...
...
...
...
...
...
[2605181]
Y
[26626520,
26627159]
-
|
havana
processed_pseudogene
<NA>
[2605182]
Y
[26626520,
26627159]
-
|
havana
exon
<NA>
[2605183]
Y
[56855244,
56855488]
+
|
havana
gene
<NA>
[2605184]
Y
[56855244,
56855488]
+
|
havana
processed_pseudogene
<NA>
[2605185]
Y
[56855244,
56855488]
+
|
havana
exon
<NA>
-------
seqinfo: 493 sequences from an unspecified genome; no seqlengths

Полученный объект DataGFF является объектом класса GRanges. Этот объект содержит координаты всех интервалов генома, представленных в загружаемых данных, а также столбцы с метаданными. Если с помощью аргумента which заданы координаты специфических интервалов генома, то объект DataGFF будет содержать только ту часть исходных данных, которая перекрывается с выбранными координатами. Если значения аргументов colnames и feature.type выбраны по умолчанию, то количество столбцов с метаданными, а также их наполнение зависит от содержания поля “attributes” исходного файла. Однако с помощью аргументов colnames и feature.type количество столбцов и содержание метаданных можно изменить. Ниже приведен пример кода для загрузки только части данных из исходного файла в формате GFF. В этом примере мы первоначально задали интересующие нас координаты генома с помощью созданного нами объекта genomicInrevals класса GRanges.

genomicInrevals = GRanges(seqnames = Rle("1"),
         ranges = IRanges(start = c(14404, 184923), end = c(29570, 200322)),
         strand = Rle("-"))

И только после этого дали команду на загрузку файла, дополнительно указав, что нас интересуют только названия и тип таких признаков как “gene”.

DataGFF = import(con = "D:/Transcriptome/Ensembl genes. GRCh38.p3, release 81, June 2015.gff3",
    format = "gff", version = "3", genome = NA, asRangedData = FALSE,
    colnames = c("Name", "biotype"),
    which = genomicInrevals, feature.type = "gene")

print(DataGFF)
GRanges object with 3 ranges and 2 metadata columns:
seqnames
ranges
strand
|
Name
biotype
<Rle>
<IRanges>
<Rle>
|
<character>
<character>
[1]
1
[  11869,
14409]
*
|
DDX11L1
transcribed_unprocessed_pseudogene
[2]
1
[  14404,
29570]
*
|
WASH7P
unprocessed_pseudogene
[3]
1
[184923,
200322]
*
|
FO538757.2
protein_coding
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

Как мы видим, полученный объект DataGFF существенно уменьшился в размере и упростился по содержанию. На этом примере мы так же видим, что загрузка файлов формата GFF с помощью функций импорта из пакета rtracklayer обладает высокой гибкостью. Кроме того, поскольку мы имеем дело с объектом класса GRanges, то это предоставляет нам широкие возможности по дальнейшему анализу импортированных данных.

Использование функции makeTxDbFromGFF из пакета GenomicFeatures.
В завершение рассмотрим еще один путь загрузки файла формата GFF с одновременным преобразованием содержимого файла в объект класса TxDb. Объекты такого класса являются специализированными контейнерами, где хранится описание транскриптов. Благодаря оптимизированной  структуре объектов класса TxDb описательная информация может быть очень быстро извлечена и подвергнута дальнейшему анализу, что является очень важным преимуществом таких объектов. Подробно о самих объектах TxDb и других способах их создания будет рассказано в разделе «Организация и создание объектов класса TxDb», а возможности по дальнейшему анализу транскриптов, хранящихся в объектах TxDb, будут описаны в разделе «Использование объектов класса TxDb».
Файл в формате GFF может быть загружен и преобразован в объект класса TxDb с помощью специализированной функции makeTxDbFromGFF из пакета GenomicFeatures. Следует обратить внимание, что функция makeTxDbFromGFF добавлена в пакет GenomicFeatures версии 1.20.2 и выше, скомпилированный под R версии 3.2.*, в более ранних же пакетах она отсутствует. Как и в случае с другими пакетами, упомянутыми в этом разделе, пакет GenomicFeatures может быть установлен с официального сайта проекта Bioconductor стандартным способом:

source("https://bioconductor.org/biocLite.R")
biocLite("GenomicFeatures")

Функция makeTxDbFromGFF имеет ряд аргументов, часть из которых является обязательной к заполнению.
q file – путь к файлу или URL-адрес файла на удаленном сервере. Путь или URL-адрес задается вместе с названием файла.
q format – формат файла. Функция работает только с форматами gff3 и gtf. По умолчанию этому аргументу присвоено значение “auto”, что позволяет в автоматическом режиме определять формат загружаемого файла.
q dataSource – краткое описание происхождения данных.
q organism – род и вид организма. Рекомендуется использовать общепринятые научные названия, например, "Homo sapiens" или "Canis familiaris", но не "human" или "my fuzzy buddy", так как эта информация может быть использована в дальнейшем анализе.
q circ_seqs – символьный вектор с перечнем кольцевых хромосом.
q chrominfo – объект типа data frame или объект класса Seqinfo, содержащий информацию о хромосомах, представленных в загружаемом файле.
q miRBaseBuild – информация о сборке пакета mirbase.db. Этот пакет основан на данных из базы miRBase (http://www.mirbase.org/), а информация о нем необходима для корректного отображения строк файла формата GFF, отведенных под микроРНК. Следует заметить, что для операционной системы Windows на официальном сайте проекта Bioconductor нет бинарного варианта пакета mirbase.db, а есть только архив с исходными кодами.
Ниже представлен код загрузки файла формата GFF с помощью функции makeTxDbFromGFF и краткое описание полученного объекта TxDb. Обращаем внимание читателей на то, что создание объекта TxDb не является быстрым процессом и требовательно к техническим характеристикам вычислительной техники. Так, при работе на компьютере, оснащенном четырехядерным процессором Intel® Corei5-4670K с тактовой частотой 3,4 GHz и 8,07 Гб оперативной памяти, в среде операционной системы Windows 7 Ultimateна загрузка и преобразование файла формата GFF3 размером 838,4 Мб, содержащего 107-й релиз базы данных NCBI RefSeq (GRCh38.p2 сборка генома человека, март 2015 года), идет чуть более 8 минут и используется 4,3 Гб оперативной памяти.

rm(list = ls())
setwd("D:/Transcriptome")
library(GenomicFeatures)
t1 = Sys.time()
TxDb = makeTxDbFromGFF(file = "D:/Transcriptome/NCBI RefSeq genes. GRCh38.p2, release 107, 12 March 2015.gff3",
    format = "auto", dataSource = "NCBI RefSeq", organism = "Homo sapiens",
    circ_seqs = DEFAULT_CIRC_SEQS, chrominfo = NULL, miRBaseBuild = NA)
t2 = Sys.time()
print(t2 – t1)
Time difference of 8.019134 mins
print(TxDb)
TxDb object:
# Db type: TxDb
# Supporting package: GenomicFeatures
# Data source: NCBI RefSeq
# Organism: Homo sapiens
# miRBase build ID: NA
# Genome: NA
# transcript_nrow: 153453
# exon_nrow: 1584842
# cds_nrow: 1232083
# Db created by: GenomicFeatures package from Bioconductor
# Creation time: 2015-09-15 15:58:53 +0300 (Tue, 15 Sep 2015)
# GenomicFeatures version at creation time: 1.20.2
# RSQLite version at creation time: 1.0.0
# DBSCHEMAVERSION: 1.1

Как уже было сказано, возможности дальнейшего анализа транскриптов, хранящихся в объектах TxDb, будут описаны в разделе «Использование объектов класса TxDb».


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

Отправить комментарий