【问题标题】:Replace numeric genotype code with DNA letter in VCF file用 VCF 文件中的 DNA 字母替换数字基因型代码
【发布时间】:2019-06-06 15:13:15
【问题描述】:

这里我有一个修改过的 VCF 文件,其中每个测序行的基因型以数字表示,对应的 SNP 字母在 SNP 列的列表中(我通过聚合 Ref 和 A​​lt 列生成了此列)。

使用 dplyr,我想用 SNP 字母替换数值。所有位置都是纯合的,SNP列表是按数字基因型的顺序排列的,所以使用数字基因型的第一个数字(加1)应该给出对应SNP字母的列表索引。像这样:

                  Ref   1st Alt   2nd Alt   3rd Alt
SNP letter:       C     A         T         G
Numeric genotype: 0     1         2         3
List index:       1     2         3         4

修改后的 VCF 数据框:

CHR   POS   SNP              Line1   Line2   Line3
01    10    c("A", "G")      0|0     1|1     0|0
01    20    c("C", "T", "A") 2|2     0|0     1|1
02    15    c("G", "T")      1|1     0|0     1|1

期望的输出:

CHR   POS   SNP              Line1   Line2   Line3
01    10    c("A", "G")      A       G       A
01    20    c("C", "T", "A") A       C       T
02    15    c("G", "T")      T       G       T

到目前为止,我已经尝试过这样的事情:

VCF %>%
    rowwise() %>%
    mutate_at(vars(4:ncol(.)), funs(str_replace(., "^(\\d)|\\d", SNP[[1]]["\\1"+1])))

但还没有成功。

提前感谢您的帮助。

结构是:

VCF <- structure(list(CHR = c("01", "01", "01"), POS = c(29463, 
29517, 29522), SNP = list(c("T", "C"), c("C", "G", "A"), c("T", 
"C")), PI548298 = c("0|0", "0|0", "1|1"), PI548488 = c("0|0", "0|0", 
"0|0"), PI548348 = c("0|0", "0|0", "1|1")), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -3L))

【问题讨论】:

  • 1) 您能否使用dput 发布您的数据,以便我们使用您实际拥有的数据?例如,SNP 真的是包含 R 表达式的字符串,而不是实际的向量或列表吗? 2) 在 3 行列中,是否存在杂合 SNP?我们可以忽略每个“行”列中的第二个数字吗?如果没有,你想怎么处理?
  • 已编辑以包含 dput。是的,所有位置都是纯合的,所以第二个数字可以忽略。

标签: r dplyr


【解决方案1】:

我们可以通过使用 substr 来获取每个 PI 列的第一个字符(因为这是我们需要的所有信息),将其转换为数字,加 1(因为 R 索引从 1 开始) ,然后使用它对SNP 列进行子集化。

通过使用rowwise,我们将这个函数单独应用于每一行并利用该行的SNP向量:

library(tidyverse)

VCF %>%
    rowwise() %>%
    mutate_at(vars(starts_with('PI')),
              list(~ SNP[as.numeric(substr(., 0, 1)) + 1]))

Source: local data frame [3 x 6]
Groups: <by row>

# A tibble: 3 x 6
  CHR     POS SNP       PI548298 PI548488 PI548348
  <chr> <dbl> <list>    <chr>    <chr>    <chr>   
1 01    29463 <chr [2]> T        T        T       
2 01    29517 <chr [3]> C        C        C       
3 01    29522 <chr [2]> C        T        C  

【讨论】:

  • 太完美了。谢谢@divibisan。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2012-06-24
  • 1970-01-01
相关资源
最近更新 更多