Использование базовой 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 ...
Использование функции 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® Core™ i5-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».
Комментариев нет:
Отправить комментарий