Get chromosome sizes from fasta file
4
6
Entering edit mode
5.9 years ago
rioualen ▴ 590

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 • 25k views
ADD COMMENT
1
Entering edit mode

A quick google search would help you. For e.x

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
21
Entering edit mode
5.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.

ADD COMMENT
0
Entering edit mode

Looks perfect! Thank you.

ADD REPLY
20
Entering edit mode
5.9 years ago
rioualen ▴ 590

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

samtools faidx input.fa
cut -f1,2 input.fa.fai > sizes.genome
ADD COMMENT
0
Entering edit mode

ALSo can be in this format at once:

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

ADD REPLY
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.

ADD REPLY
1
Entering edit mode
2.6 years ago
gbdias ▴ 130
  • 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

ADD COMMENT
0
Entering edit mode
3.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>
    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))
ADD COMMENT
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.

ADD REPLY

Login before adding your answer.

Traffic: 1747 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6