Question: Error using rsamtools to extract DNA sequences
0
gravatar for hwu12
15 months ago by
hwu1210
hwu1210 wrote:

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   
rsamtools • 285 views
ADD COMMENTlink modified 15 months ago • written 15 months ago by hwu1210
Answer: C: Error using rsamtools to extract DNA sequences
0
gravatar for hwu12
15 months ago by
hwu1210
hwu1210 wrote:
Found the problem! The indexFA should take the file path (e.g. here as  indexFA("~/Desktop/test_genome/test.fa")).
ADD COMMENTlink written 15 months ago by hwu1210
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 218 users visited in the last hour