How to create a .bed file from .fasta? [duplicate]

Question

I have some problems in creating a .bed file for hg19, so I will be able to visualize the .bed file in IGV.

The .fasta file contains rows of this form:

>chr1:0-1000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:1000-2000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:2000-3000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:3000-4000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

I tried both with biopython:

from Bio import SeqIO

with open(f1) as in_f, open(f1+'.seq','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record))

and awk:

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

Both ways gave me the same output which IGV couldn't work with:

chr1:0-1000 0   1000
chr1:1000-2000  0   1000
chr1:2000-3000  0   1000
chr1:3000-4000  0   1000
chr1:4000-5000  0   1000
chr1:5000-6000  0   1000
chr1:6000-7000  0   1000
chr1:7000-8000  0   1000

Now I tried to see what would be a valid .bed file format from UCSC and the format looks different. For example:

#bin    chrom   chromStart  chromEnd
13  chr1    39845715    39846245
19  chr1    87030756    87031903
22  chr1    111148649   111149175
22  chr1    113246065   113246628

1. Why my methods don't output the .bed in the same format?
2. How can I make the awk/biopython methods to output a .bed file that IGV can actually load.

Answer 1

Here is a fun little python script I cooked up for the occasion. It will take a standard fasta file as a command-line argument and turn it into the proper bed format. Using your example fasta:

import pandas as pd
import sys

inFasta = sys.argv[1] # take fasta as command argument

def fastaParser(fasta):
    headers = []
    with open(fasta) as f:
        header = None
        for line in f:
            if line.startswith('>'): # identifies fasta header line
                headers.append(line[1:-1]) # append all of the line that isnt >
                header = line[1:] # in reset header
    newHeader = (header.replace(':',',') for header in headers) # format to be accepted later
    newnewHeader = (header.replace('-',',') for header in newHeader) # format to accept later
    bedHead = (header.split(',') for header in newnewHeader) # separate by comma from format above
    return bedHead

fastaParsed = fastaParser(inFasta) # run the function

headers = list(fastaParsed) # go from generator to list
bedfile = pd.DataFrame(headers) # create dataframe which will output bedfile

bedfile.to_csv('fastabedfile.bed', sep='\t', index=False, header=None) # output bedfile format

to run:

[dkennetz@nodecn202  fasta_parser]$python fastaParser.py fakeFasta.fa

[dkennetz@nodecn202  fasta_parser]$more fastabedfile.bed
chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

Answer 2

I'm assuming the FASTA header contains all the information you want. You can easily generate a proper bed file with Biopython:

from Bio import SeqIO

for record in SeqIO.parse('test.fa', 'fasta'):
    lst = record.id.split(':')
    chr_, (start, end) = lst[0], lst[1].split('-')
    print('\t'.join((chr_, start, end)))

Note that native Python code is only marginally longer:

with open('test.fa') as f:
    for line in f:
        if line.startswith('>'):
            lst = line.rstrip().split(':')
            chr_ = lst[0].lstrip('>')
            start, end = lst[1].split('-')
            print('\t'.join((chr_, start, end)))

Both code snippets output bed format:

chr1    0   1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

Answer 3

If your headers are always of the form >chrZ:a-b, where chrZ is a chromosome name, a is a zero-indexed, half-open start position, and b is a zero-indexed, half-open stop position:

$ awk '($0 ~ /^>/){ print substr($0,2) }' in.fa | sed 's/[:|-]/\t/g' > out.bed

If your headers contain coordinates using a different index, then you can make adjustments to $1 and $2 in the awk statement to create zero-indexed, half-open coordinates to export to BED, depending on the initial index.

Answer 4

The UCSC format you linked to isn't a BED file, your method should never produce it. What you posted as your desired output is also not a BED file. Below is a BED file:

chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

awk and biopython were producing the appropriate thing already. You just need to zoom in on a region in IGV. Since it's not a bigBed file, IGV won't show everything at all zoom levels.

As a demonstration, below is the BED file from above in IGV:

Answer 5

Assuming that you have a reference genome sequence file (e.g., reference.fasta) available. Getting this kind of file is straightforward.

1) Index the reference genome and map your reads or FASTA sequences to it (for example with bowtie2)

# index reference genome (should be precomputed)
bowtie2-build reference.fasta reference
# map reads
bowtie2 -x  reference -U file.fasta -S file.sam
# compress SAM to a BAM (binary) file
samtools view -bS file.sam > file.bam

2) Create a BED files from the BAM file using bedtools

# extract BED file
bedtools bamtobed -i file.bam > file.bed

Alternatively, for sliding windows you can generate these from a reference sequence provided that you know the length of each chromosome (perhaps there is a way to extract these directly from reference.fasta):

# length per chromosome
samtools view -H file.bam | grep "SQ" | cut -d":" -f2-3 | sed 's/LN://' > file.chr.txt

Create sliding windows (using only the lengths per chromosome)

# create sliding windows
bedtools makewindows -g file.chr.txt -w 1000 > windows.bed