This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get chromosome sizes from fasta file

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
• 56,069 views • link
updated by Dave Carlson Dave Carlson • written by rioualen rioualen

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

• 0 views • link
written by rioualen rioualen

5 answers 

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 views • link
updated by Ram Ram • written by Matt Shirley Matt Shirley

Looks perfect! Thank you.

• 0 views • link
updated by Ram Ram • written by rioualen rioualen

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

samtools faidx input.fa
cut -f1,2 input.fa.fai > sizes.genome
• 0 views • link
updated by Ram Ram • written by rioualen rioualen

ALSo can be in this format at once:

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

• 0 views • link
written by Apex92 • Apex92

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 views • link
written by Pratik Pratik

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 views • link
written by jamesogilvie1 • jamesogilvie1

the pipe in this case is not good practice here. there is no piping being used. replace with a semicolon, or "&&"

• 0 views • link
written by cmdcolin cmdcolin

samtools faidx does not work for this purpose.

• 0 views • link
written by tomas4482 • tomas4482

  • 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 views • link
written by gbdias gbdias

Hi all,

check this: https://github.com/alejandrogzi/chromsize

Usage

Binary

Usage: chromsize --fasta <FASTA> --output <OUTPUT> [-t <THREADS>]

Arguments:
    -f, --fasta <FASTA>: FASTA file
    -o, --output <OUTPUT>: path to chrom.sizes

Options:
    -t, --threads <THREADS>: number of threads [default: your max ncpus]
    --help: print help
    --version: print version

build in rust, the simplest and fastest way to get them. Works as a binary, library and is ported to python. Here is the benchmark with a lot of tools:

benchmark

• 0 views • link
written by alejandrogzi • alejandrogzi

I guessed before checking that this tool was written in rust. I was right. :)

• 0 views • link
written by Dave Carlson Dave Carlson

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>
    load_entry_point('pyfaidx==0.5.2', 'console_scripts', 'faidx')()
  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))
• 1 views • link
updated by Ram Ram • written by aleferna • aleferna

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.

• 1 views • link
written by Matt Shirley Matt Shirley

Log in to answer this question.

More posts like this