Get chromosome sizes from fasta file
4
6
Entering edit mode
6.9 years ago
rioualen ▴ 660

Hello,

I'm wondering whether there is a program that could calculate chromosome sizes from any fasta file? The idea is to generate a tab file like the one expected in bedtools genomecov for example.

I know there's the fetchChromSize program from UCSC, but not all genomes are available over there (I need TAIR10 for instance). I've read this topic already.

I would like a tool that can deal with any genome regardless of the database. If it doesn't exist I guess it's possible to just parse fasta files, but I'd be surprised if no one else had done it before!

Cheers

genome ucsc • 32k views
1
Entering edit mode

http://www.danielecook.com/generate-fasta-sequence-lengths/

0
Entering edit mode

I did look it up on google, but thank you anyway.

24
Entering edit mode
6.9 years ago
pip install pyfaidx
faidx input.fasta -i chromsizes > sizes.genome


That's what you're looking for if you want to use bedtools genomecov, but you can transform to a BED file as well using -i bed.

0
Entering edit mode

Looks perfect! Thank you.

24
Entering edit mode
6.9 years ago
rioualen ▴ 660

Just found out faidx is available with samtools so another possibility is:

samtools faidx input.fa
cut -f1,2 input.fa.fai > sizes.genome

1
Entering edit mode

ALSo can be in this format at once:

samtools faidx genome.fa | cut -f1,2 > chromsizes

0
Entering edit mode

Tried this. Didn't work because samtools faidx genome.fa outputs to genome.fa.fai. Therefore, the command by rioualen should be the one used.

Unless maybe you use some type of flags in the command, maybe you could use pipes to run it as a one-line command as Apex92 tried.

0
Entering edit mode

This is the only line that worked for me, except the line only produced a output.fa.fai file which was identical to the chromsizes file i was looking except for extra columns (ie cut -f1,2 didnt work in this line)

samtools faidx input.fa | cut -f1,2 > chromsizes

However after getting that .fa.fai file i was able to cut out just the columns i needed (chromosome name & size) to produce the .sizes files i needed (also correction to the above the cut command needs to be given an input and output files).

cut -f1,2 Hel_final_2016_new.fa.fai > Hel_final_2016_new.fa.sizes

As a result you can run this line twice to make it work (since the first time it runs it will create the .fai file that the second use of the line will require):

samtools faidx Hel_final_2016_new.fa | cut -f1,2 Hel_final_2016_new.fa.fai > Hel_final_2016_new.fa.sizes

0
Entering edit mode

samtools faidx does not work for this purpose.

2
Entering edit mode
3.6 years ago
gbdias ▴ 150
• Another good option is to use the faSize command from UCSC tools. It works on any fasta file, not just UCSC genomes.

• You run it like this: faSize -detailed -tab file.fasta

• Default behavior is to print to STDOUT, so you can redirect the output to a file like this: faSize -detailed -tab file.fasta > output.txt

• The output will be a tab delimited file with sequence name on the first column and the length on the second.

• You can easily install faSize and other tools from UCSC using Bioconda https://anaconda.org/bioconda/ucsc-fasize

0
Entering edit mode
4.2 years ago
aleferna • 0

samtools was crashing because I included the HLA's,

HLA-A*01:01:38L 3374


Traceback (most recent call last):
File "/usr/local/bin/faidx", line 11, in <module>
File "/usr/local/lib/python2.7/dist-packages/pyfaidx/cli.py", line 197, in main
write_sequence(args)
File "/usr/local/lib/python2.7/dist-packages/pyfaidx/cli.py", line 50, in write_sequence
outfile.write(transform_sequence(args, fasta, name, start, end))
File "/usr/local/lib/python2.7/dist-packages/pyfaidx/cli.py", line 120, in transform_sequence
line_len = fasta.faidx.index[name].lenc
KeyError: 'HLA-A*01'


Finally just did:

from Bio import SeqIO
for rec in SeqIO.parse("hg38-Mix.fa","fasta"):
print rec.id+"\t"+str(len(rec.seq))

0
Entering edit mode

faidx was failing because the cli script expects UCSC-style regions encoded as contig:start-end, where start-end is optional. Since the HLA alt contains contain : characters you ran into a parsing issue. If you want to use pyfaidx in the same way as biopython you can do:

from pyfaidx import Fasta
for rec in Fasta("hg38-Mix.fa"):
printrec.name + str(len(rec)))


This will be much faster if you have a lot of sequences since none of the pyfaidx operations require reading the sequence unless you start slicing strings.