Hi all

I'm dealing with a large FASTA file that contains 2500 contings.

I'm trying to find which contings have a specific DNA sequence and get its coordinates in the specific contig.

when using vmatchPattern I get Granges object only for the 3 first contings.

I would like to get only the contings which have the specific sequence

Thank you

Avner

here are the problematic code lines :

heb_DNAset <- readDNAStringSet("GCA_027704765.1_ASM2770476v1_genomic.fna")

match <- vmatchPattern("ACTCAGGGCGCTTTAGGATG",heb_DNAset )

match 


# include your problematic code here with any corresponding output 
# please also include the results of running the following in an R session 

sessionInfo( )

suppressMessages({  
  library(Biostrings)
  library(magrittr)
})
#> Warning: Paket 'IRanges' wurde unter R Version 4.3.1 erstellt
#> Warning: Paket 'GenomeInfoDb' wurde unter R Version 4.3.1 erstellt

# Example data
dna <- c("A", "T", "C", "T")

dnaset <- 
  lapply(seq(1, 10), function(x){

    set.seed(x)
    DNAStringSet(paste(sample(dna, 50, TRUE), collapse=""))

}) %>% Reduce(c, .)
names(dnaset) <- paste0("contig", 1:length(dnaset))

dnaset
#> DNAStringSet object of length 10:
#>      width seq                                              names               
#>  [1]    50 ATCATACCTTCCAAATTTTCACA...ATTTATCATCTTTTCTTTTTAC contig1
#>  [2]    50 ACTTTTAAATAATCACTTTTCCC...TTCTATTTATTTCACTAATTTC contig2
#>  [3]    50 ATTCTTTCCTTTCTTTTATATAA...TTTTATTTTCTTCTCACCTCAC contig3
#>  [4]    50 TCCCCTCTATCTTTTTTATCCCC...TTCTACATATTATCCTTATATT contig4
#>  [5]    50 TCACCAAATCCCCTTCATTTATC...ATTCATTTTACCCTTTTTATTT contig5
#>  [6]    50 ATATTTTCATTTTTCCATATCAT...TTCACTTTATTACAATACACTC contig6
#>  [7]    50 TCCTCTCTTTTCTTTCTCTTCTT...TCATCTTCTTTTCTTACTTACA contig7
#>  [8]    50 TTCTCTCTCTAATCCTTTCTTTA...CCTCTATCCTATATTTATCTTT contig8
#>  [9]    50 CATCTCCTTCTTATATTCTTATT...CTTTTCCTTAAATCCCACATCC contig9
#> [10]    50 CATTTCTTCCCTCCTCTTTACTT...TAAAATTATATACTTTATATTC contig10

# second contig has this pattern
pattern <- "TCACTTTTCCCTCAT"             

v <- vmatchPattern(pattern, dnaset)

# list entries that are non-empty, so that have a match
w <- lengths(v)>0
w
#>  [1] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

dnaset[w]
#> DNAStringSet object of length 1:
#>     width seq                                               names               
#> [1]    50 ACTTTTAAATAATCACTTTTCCC...ATTCTATTTATTTCACTAATTTC contig2

Created on 2023-07-27 with reprex v2.0.2