Member Avatar for sudipta.mml

The quadruplex sequence of a genome looks like this Gx Ny1 Gx Ny2 Gx Ny3 Gx, where G is the Guanine base and the Ns are representing other bases. The x, y1, y2 and y3 are integer. A particular segment will be quadruplex sequence if x>=2. My question is I want to count how many such segments are present in a FASTA sequence and what are their types in terms of x, y1,y2 and y3. For that purpose I need some help from your side.

  • 4 Contributors
  • 4 Replies
  • 364 Views
  • 4 Years Discussion Span
  • Latest Post Latest Post by Matej_1

Dani AI

Generated

This thread already has the right pointers: use a FASTA-aware parser and capture match positions rather than scanning the raw file. The practical pitfalls are (1) FASTA headers/newlines — they shift file offsets vs. sequence coordinates, (2) overlapping occurrences, and (3) enforcing the quadruplex definition GxNy1GxNy2GxNy3Gx with the same x for all four G-runs. pointed you toward Biopython for parsing and showed how to get positions; below is a compact, practical pattern and workflow that ties those together and also enforces equal-length G-runs.

Read each record with SeqIO so the sequence string has no headers/newlines, compile a lookahead+backreference regex to:

  • capture the first G-run (G{2,}),
  • require the same run to repeat using backreferences,
  • capture the three loop lengths, and
  • allow overlapping matches (lookahead).
    Make the maximum allowed loop length a parameter so you can tune sensitivity.

Example snippet (adjust max_loop and allowed loop chars as needed):

from Bio import SeqIO
import re
from collections import Counter

max_loop = 12
pat = re.compile(r'(?=(G{2,})([ACGTN]{1,%d}?)(\1)([ACGTN]{1,%d}?)(\1)([ACGTN]{1,%d}?)(\1))' % (max_loop, max_loop, max_loop), re.I)

counts = Counter()
for rec in SeqIO.parse('input.fasta', 'fasta'):
    seq = str(rec.seq).upper()
    for m in pat.finditer(seq):
        x = len(m.group(1)); y1 = len(m.group(2)); y2 = len(m.group(4)); y3 = len(m.group(6))
        start = m.start(1) + 1          # 1-based inclusive
        end = m.end(7)                  # 1-based inclusive end
        counts[(x,y1,y2,y3)] += 1
        print(rec.id, start, end, x, y1, y2, y3)

Notes and troubleshooting

  • If you must work from raw file offsets (not SeqIO), build an index that maps file positions to sequence coordinates before searching.
  • For very large genomes, compile the pattern once and consider scanning in windows to limit memory/CPU.
  • Cross-check results with specialized tools such as pqsfinder (mentioned by ) if you need scoring or more biological filters.

This approach keeps coordinates correct per-sequence, reports each motif’s x,y1,y2,y3 type, and avoids missing overlapping quadruplexes.

Recommended Answers

All 4 Replies

Member Avatar for dashing.adamhughes

You may want to start using biopython as it has many classes and methods for handling sequence data, especially from a fasta-formatted file. It may already do some of the work for you.

Member Avatar for sudipta.mml

import re
fasta = open('e-coli-k12.fasta', 'r').read()
segments=re.findall('GG+[^G]+',fasta)
print segments

This script produces a pattern which starts with more than two G letters for a sequence. But I can't extract the start and end position of that pattern.

I have tried with

import re
fasta = open('e-coli-k12.fasta', 'r').read()
segments=re.compile('GG+[^G]+')
for item in segments.findall(fasta):
print item
found=re.search(item,fasta)
print found.span()

but I didn't get success. It searches a particular pattern over the whole sequence and produces multiple number for a particular pattern. But I want the only one exact start and end position which exactly corresponds to the fasta sequence. How to get the start and end position of a particular pattern.

Member Avatar for snippsat

As posted by there are made a lot of stuff that can help you with this,another eksp is pyfasta

How to get the start and end position of a particular pattern

You can use re.finditer for this.

>>> import re
>>> s = '11ATGC1111ATGC11111ATGC'
>>> p = 'ATGC'
>>> [m.start() for m in re.finditer(p, s)]
[2, 10, 19]
>>> #To also find end postion
>>> [(m.start(),m.end()) for m in re.finditer(p, s)]
[(2, 6), (10, 14), (19, 23)]
Edited by snippsat
Member Avatar for Matej_1

If you don't mind calling R from python (or vice versa), the pqsfinder package in R (http://bioconductor.org/packages/pqsfinder/) solves most of the quadruplex sequence search and manipulation. For example, it already has functions to retrieve loop lengths, positions etc.

commented: Thanks for sharing. +14
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.