I aligned my fasta files to the reference genome using bwa sampe.
I have paired end reads so first used:
bwa aln -t 2 -l 20 lyrata_genome.fa Trm_MA605_R1.fastq > Trm_MA605_R1.sai
bwa aln -t 2 -l 20 lyrata_genome.fa Trm_MA605_R2.fastq > Trm_MA605_R2.sai
Then I did the sampe alignment using:
bwa sampe \
-n 1 \
-N 2 \
lyrata_genome.fa \
Trm_MA605_R1.sai \
Trm_MA605_R2.sai \
Trm_MA605_R1.fastq \
Trm_MA605_R2.fastq > SAMPEaligned_MA605readsBWA.sam
I tried to sort the alignment and convert it to bam using picard. But, getting the error message (part of the total error message):
Exception in thread "main" htsjdk.samtools.SAMFormatException: Error parsing text SAM file. MAPQ should be 0 for unmapped read.; File SAMPEaligned_MA605readsBWA.sam; Line 673702
Line: HWI-ST588:81:C080WACXX:7:1102:6980:187157 69 4 23328260 6 84M 6 963719 0 GGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGGTTTAGGCATAAA +2AAEEEIIFFIFIEI<CEDCEDDDDDIIIBDDEIDDBDDEIIDCCDIIII=ADDDDD?@DAAAAA>?AAAA>??AAA35(54> RG:Z:MA605_C080WACXX_7 XT:A:U NM:i:3 SM:i:6 AM:i:6 X0:i:1 X1:i:55 XM:i:3 XO:i:0 XG:i:0 MD:Z:28C49A1A3
at htsjdk.samtools.SAMLineParser.reportErrorParsingLine(SAMLineParser.java:439)
at htsjdk.samtools.SAMLineParser.parseLine(SAMLineParser.java:341)
at htsjdk.samtools.SAMTextReader$RecordIterator.parseLine(SAMTextReader.java:248)
I tried to validate the sam file using picard again:
java -jar /home/everestial007/picard-tools-1.139/picard.jar ValidateSamFile \
I=SAMPEaligned_MA605readsBWA.sam \
Mode=VERBOSE \
MAX_OUTPUT=15000 \
O=validatebamMA605forGATK
But I'm getting the following error message in the output file.
ERROR: Record 46, Read name HWI-ST588:81:C080WACXX:7:1101:3935:2321, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 155, Read name HWI-ST588:81:C080WACXX:7:1101:9306:2386, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 352, Read name HWI-ST588:81:C080WACXX:7:1101:20820:2318, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 382, Read name HWI-ST588:81:C080WACXX:7:1101:2880:2595, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 562, Read name HWI-ST588:81:C080WACXX:7:1101:10717:2694, Mate negative strand flag does not match read negative strand flag of mate
ERROR: Record 603, Read name HWI-ST588:81:C080WACXX:7:1101:12062:2750, Mate negative strand flag does not match read negative strand flag of mate
I followed the bwa manual to letter but getting this kind of errors. I can not figure why is my alignment not working.
Thanks in advance.
Thanks for the info. But, should not we be considering the error message. I don't understand what "Mate negative strand flag does not match read negative strand flag of mate" means, and if I should be concerned that my alignments are working properly or not. Btw, sam flagstat reports that 86 % of the reads are aligned - so it is something cocerning.
Thanks,
It's possible that your mates are out of sync (i.e., the nth read in Trm_MA605_R1.fastq doesn't match the nth read in Trm_MA605_R2.fastq), so you might check that (you can use reformat from BBMap to do this).
Thanks for the info. I a little hesistant in trying to learn a new tool again (just takes too much of my time to set it up to getting it right), but will do it as a last resort. Do you think picard can do it? Also, if if have 86% alignment, should I be worried about it? - just wanting to weigh the pros and cons of the extra step, but I will appreciate any suggestions.
Thanks,
If the overwhelming majority are marked as "proper pairs" then I wouldn't bother with further processing, the results are likely fine.
Thanks for enlighting. Here is my samtools flagstat report,
22193026 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
19848288 + 0 mapped (89.43% : N/A)
22193026 + 0 paired in sequencing
11096513 + 0 read1
11096513 + 0 read2
18519235 + 0 properly paired (83.45% : N/A)
19450790 + 0 with itself and mate mapped
397498 + 0 singletons (1.79% : N/A)
366420 + 0 with mate mapped to a different chr
122636 + 0 with mate mapped to a different chr (mapQ>=5)
Since, 83.45% are properly paired so I think I am good to go, but just want to make sure. Also, I used picard's "CleanSam" to remove reads that mapped out of reference sequence.
That seems reasonable, you should be good to go.