среда, 22 июня 2016 г.

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

Существует несколько путей загрузки данных, хранящихся в файлах формата BED, в рабочее пространство R среды.

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

rm(list = ls())
dataBED = read.table(file"D:/Transcriptome/SRR1145838.junction.bed",
          sep = "\t")
header = read.table(file = "D:/Transcriptome/SRR1145838.junction.bed",
       comment.char = ""
       nrows = 1,
       header = FALSE,
       sep = ",",
       stringsAsFactors = FALSE)
colnames(dataBED) = c(sub("#", "", header[, 1]), gsub(" ", "", header[1, -1]))
head(dataBED)
 Chr
StartLeftBlock
EndRightBlock
Junction_Name
nSupport
Strand
StartLeftBlock
EndRightBlock
Color
nBlocks
BlockSizes
BlockStarts
 chr1
14941
15883
JUNC00202937
15
-
14941
15883
0,255,255
2
97,88
0,854
 chr1
14963
186405
JUNC00167075
30
-
14963
186405
0,255,255
2
75,89
0,171353
 chr1
16237
16634
JUNC00140845
1
-
16237
16634
0,255,255
2
73,28
0,369
 chr1
16689
16943
JUNC00100248
13
-
16689
16943
0,255,255
2
76,86
0,168
 chr1
16974
17326
JUNC00214170
6
-
16974
17326
0,255,255
2
81,94
0,258

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

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

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

Для корректной работы базовой функции import обязательным условием является правильное задание значений нескольких аргументов:
q con – путь к файлу, URL-адрес файла на удаленном сервере или объект класса BEDFile. Путь или URL-адрес задается вместе с названием файла. Архивы файлов не обязательно предварительно распаковывать, так как функция способна самостоятельно справиться с расширениями gz, bz2 и xz;
q format – формат файла. Функция работает с форматами bed, bed15, bedGraph или bedpe. Если формат был указан в названии файла, то этот аргумент можно проигнорировать;
q text – символьный вектор с данными, которые необходимо загрузить в рабочее пространство R среды если аргумент con отсутствует;
q trackLine – аргумент, используемый для загрузки строки графических атрибутов. Может принимать два логических значения: TRUE или FALSE;
q genome – референсный геном. Задается в виде UCSC-идентификатора (например, “hg38”), либо сохраняется по умолчанию как NA (если геном не известен или же нет необходимости в указании генома). Указывать конкретный геном есть смыл только в том случае, если предполагается получение сиквенс-информации о геноме из соответствующего пакета BSgenome или непосредственно от геномного браузера UCSC через on-line соединение;
q colnames – названия полей, которые необходимо подвергнуть парсингу. По умолчанию этот аргумент имеет значение NULL, при котором парсингу будут подвержены все поля файла;
q which – диапазон данных (интервалы в геноме), которые необходимо загрузить. Диапазон задается с помощью объекта класса GRanges или объектов других классов, для которых применима функция findOverlaps. По умолчанию этот аргумент имеет значение NULL, при котором будут загружены все данные из файла;
q seqinfo – объект класса Seqinfo, содержащий информацию о геноме. По умолчанию этот аргумент имеет значение NULL;
q extraCols – символьный вектор, содержащий названия и классы дополнительных полей (если присутствуют) в загружаемом BED файле. Этот аргумент позволяет работать с BED файлами, которые помимо стандартного набора обязательных и необязательных полей содержат любые дополнительные поля, введенные пользователем. Аргумент задается тем же способом, что и аргумент colClasses функции read.table.
Функция  import.bed из того же пакета rtracklayer предназначена для загрузки только BED файлов и не нуждается в аргументе formatНиже представлен код загрузки BED файла с помощью базовой функции import, а также фрагмент итогового объекта.

rm(list = ls())
library(rtracklayer)
dataBED = import(con = "D:/Transcriptome/SRR1145838.junction.bed",
    format = "bed",
    trackLine = FALSE,
    genome = "hg38",
    colnames = NULL,
    which = NULL,
    seqinfo = NULL)
