I believe the XA information is wrong in many of our reads but I can't understand why. Any help will be greatly appreciated. Essentially I want to be able to rely on the XA information (I'm not only concerned about uniq reads).
History:
- these are single-end reads from SOLiD
- bwa flags are as follows:
bwa aln -c -n 0.04 -N -q 30 -t 60 -f ${OUTFILE}.sai $INDEX ${FASTQ}
bwa samse -n 1000000 -f ${OUTFILE}.txt $INDEX ${OUTFILE}.sai ${FASTQ}
- input file that is indexed has 1314 sequences and is a total of 2.3G
- i'm getting lots of problematic reads (> 74,000) over 16 samples. Could be much more than 74,000 but those were easy to verify
- one of the problem reads are below with questions/concerns I'm hoping to be able to get assistance with
A sample problem read:
208_44_1461 0 gi|307149945|ref|NC_014501.1|_Cyanothece_sp._PCC_7822_chromosome__complete_genome 4159110 25 20M * 0 0 AAATACCGTTCAAATAATGG A!!57A?:@CHQH%=O79:< XT:A:U CM:i:0 X0:i:1 X1:i:1 XM:i:2 XO:i:0 XG:i:0 MD:Z:20 XA:Z:gi|258513327|ref|NC_013213.1|_Acetobacter_pasteurianus_IFO_3283-01_plasmid_pAPA01-040__complete_sequence,+3199,7M1I1M2I4M1D5M,2;
Some concerns with the sample read above:
1. Cyanothece seems to map, but Acetobacter does not. When looking at the Acetobacter sequence, I don't see anything that even closely resembles AAATACCGTTCAAATAATGG. Not even when taking into account the 3 insertions and 1 deletion mentioned in the CIGAR score. As you will see from the X1 score, Acetobacter is a reported suboptimal hit.
2. The Acetobacter has a start coordinate of +3199 (which makes the coordinates 3199 - 3218). However, the actual sequence of Acetobacter is only 3205 in length, so it can't go to 3218. I don't believe it's a problem with boundries because even when looking at 3199-3205 (that are in the sequence), nothing resembles AAATACCGTTCAAATAATGG. I also looked at the sequences surrounding that sequence and didn't see anything it could map to (even taking the reverse complement). My bigger fear is there are also reported reads that don't exceed the sequence boundry that don't map at all.
3. Acetobacter has a CIGAR score of 7M1I1M2I4M1D5M which I believe should have been 'filtered out' given the bwa flags shown above. While I still don't see the alignment, the CIGAR score likely should have been filtered despite it being a suboptimal hit (unless I'm not understanding the suboptimal hit thresholds correctly).
Any help will be greatly appreciated. I'm mainly trying to understand if we should continue to use BWA for reading 'non-uniq' reads, use some other method, or apply some work-arounds. Could the problems be only with sub-optimal hits? Or could it be more wide-spread and effect all repeat-hits or even uniq hits?
Thank you.
History:
- these are single-end reads from SOLiD
- bwa flags are as follows:
bwa aln -c -n 0.04 -N -q 30 -t 60 -f ${OUTFILE}.sai $INDEX ${FASTQ}
bwa samse -n 1000000 -f ${OUTFILE}.txt $INDEX ${OUTFILE}.sai ${FASTQ}
- input file that is indexed has 1314 sequences and is a total of 2.3G
- i'm getting lots of problematic reads (> 74,000) over 16 samples. Could be much more than 74,000 but those were easy to verify
- one of the problem reads are below with questions/concerns I'm hoping to be able to get assistance with
A sample problem read:
208_44_1461 0 gi|307149945|ref|NC_014501.1|_Cyanothece_sp._PCC_7822_chromosome__complete_genome 4159110 25 20M * 0 0 AAATACCGTTCAAATAATGG A!!57A?:@CHQH%=O79:< XT:A:U CM:i:0 X0:i:1 X1:i:1 XM:i:2 XO:i:0 XG:i:0 MD:Z:20 XA:Z:gi|258513327|ref|NC_013213.1|_Acetobacter_pasteurianus_IFO_3283-01_plasmid_pAPA01-040__complete_sequence,+3199,7M1I1M2I4M1D5M,2;
Some concerns with the sample read above:
1. Cyanothece seems to map, but Acetobacter does not. When looking at the Acetobacter sequence, I don't see anything that even closely resembles AAATACCGTTCAAATAATGG. Not even when taking into account the 3 insertions and 1 deletion mentioned in the CIGAR score. As you will see from the X1 score, Acetobacter is a reported suboptimal hit.
2. The Acetobacter has a start coordinate of +3199 (which makes the coordinates 3199 - 3218). However, the actual sequence of Acetobacter is only 3205 in length, so it can't go to 3218. I don't believe it's a problem with boundries because even when looking at 3199-3205 (that are in the sequence), nothing resembles AAATACCGTTCAAATAATGG. I also looked at the sequences surrounding that sequence and didn't see anything it could map to (even taking the reverse complement). My bigger fear is there are also reported reads that don't exceed the sequence boundry that don't map at all.
3. Acetobacter has a CIGAR score of 7M1I1M2I4M1D5M which I believe should have been 'filtered out' given the bwa flags shown above. While I still don't see the alignment, the CIGAR score likely should have been filtered despite it being a suboptimal hit (unless I'm not understanding the suboptimal hit thresholds correctly).
Any help will be greatly appreciated. I'm mainly trying to understand if we should continue to use BWA for reading 'non-uniq' reads, use some other method, or apply some work-arounds. Could the problems be only with sub-optimal hits? Or could it be more wide-spread and effect all repeat-hits or even uniq hits?
Thank you.
Comment