select 唯一映射到 SNP 的基因

Question

我有一个 vcf，想要 select 100 个基因，并且对于每个基因，select 一个 SNP？从技术上讲，如果我们考虑一个基因，它有许多 rsid 映射到它。对于我的分析，我需要 select 100 个基因，每个基因 select 只有一个 SNP，忽略其他基因，并有一个最终的 vcf 文件，其中只有 100 个基因唯一映射到 1 个 rsid。有没有可以做到这一点的软件包？

22 16050075 rs587697622 A G 100.0 PASS
AC=1;AF=0.000199681;AN=5008;NS=2504;DP=8012;EAS_AF=0;AMR_AF=0;AFR_AF=0;EUR_AF=0;SAS_A F=0.001;AA=.|||;VT=SNP;ANN=G|intergenic_region|MODIFIER|CHR_START-LA16c-4G1.3|CHR_START-ENSG00000233866|intergenic_region|CHR_START-ENSG00000233866|||n.16050 075A>G|||||| GT

这就是我想要做的：图书馆（vcfR）

# Import VCF
my.vcf <- read.vcfR('my.vcf.gz')
# Combine CHROM thru FILTER cols + INFO cols
my.vcf.df <- cbind(as.data.frame(getFIX(my.vcf)), INFO2df(my.vcf))
data_new <- my.vcf.df[!duplicated(my.vcf.df[ ,"ANN"]),]

这是我的数据集的样子：

如您所见，ANN 列有多个 ENSG00000233866 值。我想根据作为基因 ID 的 ANN 值过滤 df。因此，对于 ENSG00000233866，我只想 select 只有一个 rsid 或 rs 值，而忽略其他值。我想对其他基因做同样的事情。

Answer 1

考虑 ANN 中 pipe-delimited 列的 str.split 的 gene 大约是第四个管道。然后运行基因分组的顺序计数 ave。然后subset组号为1。下面请不要使用新管道，|> in R 4.1.0+:

data_new <- within(
  data_new, {
    GENE <- gsub(
       "CHR_START-", "", 
       sapply(strsplit(ANN, "|"), "[", 4)
    )
    GENE_GROUP_SEQ <- ave(1:nrow(data_new), GENE, FUN=seq_along)
  }
) |> subset(
  GENE %in% unique(GENE)[1:100] &  # SUBSET FOR FIRST 100 GENES
  GENE_GROUP_SEQ == 1              # SUBSET FOR FIRST OF EACH GENE GROUP
) |> `row.names<-`(NULL)           # RESET ROW NAMES

对于低于 4.1.0 的 R 版本，断开管道以单独调用：

data_new <- within(
  data_new, {
    GENE <- gsub(
       "CHR_START-", "", 
       sapply(strsplit(ANN, "|"), "[", 4)
    )
    GENE_GROUP_SEQ <- ave(1:nrow(data_new), GENE, FUN=seq_along)
  }
)

data_sub <- subset(
  data_new,
  GENE %in% unique(GENE)[1:100] &           # SUBSET FOR FIRST 100 GENES
  GENE_GROUP_SEQ == 1                       # SUBSET FOR FIRST OF EACH GENE GROUP
)

dat_sub <- `row.names<-`(data_sub, NULL)    # RESET ROW NAMES

select 唯一映射到 SNP 的基因

select uniquely mapped genes to SNP

code-generation

r

filter

genome

vcf-variant-call-format