我最近开始研究单倍型数据，我正在处理来自 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”必须属于“单倍型”类
如果我打印出 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 文件的资料也很少。有没有人知道一个很好的资源(网页或书籍)来获取有关如何从 VCF 文件中操作单倍型的信息，它甚至不必用于 Pegas，任何 R 库都可以，或者 Python ......任何真的。
谢谢你的帮助，彼得
最佳答案
这是一个预期的结果:目前 haploNet() 仅适用于从 DNA seqs 生成的类“单倍型”(类“DNAbin”)。 read.vcf() 的输出属于“loci”类，haplotype() 是适用于这两个类的通用函数。
如果您只处理 SNP，您可以通过以下方式避免这种情况:
class(h) <- NULL
h <- as.DNAbin(h)
(最终)目标是让 haploNet() 也适用于“haplotype.loci”类(仍在开发中)和其他类。
干杯，伊曼纽尔
关于R:(飞马)单倍型问题 - (错误: 'h' 必须属于 'haplotype' 类)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/37015751/