lowercase fasta coordinates to GFF

Question

[Biomonika (Noolean)]

I have fasta file which has been soft-masked to lowercase. Now, I would like to get GFF (BED) file with coordinates of all soft-masked sequences. I was wondering if such a tool is already available. Thank you.

Answer 1

[Pierre Lindenbaum]

Using gnu flex? (not fully tested, should work)

[Biomonika (Noolean)]

Thanks for an answer! I haven't used flex before, could you please explain in greater detail how to install/run this code? (flex: can't open jeter.l)

[Pierre Lindenbaum]

in the current context GNU flex runs like a awk for C.

for example

^>.*\n

is a pattern for the a fasta header.

[Biomonika (Noolean)]

Perfect! Thank you very much. I will just add note for the others that code reports are 0-based.

[Pierre Lindenbaum]

sorry I renamed the file when copy+paste, that should be flex code.l not flex jeter.l

[ramvizard]

Hi,

I'm trying to do the same but while testing the above script, I figured that I'm not getting the coordinates of the softmasked region at the end of the sequence.

For example, testing the script on the following three cases,

>QWEQ
qweqertADEADASDAWDW
>ASAD
ASQWEQQqweqeASAasa
>adad
asasasdaa

I get back,

QWEQ    0       7
ASAD    7       12
ASAD    15      18

But how would I change the flex script to give back the coordinates of the sequence "adad"?

Thanks,

Cheers,

[Pierre Lindenbaum]

ah yes, I forgot the trailing lowercase letters. I've added a <<EOF>> condition.

[Biostar]

This is very useful, but how do I run it with a locally stored .fasta file? I tried the following and I got an empty co-ordinate file and it doesn't finish running either. Also, is there a way to use multiple threads to run this?

flex code.l && \

gcc lex.yy.c && \

./a.out genome_softmasked.fasta > co-ords.gff

Answer 2

[bcontreras]

You can do this in Perl in two steps:

1) convert your FASTA file to 1 sequence per line:

perl -lne 'if(/^(>.*)/){ $head=$1 } else { $fa{$head} .= $_ } END{ foreach $s (sort(keys(%fa))){ print "$s\n$fa{$s}\n" }}'  infile.fasta > infile.1line.fasta

2) extract lower-case segments from that FASTA file:

perl -lne 'if(/^>(\S+)/){ $n=$1} else { while(/([a-z]+)/g){ printf("%s\t%d\t%d\n",$n,pos($_)-length($1),pos($_)) } }' infile.1line.fasta > infile.bed

In case you only want sements of some min length, 60 in the example:

perl -lne 'if(/^>(\S+)/){ $n=$1} else { while(/([a-z]{60,})/g){ printf("%s\t%d\t%d\n",$n,pos($_)-length($1),pos($_)) } }' infile.1line.fasta > infile.60.bed