Проанализировать INFO файла VCF на фрейме данных R - PullRequest
Новая
       38

Проанализировать INFO файла VCF на фрейме данных R

1 голос
/ 04 марта 2020

Я пытаюсь создать фрейм данных из файла vcf, включая только некоторые элементы из поля INFO. Проблема в том, что значения этих элементов не всегда находятся в одной и той же позиции, поэтому, когда я загружаю поле VCF и разделенное INFO, я получаю эти указанные c элементы в разных столбцах.

Например: 

Pos         Score       Strand     Length     
CIPOS=0     SCORE=1     STRAND=+   LEN=634
SCORE=89    STRAND=-  LEN=567      UTR=+
CIPOS=9     SCORE=1     STRAND=+   LEN=0
CIPOS=8     SCORE=1     STRAND=+   LEN=1
STRAND=+    LEN=555     UTR=+      B

Как видите, некоторые строки сдвинуты, потому что в vcf отсутствует символ отсутствия какого-либо элемента INFO, а информация о поле читается как строка, поэтому при разбиении я не делаю знать, как сказать R записать NA в соответствующую строку каждого столбца.

Есть ли способ записать каждое значение "SCORE =" в столбце Score, каждое значение "STRAND =" в столбце Strand, и т. д. c?

Заранее спасибо!

Ответы [ 2 ]

1 голос
/ 04 марта 2020

Для этого предусмотрены пакеты, например, VariantAnnotation от Bioconductor. После прочтения в файле vcf информация упаковывается в файл data.frame, и вы можете оценить его, как показано ниже: 

library(VariantAnnotation)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
info(vcf)

DataFrame with 10376 rows and 22 columns
                 LDAF   AVGPOST       RSQ     ERATE     THETA
            <numeric> <numeric> <numeric> <numeric> <numeric>
rs7410291      0.3431     0.989    0.9856     0.002     5e-04
rs147922003    0.0091    0.9963    0.8398     5e-04    0.0011
rs114143073    0.0098    0.9891    0.5919     7e-04     8e-04
rs141778433    0.0062     0.995    0.6756     9e-04     3e-04
rs182170314    0.0041    0.9981    0.7909     7e-04     4e-04
...               ...       ...       ...       ...       ...
rs187302552     9e-04    0.9992    0.5571     3e-04    0.0026
rs9628178      0.0727    0.9997    0.9983     3e-04    0.0011
rs5770892      0.3678    0.9868    0.9784    0.0021     7e-04
rs144055359    0.0011    0.9987    0.5323     5e-04     4e-04
rs114526001    0.0543    0.9622    0.7595    0.0021     5e-04
                    CIEND         CIPOS       END        HOMLEN
            <IntegerList> <IntegerList> <integer> <IntegerList>
rs7410291           NA,NA         NA,NA        NA              
rs147922003         NA,NA         NA,NA        NA              
rs114143073         NA,NA         NA,NA        NA              
rs141778433         NA,NA         NA,NA        NA              
rs182170314         NA,NA         NA,NA        NA              
...                   ...           ...       ...           ...
rs187302552         NA,NA         NA,NA        NA              
rs9628178           NA,NA         NA,NA        NA              
rs5770892           NA,NA         NA,NA        NA              
rs144055359         NA,NA         NA,NA        NA              
rs114526001         NA,NA         NA,NA        NA              
                     HOMSEQ     SVLEN      SVTYPE            AC
            <CharacterList> <integer> <character> <IntegerList>
rs7410291                          NA          NA           751
rs147922003                        NA          NA            20
rs114143073                        NA          NA            20
rs141778433                        NA          NA            12
rs182170314                        NA          NA             8
...                     ...       ...         ...           ...
rs187302552                        NA          NA             1
rs9628178                          NA          NA           158
rs5770892                          NA          NA           801
rs144055359                        NA          NA             1
rs114526001                        NA          NA           113
                   AN          AA        AF    AMR_AF    ASN_AF
            <integer> <character> <numeric> <numeric> <numeric>
