14
$\begingroup$

I have a FASTA file:

> Sequence_1
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC...
> Sequence_2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC....
....

I want to generate a BED file for each sequence like:

Sequence_1 0 1500
Sequence_2 0 1700

The BED regions will simply be the size of the sequences.

Q: I did that before with a one-line command. I don't remember what that is, it was on Biostars. I can't find the post now. What's the simplest way to do the conversion?

$\endgroup$
  • $\begingroup$ Do you mean this biostars post? biostars.org/p/15476 $\endgroup$ – Greg May 18 '17 at 0:42
  • $\begingroup$ @Greg The post doesn't give me one-line command to do what I want. I'm 100% sure I saw it somewhere on biostars a while ago, but I don't remember the name of the program that can do it. $\endgroup$ – SmallChess May 18 '17 at 0:43
  • $\begingroup$ What about this one? biostars.org/p/15476/#189089 $\endgroup$ – Greg May 18 '17 at 0:49
  • $\begingroup$ @Greg It was a Python program. I don't remember the name and thus no idea how to search. $\endgroup$ – SmallChess May 18 '17 at 0:50
  • 5
    $\begingroup$ Would we need to create a new tag, maybe something like "format-conversion", "file-conversion"? As I expect there will be many more convert from "X to Y format" type of questions. $\endgroup$ – zx8754 May 18 '17 at 7:13
12
$\begingroup$

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:

bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

-c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.

$\endgroup$
6
$\begingroup$

You could adapt this awk one-liner. Note that it assumes that sequence IDs are not longer than 100 characters and that there is no description following the sequence ID on the header line.

cat myseqs.fasta | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t0\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

Otherwise, any Bio* library (Perl, Python, Ruby) provides FASTA format parsers which will extract sequence IDs and lengths.

I'd point out that whilst this resembles BED it is not, strictly-speaking, since BED maps to coordinates on a chromosome or some longer sequence object.

$\endgroup$
6
$\begingroup$

It's good practice to have your FASTA indexed, so you can leverage the .fai you are likely to already have. If not, you can just generate the index with samtools and use some awk to make your BED:

samtools faidx $fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed

This will maintain tab separation but you can drop the BEGIN statement to use spaces. The BED spec only requires "whitespace" for the simple BED format.

$\endgroup$
4
$\begingroup$

Here is an approach with BioPython. The with statement ensures both the input and output file handles are closed and a lazy approach is taken so that only a single fasta record is held in memory at a time, rather than reading the whole file into memory, which is a bad idea for large input files. The solution makes no assumptions about the sequence ID lengths or the number of lines that the sequences are spread across:

from Bio import SeqIO

with open('sequences.fasta') as in_f, open('sequences.bed','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record)))
$\endgroup$
4
$\begingroup$

Inspired by this answer to a related question on read length distributions, you could do this with Biopython:

from Bio.SeqIO import parse
with open("regions.bed", "w") as bed:
    for record in parse("regions.fasta", "fasta"):
        print(record.id, 0, len(record.seq), sep="\t", file=bed)
$\endgroup$
  • 1
    $\begingroup$ I love Python. I prefer this answer simply because it is in Python and that is what we should be using in 2017. But why does this look so much like Perl? Readability counts. I wanted to propose a cleaner snippet, but it looks like comments do not allow code. Would it be considered good tone to edit your answer? $\endgroup$ – Kirill G May 18 '17 at 2:20
  • 1
    $\begingroup$ @SamStudio8 On the surface level, I agree. At the same time, if a person has a Python REPL constantly open (which is a lot of people), this has minimal overhead, and, in my humble opinion, quick algorithmic thinking should be promoted at all times in this profession. Moreover, if this conversion is part of a Snakemake pipeline, then it may be even more natural than calling shell() or writing a shell block. $\endgroup$ – Kirill G May 18 '17 at 12:20
  • 1
    $\begingroup$ @KirillG I usually have a REPL open, but I think I'd still favour using awk (or such). Years ago was a different story, but learning how to use the basics of awk has been one of the most productivity improving things I have done since switching to use Linux. But sure, I definitely agree if this is part of a larger script or pipeline, I wouldn't dare call shell. Though I'd probably use pysam over biopython, but that's just personal preference :) $\endgroup$ – Sam Nicholls May 18 '17 at 12:23
  • 2
    $\begingroup$ @KirillG I wrote a BioPython solution that is arguably more readable/pythonic and doesn't put the whole input file into memory bioinformatics.stackexchange.com/a/106/104 $\endgroup$ – Chris_Rands May 18 '17 at 13:08
  • 1
    $\begingroup$ @Chris_Rands' answer is much better than mine, please upvote it! $\endgroup$ – Scott Gigante May 18 '17 at 13:52
3
$\begingroup$

We have many excellent answers! This will be an excellent reference for future users.

I found what exactly what I was asking in my question:

https://www.biostars.org/p/191052/

$ pip install pyfaidx  
$ faidx --transform bed test.fasta > test.bed

This is the one-line command I was asking. The other answers also work, but I want to accept my own answer.

$\endgroup$
  • 1
    $\begingroup$ Developer of pyfaidx here. I was just about to write this answer and lament that I didn't pick a more pronounceable or memorable package name. $\endgroup$ – Matt Shirley May 22 '17 at 14:15
0
$\begingroup$

If the FASTA sequences are all on a single line, then the following perl one-liner should work:

cat myseqs.fasta | perl -ne 'if(/^>([^ ]+)/){print $1} else {print " 0 ",length,"\n"}'

Explanation:

  • If the line begins with a '>', then print everything up to the first space (but don't put a line break at the end)
  • Otherwise, print " 0 ", followed by the length of the line, followed by a line break
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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