Загрузка файла формата 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”).