Подробно структура
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
|
Всю необходимую
информацию можно получить и через слоты созданного объекта, но это не меняет
принципиально ситуацию. Дальнейшая же работа со всеми этими данными зависит от
тех задач, которые стоят перед исследованием.
Комментариев нет:
Отправить комментарий