3
$\begingroup$

I am currently working on a bioinformatics problem where I need to lookup and count the location and count of occurences of 4000-ish 5 character long patterns in each sequence of a fasta file of 700GB. This fasta file has aligned sequences of one virus. To start with, I want to just look up 6 patterns.

The characters are limited to ATGC for the patterns and the fasta.

The input to the code is just the patterns I want to look up and the aligned sequences of the virus.

The output should be a file which is of a matrix format with #rows = num of patterns, # columns = #num of characters in any given fasta sequence (should be same for all sequences in the fasta as it is an aligned file). Each value of the matrix should be the % of sequences that have that pattern starting in that coordinate.

For example-
input:

seq1 - ATGCCAT
seq2 - ATG-GAT
seq3 - GGAATTC
seq4 - -ATATCC

pattern #1 - AT
pattern #2 - GC

NUMBER OF SEQUENCES = 4
'-' = any char which has not been sequenced/ not in ATGC total number of sequences at any particular coordinate should be = number of seq - num of invalid sequences at that coordinate

output matrix:

Pattern 1 2 3 4 5 6
AT 2/3 1/4 0 2/3 0 2/4
GC 0 0 1/3 0 0 0

let me explain the value 2/3 in position 1*1 of the matrix - in the first coordinate of 2 out of the 4 sequences, the pattern AT starts. 1 sequence has an invalid - character and hence will be excluded from the total number of considered sequences.

I am a newbie to bioinformatics applications and hence am unaare if this is common problem but I am yet to discover a straghtforard solution through my research so far. I recently came across the aho-corasick algortihm for multi-pattern matching with strings through large text files as a very memory efficient and quick means of solving such issues. I am interested in scaling this search operation to many more patterns in the future as well. I am looking for advice on the following:

  1. Are there any well- established tool for such and operation already? if so, what are they?
  2. Am I right in attempting to use aho-corasick for this application?
  3. How can I efficiently apply this algorithm to a largned fasta with aligned sequenced? can I read the fasta in chunks to make it quicker?
  4. Should I code this in something other than Python? (I have been using Python so far but I am open to others as well).

Open to any suggestions as I am new to big data processing as well

$\endgroup$
4
  • 2
    $\begingroup$ This question needs more context about the overall project, because it seems like an odd task to carry out. What is the purpose of doing this? What will you do with these match results once the answer has been found? $\endgroup$
    – gringer
    Commented May 30 at 20:32
  • 1
    $\begingroup$ Wait ... 700GB ... I have not read the question in detail, but I have produced solutions involving industrial data engineering (I work in industry), but you'd need to be an expert and it's an entire project. Thus it is a full scale parallelisation OR a single CPU with massive RAM (AWS stuff). I'll think if there's an easy way in ... BTW is 700GB zipped or unzipped (most fasta/fastq are zipped))? $\endgroup$
    – M__
    Commented May 30 at 21:04
  • $\begingroup$ The unzipped ones is 700GB $\endgroup$ Commented Jun 7 at 3:39
  • $\begingroup$ The purpose is to calculate conservation scores $\endgroup$ Commented Jun 7 at 3:39

2 Answers 2

1
$\begingroup$

I wonder if FIMO might be something that you could use. The online version doesn't have the capacity 70 deal with 700GB of data, but you could try it out with a subsample to see if the output is something that would fit your needs. https://meme-suite.org/meme/tools/fimo The MEME suite (including FIMO) is also available via conda, or Docker etc.

$\endgroup$
1
$\begingroup$

Maybe a job for seqkit locate. Here -m 0means no mismatches allowed. I do not understand really your output format, so maybe start with seqkit and see whether you can tune it for your problem. It allows mismatches by some parameters.

$ cat foo.fa 
>entry1
TACTACGATCACTGACGTAGCACATTGGATTTAGAGTACGATCGAT
>entry2
TAGCACTATCGACGTAGCTACGTACGTACTGATCGATCTGACACGA
>entry3
TGGTGTGTTGGTGTGTTGAAAAAAAAAAAAAAAAAAAAAAAAAAAA

$ cat pattern.fa 
>pattern1
TAGCAC

$ ./seqkit locate -m 0 -f pattern.fa foo.fa 
seqID   patternName pattern strand  start   end matched
entry1  pattern1    TAGCAC  +   18  23  TAGCAC
entry2  pattern1    TAGCAC  +   1   6   TAGCAC

https://bioinf.shenwei.me/seqkit/usage/#locate

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.