dataBED
GRanges object with 245354 ranges and 5 metadata columns:
seqnames
ranges
strand
|
name
score
itemRgb
thick
blocks
<Rle>
<IRanges>
<Rle>
|
<character>
<numeric>
<character>
<IRanges>
<IRangesList>
[1]
chr1
[14942,
15883]
-
|
JUNC00202937
15
#00FFFF
[14942,15883]
[1,97][855,942]
[2]
chr1
[14964,
186405]
-
|
JUNC00167075
30
#00FFFF
[14964,186405]
[1,75][171354,171442]
[3]
chr1
[16238,
16634]
-
|
JUNC00140845
1
#00FFFF
[16238,16634]
[1,73][370,397]
[4]
chr1
[16690,
16943]
-
|
JUNC00100248
13
#00FFFF
[16690,16943]
[1,76][169,254]
[5]
chr1
[16975,
17326]
-
|
JUNC00214170
6
#00FFFF
[16975,17326]
[1,81][259,352]
-------
seqinfo: 455 sequences (1 circular) from hg38 genome

Полученный объект dataBED является объектом класса GRanges. Этот объект содержит координаты всех интервалов генома, представленных в BED файле, а также столбцы с метаданными. Если же с помощью аргумента which были бы заданы координаты специфических интервалов генома, то объект dataBED содержал бы только ту часть исходных данных, которая перекрывалась бы с выбранными координатами. Дополнительную гибкость при загрузке данных из BED файла c помощью функции import придают такие аргументы, как trackLine, colnames и extraCols. Учитывая также, что мы имеем дело с объектом класса GRanges, который предоставляет широкие возможности по дальнейшему анализу импортированных данных, загрузка BED файла с помощью функции import выглядит куда предпочтительнее, чем с помощью функции read.table.
Следует обратить внимание, что при загрузке данных BED файла с помощью функции import данные загружаются не “как есть”, а в модифицированном виде (для этого следует сравнить объект dataBED, созданный нами двумя разными способами):
1) геномные координаты указаны в “правильном” виде – первый нуклеотид хромосомальной ДНК имеет координату 1, а не 0;
2) RGB кодировка цвета заменена на шестнадцатеричную кодировку;
3) два поля thickStart и thickEnd объединены в один столбец метаданных thick класса IRanges;
4) три поля blockCount, blockSizes и blockStarts сжаты в один столбец метаданных blocks класса CompressedIRangesList.
На последней особенности остановимся подробнее. Формально столбик blocks является списком интервалов:

dataBED$blocks
IRangesList of length 245354
[[1]]
IRanges object with 2 ranges and 0 metadata columns:

start
end
width

<integer>
<integer>
<integer>
 [1]
1
97
97
 [2]
855
942
88
[[2]]
IRanges object with 2 ranges and 0 metadata columns:

start
end
width

<integer>
<integer>
<integer>
 [1]
1
75
75
 [2]
171354
171442
89
[[3]]
IRanges object with 2 ranges and 0 metadata columns:

start
end
width

<integer>
<integer>
<integer>
 [1]
1
73
73
 [2]
370
397
28
...
<245351 more elements>

Длина такого списка равна количеству строк в BED файле – во всем файле или же в той его части, которая перекрывает запрошенный через аргумент which участок генома. В нашем случае аргумент which не был задан, поэтому были загружены все 245354 строки BED файла. Каждый элемент списка является объектом класса IRanges, длина которого (количество включенных интервалов) равна соответствующему значению поля blockCount из BED файла (в нашем случае это всегда 2). По сути же каждый элемент списка содержит координаты и длины блоков перекрывающихся фрагментов, которые принадлежат ридам одного и того же сплайсингового события (более подробно о блоках написано в разделе “Парадигма картирования seed-and-vote: итоговые файлы в формате BED”).

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

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