Entering edit mode
Hi,
I am testing how to use rsamtools to extract DNA sequences from a fasta file. The testing fasta (as test.fa) is very simple:
>scaffold_0 AAACGCGTTAGAGGCCTAACACGAGGCAGCTTCATGCTGCAGGTCACTCTAGAGGACCCAGCGATCTGCT
I used GRanges to create a GRanges object:
gr <- GRanges(seqnames = "scaffold_0", strand = c("+"), ranges = IRanges(start=4, end = 10))
and load/index the fasta and finally use getSeq to extract the sequence:
FA <- FaFile("~/Desktop/test_genome/test.fa") indexFa(FA) getSeq(FA,gr)
However, the getSeq function gives me an error:
Error in value[[3L]](cond) : record 1 (scaffold_0:4-10) failed file: /Users/myhome/Desktop/test_genome/test.fa
Any thoughts? Thanks!!!
> sessionInfo() R version 3.4.3 (2017-11-30) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS High Sierra 10.13.4 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats4 stats graphics grDevices [6] utils datasets methods base other attached packages: [1] Rsamtools_1.30.0 Biostrings_2.46.0 [3] XVector_0.18.0 GenomicFeatures_1.30.3 [5] AnnotationDbi_1.40.0 Biobase_2.38.0 [7] GenomicRanges_1.30.3 GenomeInfoDb_1.14.0 [9] IRanges_2.12.0 S4Vectors_0.16.0 [11] BiocGenerics_0.24.0 biomaRt_2.34.2 loaded via a namespace (and not attached): [1] httr_1.3.1 RMySQL_0.10.14 [3] bit64_0.9-7 splines_3.4.3 [5] Formula_1.2-3 assertthat_0.2.0 [7] latticeExtra_0.6-28 blob_1.1.1 [9] GenomeInfoDbData_1.0.0 yaml_2.1.19 [11] progress_1.1.2 pillar_1.2.2 [13] RSQLite_2.1.1 backports_1.1.2 [15] lattice_0.20-35 digest_0.6.15 [17] RColorBrewer_1.1-2 checkmate_1.8.5 [19] colorspace_1.3-2 htmltools_0.3.6 [21] Matrix_1.2-14 plyr_1.8.4 [23] DESeq2_1.18.1 XML_3.98-1.11 [25] genefilter_1.60.0 zlibbioc_1.24.0 [27] xtable_1.8-2 scales_0.5.0 [29] BiocParallel_1.12.0 htmlTable_1.11.2 [31] tibble_1.4.2 annotate_1.56.2 [33] ggplot2_2.2.1 SummarizedExperiment_1.8.1 [35] nnet_7.3-12 lazyeval_0.2.1 [37] survival_2.42-3 magrittr_1.5 [39] memoise_1.1.0 foreign_0.8-70 [41] tools_3.4.3 data.table_1.10.4-3 [43] prettyunits_1.0.2 matrixStats_0.53.1 [45] stringr_1.3.1 munsell_0.4.3 [47] locfit_1.5-9.1 cluster_2.0.7-1 [49] DelayedArray_0.4.1 compiler_3.4.3 [51] rlang_0.2.0 grid_3.4.3 [53] RCurl_1.95-4.10 rstudioapi_0.7 [55] htmlwidgets_1.2 bitops_1.0-6 [57] base64enc_0.1-3 gtable_0.2.0 [59] DBI_1.0.0 R6_2.2.2 [61] GenomicAlignments_1.14.2 gridExtra_2.3 [63] rtracklayer_1.38.3 knitr_1.20 [65] bit_1.1-13 Hmisc_4.1-1 [67] stringi_1.2.2 Rcpp_0.12.16 [69] geneplotter_1.56.0 rpart_4.1-13 [71] acepack_1.4.1