我有一个对应于我的参考基因组的 fasta 文件和一个对应于我的数据的 SNP 调用的 vcf 文件。我想从我的 fasta 中获取每个 SNP 的序列。为此，使用 RI 加载 vcf 文件并使用以下命令从中提取基因组范围对象：
vcf.fn<-"SNPsAcrossAlltheIndividuals.vcf"
vcf <- readVcf(vcf.fn, verbose=FALSE)
SNPrange <- vcf@rowRanges
我将 SNP 的位置扩展到每侧一个碱基，但我不会在这里考虑它，因为它会给我的问题增加更多偏见。之后也使用 R，我使用 Rsamtools 包使用以下命令读取 fasta：
library(Rsamtools)
file_path <- "F1.fasta"
indexFa(file_path)
fa = FaFile("F1.fasta")
我使用以下命令检查了所有 SNP（或扩展窗口）是否没有超出我的 fasta 边界：
gr = as(seqinfo(fa), "GRanges")
findOverlaps(gr, SNPrange)
最后，我运行命令 do 使用我的 fasta 文件从 SNPrange 获取序列。但是我收到以下错误：
seq_ <-getSeq(fa, SNPrange)
Error in value[[3L]](cond) :
record 12177 (chr7:88167221-88167221) failed
file: F1.fasta
我注意到其他人也有同样的问题，但没有一个解决方案，所以我试图用我的方式解决它。我试图分别为每个染色体获取序列：
chr1<- gr[seqnames(gr) == "chr1" ]
chr2<- gr[seqnames(gr) == "chr2" ]
chr3<- gr[seqnames(gr) == "chr3" ]
...
seq1 <-getSeq(fa, chr1)
seq2 <-getSeq(fa, chr2)
seq3 <-getSeq(fa, chr3)
...
并且有效，但是有一些染色体存在完全相同的问题：
seq7 <-getSeq(fa, chr7)
Error in value[[3L]](cond) : record 993 (chr7:88167220-88167222) failed
file: F1.fasta
我怀疑问题可能是我要从中提取字符串的位置是我的 fasta 文件中的“N”。所以我试图在我的 fasta 文件中找到其中一个位置，其中 R 向我显示了一个错误。令我惊讶的是，它们不是“Ns”，而是杂合碱基。然而，当我通过染色体对我的数据进行子集化时，该算法能够识别出一个 Y（C/T）和其他杂合碱基，也就是说，它对简并碱基没有问题。所以我认为问题在于算法而不是我的数据。我在 bash 中使用以下命令从 fasta 文件的所需位置提取序列：
samtools faidx F1.fasta chr7:88167220-88167222"
>chr7:88167220-88167222
> CRA
这是我的 sessinInfo
> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=Swedish_Sweden.1252 LC_CTYPE=Swedish_Sweden.1252 LC_MONETARY=Swedish_Sweden.1252
[4] LC_NUMERIC=C LC_TIME=Swedish_Sweden.1252
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] readr_1.3.1 limma_3.44.3
[3] ggplot2_3.3.2 stringr_1.4.0
[5] vcfR_1.12.0 adegenet_2.1.3
[7] ape_5.4-1 ade4_1.7-15
[9] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0 GenomicFeatures_1.40.1
[11] AnnotationDbi_1.50.3 VariantAnnotation_1.34.0
[13] Rsamtools_2.4.0 Biostrings_2.56.0
[15] XVector_0.28.0 SummarizedExperiment_1.18.2
[17] DelayedArray_0.14.1 matrixStats_0.56.0
[19] Biobase_2.48.0 GenomicRanges_1.40.0
[21] GenomeInfoDb_1.24.2 IRanges_2.22.2
[23] S4Vectors_0.26.0 BiocGenerics_0.34.0
loaded via a namespace (and not attached):
[1] colorspace_1.4-1 seqinr_3.6-1 deldir_0.1-29 ellipsis_0.3.1
[5] class_7.3-17 rstudioapi_0.11 farver_2.0.3 bit64_4.0.5
[9] fansi_0.4.1 xml2_1.3.2 codetools_0.2-16 splines_4.0.2
[13] memuse_4.1-0 cluster_2.1.0 dbplyr_2.0.0 shiny_1.5.0
[17] compiler_4.0.2 httr_1.4.2 assertthat_0.2.1 Matrix_1.2-18
[21] fastmap_1.0.1 cli_2.1.0 later_1.1.0.1 htmltools_0.5.0
[25] prettyunits_1.1.1 tools_4.0.2 igraph_1.2.5 coda_0.19-4
[29] gtable_0.3.0 glue_1.4.2 GenomeInfoDbData_1.2.3 reshape2_1.4.4
[33] dplyr_1.0.2 rappdirs_0.3.1 gmodels_2.18.1 Rcpp_1.0.5
[37] raster_3.3-13 vctrs_0.3.4 spdep_1.1-5 gdata_2.18.0
[41] nlme_3.1-148 rtracklayer_1.48.0 pinfsc50_1.2.0 mime_0.9
[45] lifecycle_0.2.0 gtools_3.8.2 XML_3.99-0.3 LearnBayes_2.15.1
[49] zlibbioc_1.34.0 MASS_7.3-51.6 scales_1.1.1 BSgenome_1.56.0
[53] hms_0.5.3 promises_1.1.1 expm_0.999-5 curl_4.3
[57] memoise_1.1.0 biomaRt_2.44.4 stringi_1.5.3 RSQLite_2.2.1
[61] e1071_1.7-3 permute_0.9-5 boot_1.3-25 BiocParallel_1.22.0
[65] spData_0.3.8 rlang_0.4.7 pkgconfig_2.0.3 bitops_1.0-6
[69] lattice_0.20-41 purrr_0.3.4 sf_0.9-6 labeling_0.4.2
[73] GenomicAlignments_1.24.0 bit_4.0.4 tidyselect_1.1.0 plyr_1.8.6
[77] magrittr_1.5 R6_2.5.0 generics_0.1.0 DBI_1.1.0
[81] withr_2.3.0 mgcv_1.8-31 pillar_1.4.6 units_0.6-7
[85] RCurl_1.98-1.2 sp_1.4-2 tibble_3.0.3 crayon_1.3.4
[89] KernSmooth_2.23-17 BiocFileCache_1.12.1 progress_1.2.2 grid_4.0.2
[93] blob_1.2.1 vegan_2.5-6 digest_0.6.25 classInt_0.4-3
[97] xtable_1.8-4 httpuv_1.5.4 openssl_1.4.3 munsell_0.5.0
[101] viridisLite_0.3.0 askpass_1.1