BWA index is not working
0
0
Entering edit mode
3.3 years ago

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.0k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

I see no error here.

what is the output of

ls -la ref_genome
ADD REPLY
0
Entering edit mode
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
ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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!

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 1289 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