This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BWA index is not working

I am trying to align SRA to a mitochondrial reference genome that has been pulled in fasta form from NCBI. The output I keep getting is not normal, its going through it in 0.023 seconds.

[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 0.00 seconds elapse.
[bwa_index] Update BWT... 0.00 sec
[bwa_index] Pack forward-only FASTA... 0.00 sec
[bwa_index] Construct SA from BWT and Occ... 0.00 sec
[main] Version: 0.7.5a-r405
[main] CMD: bwa index -a is ref_genome/sequence.fasta
[main] Real time: 0.009 sec; CPU: 0.023 sec

then nothing else works beyond this. I will get an error about not being able to retrieve the indexed file..

I dont know if this is a good enough explanation and I am definitely new to this, so thank you for your patience.

I can follow up with any other info if necessary.

alignment genome
• 2,991 views • link
written by samantharenee98 • samantharenee98

Hello, please show the code you run, without that and the error message this is not reproducible.

• 0 views • link
written by ATpoint ATpoint

So sorry.

bwa index ref_genome/sequence.fasta and bwa index sequence.fasta

and ive also done 

bwa index -p ref_genome/sequence.fasta -a bwtsw ref_genome/sequence.fasta

Output: 

[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=33052, availableWord=65536

and 

bwa index -a is ref_genome/sequence.fasta

Output: 

[bwa_index] Pack FASTA... 0.00 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 0.00 seconds elapse.
[bwa_index] Update BWT... 0.00 sec
[bwa_index] Pack forward-only FASTA... 0.00 sec
[bwa_index] Construct SA from BWT and Occ... 0.00 sec
[main] Version: 0.7.5a-r405
[main] CMD: bwa index -a is ref_genome/sequence.fasta
[main] Real time: 0.009 sec; CPU: 0.023 sec
• 0 views • link
updated by Istvan Albert Istvan Albert • written by samantharenee98 • samantharenee98

I see no error here.

what is the output of 

ls -la ref_genome
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum 
total 144
drwxr-xr-x@ 8  staff    256 Jan 14 18:19 ./
drwxr-xr-x@ 8  staff    256 Jan 14 18:19 ../
-rw-r--r--@ 1   staff  16830 Jan 14 13:02 sequence.fasta
-rw-r--r--  1  staff     10 Jan 14 18:19 sequence.fasta.amb
-rw-r--r--  1  staff     88 Jan 14 18:19 sequence.fasta.ann
-rw-r--r--  1   staff  16624 Jan 14 17:49 sequence.fasta.bwt
-rw-r--r--  1   staff   8265 Jan 14 18:19 sequence.fasta.pac
-rw-r--r--  1   staff   8312 Jan 14 17:49 sequence.fasta.sa
• 0 views • link
updated by i.sudbery i.sudbery • written by samantharenee98 • samantharenee98

Yes, that is the normal indexing output, but you say:

I will get an error about not being able to retrieve the indexed file..

so there must be more code that you run, so the actual alignment?

• 0 views • link
written by ATpoint ATpoint

This is what ive done after that indexing. 

$ bwa mem data/ref_genome/sequence.fasta data/fundhe.fastq > results/sam/fundhe.sam

[main] Version: 0.7.5a-r405
[main] CMD: bwa mem data/ref_genome/sequence.fasta data/fundhe.fastq
[main] Real time: 0.004 sec; CPU: 0.021 sec

$ head -n 10 fundhe.sam

@SQ SN:NC_012312.1  LN:16526

$ samtools view -bSo fundhe.bam fundhe.sam

No output

$ samtools flagstat fundhe.bam

Output:

0 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 mapped (N/A : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Also just want to thank you all for you fast responses!

• 0 views • link
updated by i.sudbery i.sudbery • written by samantharenee98 • samantharenee98

First things I'd check are that the formats of your fastqs and fastqs are okay.

• 0 views • link
written by swbarnes2 swbarnes2

Odd, what is the output of head data/fundhe.fastq and wc -l data/fundhe.fastq? Looks like either the files are empty or corrupted because the first line of flagstat would also show unmapped reads and I guess if the index was corrupted bwa would complain as well, so it is probably the fastq being off as @swbarnes2 says.

• 0 views • link
written by ATpoint ATpoint

0 answers

No answers yet.

Log in to answer this question.

More posts like this