This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Bwa Index On All Human Grch37 Sequences

Hi,

Is there any reason why "bwa index -a bwtsw" would fail on all the current toplevel sequences in the human GRCh37 assembly, including the patches?

wget ftp://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.61.dna.toplevel.fa.gz

~/bwa-0.5.9rc1/bwa index -a bwtsw Homo_sapiens.GRCh37.61.dna.toplevel.fa.gz

After a good while, it fails with: BWTIncConstructFromPacked: Cannot determine file length!

Cheers

bwa human
• 17,602 views • link
updated by Pablo Marin-Garcia Pablo Marin-Garcia • written by 2184687-1231-83- 2184687-1231-83-

@Pierre: yes, it does.

• 0 views • link
written by lh3 lh3

By the way, following lh3 recomendatations (in this BioStar question)

I will use the 1kgenomes reference (and it is already indexed), that does not have the haplotype regions (prevent some variation calls)

wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.fai
• 0 views • link
updated by Ram Ram • written by Pablo Marin-Garcia Pablo Marin-Garcia

No, you should NOT include allelic sequences. You lose calls rather than gain them. You need special treatment to handle MHC and the chr17 inversion.

• 0 views • link
written by lh3 lh3

2 answers

Bwa accepts gzipped fasta/fastq files in indexing and alignment. It automatically tells compression and whether the input is fasta or fastq. The problem is most likely to be caused by the use of the "toplevel" fasta, which is longer than 4GB concatenated. Please use the reference genome from the 1000 genomes project.

• 1 views • link
written by lh3 lh3

@Albert, I don't know if this would be a problem with bwa, but we aware that ensembl fasta files had two Y chromosome fasta entries with the same name 'Y':

$  grep '^>Y' Homo_sapiens.GRCh37.62.dna.toplevel.fa
 >Y dna:chromosome chromosome:GRCh37:Y:1:10000:1
 >Y dna:chromosome chromosome:GRCh37:Y:2649521:59034049:1

And many programs don't like that at all.

• 1 views • link
updated by Ram Ram • written by Pablo Marin-Garcia Pablo Marin-Garcia

Log in to answer this question.

More posts like this