将fasta文件组织成数据框或表格

Question

我需要整理本网站的fasta文件 https://www.uniprot.org/uniprot/?query=cytokines&sort=score 到一个数据框或表格中，其中包含列 ID、名称物种和每个蛋白质的序列请帮助我使用 R !!

Answer 1

我使用网页上的界面下载了fasta文件。然后我使用了Biocuconductor包Biostrings

install.packages("BiocManager")
BiocManager::install("Biostrings")

将Fasta文件读为“Aastringset”

aa <- readAAStringSet("uniprot-cytokines.fasta.gz")

names(aa)是在Fasta文件中报告的每个序列的名称。

> head(names(aa), 3)
[1] "sp|Q8IV20|LACC1_HUMAN Purine nucleoside phosphorylase LACC1 OS=Homo sapiens OX=9606 GN=LACC1 PE=1 SV=1"          
[2] "sp|Q9H257|CARD9_HUMAN Caspase recruitment domain-containing protein 9 OS=Homo sapiens OX=9606 GN=CARD9 PE=1 SV=2"
[3] "sp|Q8BZT5|LRC19_MOUSE Leucine-rich repeat-containing protein 19 OS=Mus musculus OX=10090 GN=Lrrc19 PE=1 SV=1"  

我使用 dplyr 和 tidyr 包通过正则表达式提取相关字段

fields <-
    names(aa) |>
    dplyr::as_tibble() |>
    tidyr::extract(
        "value",
        into = c(
            "xx", "entry", "entry_name", "protien_names",
            "organism", "organism_id", "gene_names"
        ),
        regex = "^(.+)\\|(.+)\\|([^ ]+) (.*) OS=(.*) OX=(.*) GN=(.*) PE=.*$"
    )

导致

> fields
# A tibble: 17,166 × 7
   xx    entry  entry_name  protien_names        organism organism_id gene_names
   <chr> <chr>  <chr>       <chr>                <chr>    <chr>       <chr>     
 1 sp    Q8IV20 LACC1_HUMAN Purine nucleoside p… Homo sa… 9606        LACC1     
 2 sp    Q9H257 CARD9_HUMAN Caspase recruitment… Homo sa… 9606        CARD9     
 3 sp    Q8BZT5 LRC19_MOUSE Leucine-rich repeat… Mus mus… 10090       Lrrc19    
 4 sp    P10145 IL8_HUMAN   Interleukin-8        Homo sa… 9606        CXCL8     
 5 sp    P09429 HMGB1_HUMAN High mobility group… Homo sa… 9606        HMGB1     
 6 sp    Q9NZH7 IL36B_HUMAN Interleukin-36 beta  Homo sa… 9606        IL36B     
 7 sp    Q8IUC6 TCAM1_HUMAN TIR domain-containi… Homo sa… 9606        TICAM1    
 8 sp    Q80UF7 TCAM1_MOUSE TIR domain-containi… Mus mus… 10090       Ticam1    
 9 sp    O35718 SOCS3_MOUSE Suppressor of cytok… Mus mus… 10090       Socs3     
10 sp    P47811 MK14_MOUSE  Mitogen-activated p… Mus mus… 10090       Mapk14    
# … with 17,156 more rows

然后我添加了之前创建的 aa 对象中的序列

> dplyr::bind_cols(fields, sequence = as.character(aa))
# A tibble: 17,166 × 8
   xx    entry entry_name protien_names organism organism_id gene_names sequence
   <chr> <chr> <chr>      <chr>         <chr>    <chr>       <chr>      <chr>   
 1 sp    Q8IV… LACC1_HUM… Purine nucle… Homo sa… 9606        LACC1      MAEAVLI…
 2 sp    Q9H2… CARD9_HUM… Caspase recr… Homo sa… 9606        CARD9      MSDYEND…
 3 sp    Q8BZ… LRC19_MOU… Leucine-rich… Mus mus… 10090       Lrrc19     MKVTRFM…
 4 sp    P101… IL8_HUMAN  Interleukin-8 Homo sa… 9606        CXCL8      MTSKLAV…
 5 sp    P094… HMGB1_HUM… High mobilit… Homo sa… 9606        HMGB1      MGKGDPK…
 6 sp    Q9NZ… IL36B_HUM… Interleukin-… Homo sa… 9606        IL36B      MNPQREA…
 7 sp    Q8IU… TCAM1_HUM… TIR domain-c… Homo sa… 9606        TICAM1     MACTGPS…
 8 sp    Q80U… TCAM1_MOU… TIR domain-c… Mus mus… 10090       Ticam1     MDNPGPS…
 9 sp    O357… SOCS3_MOU… Suppressor o… Mus mus… 10090       Socs3      MVTHSKF…
10 sp    P478… MK14_MOUSE Mitogen-acti… Mus mus… 10090       Mapk14     MSQERPT…
# … with 17,156 more rows

要获得所需的结果（知道信息都在 aa 对象中，使用它可能比完成这些步骤更有意义）。

这似乎不太对，因为有行全是NA，可能反映的名字与正则表达式不匹配……

> fields |> dplyr::count(xx, sort = TRUE)
# A tibble: 3 × 2
  xx        n
  <chr> <int>
1 tr    14644
2 sp     1336
3 NA     1186

hmm ...看起来它们缺少GN=（gene_names）元素

> names(aa)[ids]
[1] "tr|G3QVN6|G3QVN6_GORGO C-C motif chemokine OS=Gorilla gorilla gorilla OX=9595 PE=3 SV=1"     
[2] "tr|K9IQA0|K9IQA0_DESRO C-C motif chemokine (Fragment) OS=Desmodus rotundus OX=9430 PE=2 SV=1"
[3] "tr|K9IFY6|K9IFY6_DESRO C-C motif chemokine OS=Desmodus rotundus OX=9430 PE=2 SV=1" 

我不确定如何最好地处理这个问题，例如，分别解析带有和不带有基因名称的记录...

idx <- grepl("GN=", names(aa))
aa_with_gene_names <- aa[ idx ]
aa_without_gene_names <- aa[ !idx ]