R: (Pegas) problems with haplotypes - (error: 'h' must be of class 'haplotype')

Question

我最近开始研究单倍型数据，我正在处理来自 1000 基因组项目的数据，并尝试使用 R 中的 Pegas 包对其进行操作。到目前为止，我已经走到了这一步：

library(pegas)
a <- "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502"
b <- "ALL.chrY.phase3_integrated_v1b.20130502.genotypes.vcf.gz"
url <- paste(a, b, sep = "/")
download.file(url, "chrY.vcf.gz")
(info <- VCFloci("chrY.vcf.gz"))
SNP <- is.snp(info)
X.SNP <- read.vcf("chrY.vcf.gz", which.loci = which(SNP))
h <- haplotype(X.SNP, 6020:6030)
net <- haploNet(h)
plot(net)

我想绘制一个单倍型网，但它没有执行。我收到以下消息：'h' 必须属于 class 'haplotype'

如果我打印出 h 我得到：

> h
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
. "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"   "C"   "C"   "C"   "C"   "C"   "C"   "T"   "C"   "C"  
. "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"   "G"   "G"   "G"   "G"   "G"   "G"   "G"   "A"   "G"  
. "C"  "C"  "C"  "C"  "C"  "C"  "T"  "C"  "C"  "C"   "C"   "C"   "C"   "C"   "C"   "C"   "C"   "C"   "C"  
. "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"   "C"   "T"   "T"   "T"   "T"   "T"   "T"   "T"   "T"  
. "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"  "G"   "G"   "G"   "G"   "A"   "G"   "G"   "G"   "G"   "G"  
. "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"   "T"   "C"   "T"   "T"   "T"   "T"   "T"   "T"   "T"  
. "A"  "A"  "A"  "A"  "A"  "A"  "A"  "A"  "A"  "C"   "A"   "A"   "A"   "A"   "A"   "A"   "A"   "A"   "A"  
. "G"  "G"  "G"  "."  "G"  "G"  "G"  "G"  "G"  "G"   "G"   "G"   "A"   "G"   "G"   "G"   "G"   "G"   "G"  
. "."  "T"  "C"  "T"  "T"  "C"  "T"  "."  "."  "."   "T"   "T"   "T"   "T"   "C"   "T"   "T"   "T"   "T"  
. "."  "A"  "."  "A"  "."  "C"  "A"  "A"  "C"  "."   "A"   "A"   "A"   "A"   "A"   "C"   "A"   "A"   "A"  
. "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"  "T"   "T"   "T"   "T"   "T"   "T"   "T"   "T"   "T"   "C"  
attr(,"class")
[1] "haplotype.loci"
attr(,"freq")
 [1]   18 1142    2    5   25    6    1    4    2    1    2    5    1    9    1    3    1    4    1

它显然分配了 19 个单倍型。数据的呈现方式一定有问题。有什么建议吗？此外，关于 Pegas 以及如何使用 Pegas 处理 VCF 文件的内容很少 material。有没有人知道一个很好的资源（网页或书籍）来获取有关如何从 VCF 文件中处理单倍型的信息，它甚至不必用于 Pegas，任何 R 库都可以，或者 Python.. .真的什么都可以。

谢谢你的帮助，彼得

Answer 1

这是预期的结果：目前 haploNet() 仅适用于从 DNA 序列 (class "DNAbin") 生成的 class "haplotype"。 read.vcf() 的输出是 class "loci" 并且 haplotype() 是一个通用函数，适用于 classes.

如果您只处理 SNP，您可以通过以下方式避免这种情况：

class(h) <- NULL
h <- as.DNAbin(h)

（最终）目标是让 haploNet() 也与 class "haplotype.loci"（仍在开发中）和其他人一起工作。

干杯，伊曼纽尔

Answer 2

我知道这是一个旧问题 post，但如果其他人遇到同样的问题，我已经找到了解决该问题的方法。使用 pacakage "vcfR" 您可以使用 read.vcfR() 读取 vcf，然后使用 vcfR2DNAbin() 将其转换为 DNAbin。在 DNAbin 上使用 haplotype() 会导致 class "haplotype" 而不是 "haplotype.loci".

R: (Pegas) problems with haplotypes - (error: 'h' must be of class 'haplotype')

R: (Pegas) problems with haplotypes - (error: 'h' must be of class 'haplotype')

r

bioinformatics

vcf-variant-call-format