Ошибка при извлечении нескольких последовательностей из файла .fasta с использованием объекта геномного диапазона в качестве ссылки
У меня есть файл fasta, соответствующий моему эталонному геному, и файл vcf, который соответствует вызову SNP моих данных. Я хотел бы получить последовательность каждого из SNP из моего фаста. Для этого с помощью R I загрузил файл 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 (или расширенные окна) не выходят за пределы моего фаста, используя эту команду:
gr = as(seqinfo(fa), "GRanges")
findOverlaps(gr, SNPrange)
И, наконец, я запускаю команду, чтобы получить последовательность из SNPrange, используя мой файл fasta. Однако я получил следующую ошибку:
seq_ <-getSeq(fa, SNPrange)
Error in value[[3L]](cond) :
record 12177 (chr7:88167221-88167221) failed
file: F1.fasta
Я заметил, что у других людей была такая же проблема, но ни у кого не было решения, поэтому я попытался решить ее по-своему. Я попытался получить Seq для каждой хромосомы отдельно:
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
Я подозреваю, что проблема могла быть в том, что позиции, из которых я хочу извлечь строку, были "N" в моем файле fasta. Итак, я попытался найти одну из этих позиций в моем файле fasta, где R показал мне ошибку. И, к моему удивлению, это были не "Ns", а гетерозиготные основания. Однако, когда я разбил свои данные на хромосомы, алгоритм смог идентифицировать Y (C / T) и другие гетерозиготные основания, то есть у него нет проблем с вырожденными основаниями. Так что я думаю, что проблема в алгоритме, а не в моих данных. Я использовал следующую команду в bash, чтобы извлечь последовательность из желаемой позиции файла fasta:
samtools faidx F1.fasta chr7:88167220-88167222"
>chr7:88167220-88167222
> CRA
А вот и моя сессияИнфо
> 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