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

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

Подробно структура VCF файлов анализировалась нами ранее в разделе “Парадигма картирования seed-and-vote: итоговые файлы в формате VCF” в контексте работы элайнера Subjunc. Теперь мы рассмотрим процедуру загрузки содержимого таких файлов в рабочее пространство R среды. Для этого мы возьмем один из файлов с расширением indel, генерируемых элайнером Subjunc, а также воспользуемся функцией readVcf из R/Bioconductor библиотеки VariantAnnotation. Эта функция имеет четыре аргумента:
q file – файл формата VCF. Файл может быть задан в виде символьного вектора (вместе с указанием пути к файлу), или же как объект класса TabixFile (см. ниже). Последний вариант подачи файла VCF на функцию readVcf является предпочтительным, так как он позволяет ускорить загрузку файла в рабочее пространство R среды, что особенно актуально при работе с большими исходными файлами. Более того, если планируется загрузить не весь файл, а лишь его часть, указанную с помощью аргумента param (см. далее), то путь загрузки через создание объекта класса TabixFile является единственно правильным и возможным.
q genome – идентификатор генома. Геном может быть указан в виде простого символьного вектора или же как объект класса Seqinfo.
q param – комплексный аргумент, который может быть задан как объект класса ScanVcfParam (см. ниже), GRanges, GRangesList, RangedData или RangesList. Этот аргумент позволяет загружать в рабочее пространство R среды не весь файл формата VCF, а лишь какую-то его часть.
q row.namzses – логический аргумент, который может иметь одно из двух значений: либо TRUE (по умолчанию), либо FALSE. В первом случае каждой строке данных файла формата VCF будет присвоено имя, что можно использовать в дальнейшей работе с таким файлом после его загрузки.
Объект класса TabixFile является активной ссылкой на индексированный файл формата VCF. Такого рода объект может быть создан с помощью одноименной функции TabixFile из библиотеки Rsamtools. Функция TabixFile имеет три аргумента, но обязательными являются только первых два:
q file – символьный вектор с названием индексированного файла и пути доступа к нему (путь может быть удаленным, наподобие http://... или ftp://...).
q index – символьный вектор с названием файла с индексами и пути доступа к нему.
q yieldSize – количество записей, загружаемых одномоментно из исходного файла в рабочее пространство R среды.
Наиболее гибким способом задания аргумента param функции readVcf является конструирование объекта класса ScanVcfParam, который является хранилищем информации о том, какую часть файла формата VCF следует загружать в рабочее пространство R среды. Такого рода объекты создаются с помощью одноименной функции ScanVcfParam, входящей в состав библиотеки VariantAnnotation. Функция ScanVcfParam имеет 6 основных аргументов:
q fixed – символьный вектор с названием(-ми) фиксированного (обязательного и постоянного по содержанию) поля (столбца) файла формата VCF, которое будет загружено в рабочее пространство R среды. Если этот аргумент не задан, то по умолчанию будут загружены все фиксированные поля (CHROM, POS, ID, REF, ALT, QUAL и FILTER). Если ему присвоено значение NA, то будет загружено только поле REF.
q info – символьный вектор с названием атрибута(-ов) поля INFO (если имеется) файла формата VCF, который будет загружен в рабочее пространство R среды. Если этот аргумент не задан, то будут загружены все атрибуты поля INFO, а если ему присвоено значение NA, то ничего из этого поля загружаться не будет. Чтобы узнать, какие атрибуты поля INFO доступны, можно воспользоваться функцией scanVcfHeader, подав на нее название файла в формате VCF.
q geno – символьный вектор с названием атрибута(-ов) поля GENO (если имеется) файла формата VCF, который будет загружен в рабочее пространство R среды. Если этот аргумент не задан, то будут загружены все атрибуты поля GENO, а если ему присвоено значение NA, то ничего из этого поля загружаться не будет. Как и в случае с полем INFO, чтобы узнать, какие атрибуты поля GENO доступны, можно воспользоваться функцией scanVcfHeader, подав на нее название VCF файла.
q trimEmpty – логический аргумент, который контролирует загрузку и вывод поля GENO в тех случаях, когда оно целиком или частично не заполнено в файле формата VCF. По умолчанию этому аргументу присвоено значение TRUE – пустое поле GENO игнорируется.
q samplesсимвольный вектор с названием образца(-ов) (если имеется) из файла формата VCF, данные по которому будут загружены в рабочее пространство R среды. Если этот аргумент не задан, то будут загружены данные по всем образцам, имеющимся в файле. Чтобы получить список имеющихся в файле формата VCF образцов можно предварительно сделать запрос samples(scanVcfHeader()), включив в него название интересующего нас файла.
q which – аргумент, определяющий геномные интервалы, по которым будет сделан запрос в файл формата VCF. Задается с помощью объекта класса GRanges или RangedData. Если этот аргумент не задан, то будут загружены все данные из файла.
Теперь мы можем выразить все описанное выше в виде одного консолидированного кода и загрузить нужный нам файл в рабочее пространство R среды:

rm(list = ls())
library(VariantAnnotation)
indels = readVcf(file = “D:/Transcriptome/SRR1145838.indel", genome = "hg38", row.names = FALSE)
indels
class: CollapsedVCF
dim: 454119 0
rowRanges(vcf):
  GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
info(vcf):
  DataFrame with 3 columns: INDEL, DP, SR
info(header(vcf)):

Number
Type
Description
INDEL
0
Flag
Indicates that the variant is an INDEL.
DP
1
Integer
Raw read depth
SR
1
String
Number of supporting reads for variants
geno(vcf):
  SimpleList of length 0:

В итоге мы получили объект indels класса CollapsedVCF, содержащий всю информацию из нашего файла с расширением indel, причем эта информация быстро извлекается и может подвергаться дальнейшей обработке. Так, что бы посмотреть содержимое всех фиксированных полей, достаточно воспользоваться простой командой:

rowRanges(indels)
GRanges object with 9 ranges and 5 metadata columns:
seqnames
ranges
strand
|
paramRangeID
REF
ALT
QUAL
FILTER
<Rle>
<IRanges>
<Rle>
|
<factor>
<DNAStringSet>
<DNAStringSetList>
<numeric>
<character>
[1]
chr21
[8395527,
8395529]
*
|
<NA>
GGG
GG
250
.
[2]
chr21
[8212492,
8212494]
*
|
<NA>
GGG
GG
250
.
[3]
chr21
[8253746,
8253747]
*
|
<NA>
TC
TCC
250
.
[4]
chr21
[8253746,
8253747]
*
|
<NA>
TC
TCGTCC
250
.
[5]
chr21
[8256695,
8256698]
*
|
<NA>
GGCG
GG
250
.
-------
seqinfo: 195 sequences from hg38 genome; no seqlengths

Аналогичным образом мы можем получить доступ к содержимому поля INFO:

info(indels)
DataFrame with 454119 rows and 3 columns

INDEL
DP
SR

<logical>
<integer>
<character>
1
TRUE
33888
16694
2
TRUE
99008
34531
3
TRUE
296
2
4
TRUE
22945
80
5
TRUE
52528
228

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

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

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