Ответ 1
попробуйте следующее:
> sub(".*?GN=(.*?);.*", "\\1", a)
[1] "NOC2L"
У меня вопрос об извлечении части строки. Например, у меня есть строка вроде этого:
a <- "DP=26;AN=2;DB=1;AC=1;MQ=56;MZ=0;ST=5:10,7:2;CQ=SYNONYMOUS_CODING;GN=NOC2L;PA=1^1:0.720&2^1:0"
Мне нужно извлечь все между GN=
и ;
. Здесь будет NOC2L
.
Возможно ли это?
Примечание: Это форма столбца INFO
Формат файла VCF. GN - это имя гена, поэтому мы хотим извлечь имя гена из столбца INFO
.
попробуйте следующее:
> sub(".*?GN=(.*?);.*", "\\1", a)
[1] "NOC2L"
Предполагая, что точки с запятой разделяют ваши элементы, а знаки равенства встречаются исключительно между парами ключ/значение, метод не строго-регулярного выражения будет:
bits <- unlist(strsplit(a, ';'))
do.call(rbind, strsplit(bits, '='))
[,1] [,2]
[1,] "DP" "26"
[2,] "AN" "2"
[3,] "DB" "1"
[4,] "AC" "1"
[5,] "MQ" "56"
[6,] "MZ" "0"
[7,] "ST" "5:10,7:2"
[8,] "CQ" "SYNONYMOUS_CODING"
[9,] "GN" "NOC2L"
[10,] "PA" "1^1:0.720&2^1:0"
Тогда это просто вопрос выбора соответствующего элемента.
Один из способов:
gsub(".+=(\\w+);.+", "\\1", a, perl=T)
Я уверен, что есть более элегантные способы сделать это.
a <- "DP=26;AN=2;DB=1;AC=1;MQ=56;MZ=0;ST=5:10,7:2;CQ=SYNONYMOUS_CODING;GN=NOC2L;PA=1^1:0.720&2^1:0"
m = regexpr("GN.*;",a)
substr(a,m+3,m+attr(m,"match.length")-2)
Поскольку строка поступает из файла VCF, мы можем использовать VariantAnnotation package:
library(VariantAnnotation)
# read dummy VCF file
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
# see first 5 variables for info column
info(vcf)[1:3, 1:5]
# DataFrame with 3 rows and 5 columns
# LDAF AVGPOST RSQ ERATE THETA
# <numeric> <numeric> <numeric> <numeric> <numeric>
# rs7410291 0.3431 0.9890 0.9856 2e-03 0.0005
# rs147922003 0.0091 0.9963 0.8398 5e-04 0.0011
# rs114143073 0.0098 0.9891 0.5919 7e-04 0.0008
# Now extract one column, e.g.: LDAF
info(vcf)[1:3, "LDAF"]
# [1] 0.3431 0.0091 0.0098
В приведенном выше примере объект VCF отсутствует столбец "GN", но идея такая же, поэтому в вашем случае ниже должно работать:
# extract gene name
info(vcf)[, "GN"]