Extracting 22-mers ending in GG from FASTA using sed/awk?

[Genomics101]

I feel like there must be a simple one-liner with see/awk that can do this, but can't think of it and I hope y'all can help be out.

I have a set of FASTA files for gene sequences, from GenBank. I need to extract all the instances of "GG," plus the twenty bases up stream. The output would be a list of 22-mers all ending in GG (with line breaks after each), all the instances of this in each gene.

Any help is greatly appreciated!

[maasha]

Let me see if I understand this correctly. You can use Biopieces (www.biopieces.org) like this (or modify a bit):

read_fasta -i genes.fna | # read in your sequences
split_seq -w 22 | # split the sequences in 1-base overlapping 22-mers
grab -r 'GG$' -k SEQ | # grab all k-mers ending in GG
write_fasta -xo 22-mers_terminal_GG.fna

[RickBioinf]

I think he means that you first find the 'GG' and then take the 20 bases before the 'GG' and the 'GG' itself, so not splitting the sequences in 22-mers before hand.

[JamieHeather]

Do you expect there to be more than one instance per line?

[JamieHeather]

Wait, never mind, remembered the flag needed:

grep -io ......................GG input.fasta > output

(-i = case insensitive, -o = just outputs match)

Edit - no wait, this won't output overlapping matches, never mind!

You might have to go a little higher level than bash tools.

[dariober]

Originally Posted by Genomics101

I have a set of FASTA files for gene sequences, from GenBank. I need to extract all the instances of "GG," plus the twenty bases up stream. The output would be a list of 22-mers all ending in GG (with line breaks after each), all the instances of this in each gene.

Not exactly a one liner but it might do the trick...

## Test fasta file
cat test.fa 
>seq1
GGGGAATTAGCTCAAGCGGTAGAGCGCTCCCTTAGCATGCGAGAGGTAGCGGGATCGACG
CCCCCATTCTCTA
>seq2
GGGGGATTAGCTCAAGCGGTAGGGTGCCTGCTTAGCATGCAAGAGGTAGCAGGATCGACG
CCTGCATTCTCCA

python -c "
import re
fin= open('test.fa')
seqname= fin.readline().strip().lstrip('>')
faseq= []
while True:
    line= fin.readline().strip()
    if line.startswith('>') or line == '':
        mpat= re.finditer(r'(?=(.{20}GG))', ''.join(faseq))
        for m in mpat:
            print(seqname + '\t' + m.group(1))
        seqname= line.lstrip('>')
        faseq= []
    else:
        faseq.append(line)
    if line == '':
        break
"
## Output:
seq1	CGCTCCCTTAGCATGCGAGAGG
seq1	CTTAGCATGCGAGAGGTAGCGG
seq1	TTAGCATGCGAGAGGTAGCGGG
seq2	GGGGATTAGCTCAAGCGGTAGG
seq2	GGGATTAGCTCAAGCGGTAGGG
seq2	TGCCTGCTTAGCATGCAAGAGG
seq2	TTAGCATGCAAGAGGTAGCAGG

Dario

[maasha]

OK, another try with Biopieces - this time from GenBank files ...

read_genbank -i test.gb -f gene | # read in genbank entries with gene feature key
add_ident -k SEQ_NAME | # add an identifier
patscan_seq -cp '20 ... 20 GG' | # scan for pattern
write_tab -xck S_ID,MATCH,S_BEG,S_END,STRAND # output table