rs7410291        2184           N      0.34       0.2      0.19
rs147922003      2184           c      0.01      0.01        NA
rs114143073      2184           G      0.01    0.0028      0.02
rs141778433      2184           C      0.01      0.01        NA
rs182170314      2184           C    0.0037      0.01        NA
...               ...         ...       ...       ...       ...
rs187302552      2184           a     5e-04        NA    0.0017
rs9628178        2184           a      0.07      0.03      0.01
rs5770892        2184           a      0.37      0.32      0.38
rs144055359      2184           g     5e-04        NA        NA
rs114526001      2184           g      0.05      0.01      0.01
               AFR_AF    EUR_AF          VT       SNPSOURCE
            <numeric> <numeric> <character> <CharacterList>
rs7410291        0.83      0.22         SNP          LOWCOV
rs147922003      0.02      0.01         SNP          LOWCOV
rs114143073      0.01      0.01         SNP          LOWCOV
rs141778433      0.02        NA         SNP          LOWCOV
rs182170314      0.01        NA         SNP          LOWCOV
...               ...       ...         ...             ...
rs187302552        NA        NA         SNP          LOWCOV
rs9628178        0.17      0.08         SNP          LOWCOV
rs5770892        0.59      0.23         SNP          LOWCOV
rs144055359        NA    0.0013         SNP          LOWCOV
rs114526001      0.16      0.04         SNP          LOWCOV

В этом примере вы можете выполнить преобразование в файл data.frame и оценить столбцы. нет информации о цепях, но она должна работать, если у вас есть цепочка: 

df = as.data.frame(info(vcf))
df$CIPOS
0 голосов
/ 04 марта 2020

Я не уверен, что это самый простой способ сделать это, но он должен работать.

Прежде всего, попробуйте опубликовать свои данные, используя dput(yourdata). Тогда людям будет проще работать с вашим конкретным c примером.

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

Данные: 

a <- data.frame(Pos = c("CIPOS=0", 
                            "SCORE=89",
                            "CIPOS=9",
                            "CIPOS=8", 
                            "STRAND=+"),
                    Score = c("SCORE=1",
                              "STRAND=-", 
                              "SCORE=1", 
                              "SCORE=1", 
                              "LEN=555"),
                    Strand = c("STRAND=+",
                               "LEN=567",  
                               "STRAND=+",
                               "STRAND=+",
                               "UTR=+")
    )

dput(a)
    #> structure(list(Pos = structure(c(1L, 4L, 3L, 2L, 5L), .Label = c("CIPOS=0", 
    #> "CIPOS=8", "CIPOS=9", "SCORE=89", "STRAND=+"), class = "factor"), 
    #>     Score = structure(c(2L, 3L, 2L, 2L, 1L), .Label = c("LEN=555", 
    #>     "SCORE=1", "STRAND=-"), class = "factor"), Strand = structure(c(2L, 
    #>     1L, 2L, 2L, 3L), .Label = c("LEN=567", "STRAND=+", "UTR=+"
    #>     ), class = "factor")), class = "data.frame", row.names = c(NA, 
    #> -5L))

А затем, если вы запустите этот код: 

library(tidyverse)

    a %>% mutate(sample = row_number()) %>% 
      pivot_longer(-sample) %>%
      # here you have to unselect the name column, since this is actually wrong
      select(-name) %>% 
      # separate the strings that contain the data
      separate(value, into = c("group", "value"), sep = "=") %>% 
      # put it back to wide format
      pivot_wider( names_from = group, 
                   values_from = value)
    #> # A tibble: 5 x 6
    #>   sample CIPOS SCORE STRAND LEN   UTR  
    #>    <int> <chr> <chr> <chr>  <chr> <chr>
    #> 1      1 0     1     +      <NA>  <NA> 
    #> 2      2 <NA>  89    -      567   <NA> 
    #> 3      3 9     1     +      <NA>  <NA> 
    #> 4      4 8     1     +      <NA>  <NA> 
    #> 5      5 <NA>  <NA>  +      555   +

Создано в 2020 г. -03-04 * Представить пакет (v0.3.0)

...