parameter set of bwa-samtools for snp calling

[heiya]

hi, everyone! I use bwa-samtools for my rice genome snp calling recently, However I found in the "mpileup" step, the snp calling result is significantly different according to the related parameter set.eg. no -6A, the variation is 66041, however, used -6A ,the variation is 181319(triple than before), is that due to the anomalous reads?
The result let me confused on how to set the parameter to get the accurate result?

[heiya]

I have tested, that's due to the -6, it means Illumina 1.3+ format. why?
I have check my fastq file, it is Illumina 1.5+ format, should I add -6?
All the reply is appreciated!

[vv85]

The answer is yes, you should use the -6 parameter if your reads are in Illumina 1.3+ format. Also when aligning with bwa did you use the -I option? If not your mapping qualities may not be accurate and that may affect the SNP calling

[heiya]

The answer is yes, you should use the -6 parameter if your reads are in Illumina 1.3+ format. Also when aligning with bwa did you use the -I option? If not your mapping qualities may not be accurate and that may affect the SNP calling

Could you say more in detail？

[vv85]

For more info on the two different fastq formats you can use the wiki page:

http://en.wikipedia.org/wiki/FASTQ_format

So for example, let's say you have this:
@HWUSI-232323
ATTTTA
+
BBBBBB

If the format is Phred+33 then you have a good read. If the format is Phred+64 then you have a bad read. Mapping qualities depend on base qualities so by not specifying the actual format mapping qualities will not reflect correctly what's in the data. Since SNP callers use this MAPQ for variant calls, this leads to differences.

[haojam]

parameter set of bwa-samtools for snp calling

Hi,
Everyone, I use bwa-samtools for human genome SNPs calling, but the SRA READS are performed using the Illumina Genome Analyzer Pipeline v.1.4.0 software suite. Should I need to run BWA alignment without -I or -L or -J options for diferent illumina versions. When I ran with -I option there was error saying sequence and quality are inconsistent. I would be glad and appreciate for your consideration.

Rocky

[vv85]

I wouldn't think that specific problem is related to the parameter. Do you get a different result without using the -I option? At which step does the error occur?

[heiya]

I wouldn't think that specific problem is related to the parameter. Do you get a different result without using the -I option? At which step does the error occur?

hi,vv85
I had checked the result of using -I and without using -I when run bwa, there is no difference in sam format report, except that when using -I , it converted the quality score of Illumina format into sanger format in the reads quality column. However, without using -I, it didn't do this job.

[vv85]

That's more or less what's expected, depending on the overall quality of your reads some mapping qualities could be changed too. With these reads now should be able to use samtools without the -6 and get (somewhat) confident SNPs

[haojam]

When I try running with -I option at BWA alignment and converted to sam file then to bam file error occurs.

These are the comment
SEQUENCE ALIGNMENT
$ ./bwa aln -I -t 4 /reference_genome/hg18bwaidx ~/SRR000001_1.fastq > ~/SRR000001_1.sai
$ ./bwa aln -I -t 4 /reference_genome/hg18bwaidx ~/SRR000001_2.fastq > ~/SRR000001_2.sai


PARIED-END MAPPING
$ ./bwa sampe /reference_genome/hg18bwaidx ~/bwa-0.5.9/SRR000001_1.sai ~/bwa-0.5.9/SRR000001_2.sai ~/SRR000001_1.fastq ~/SRR000001_2.fastq > ~/SRR000001.sam


CONVERT TO BAM FILES
$ ./samtools view -b -h -S /SRR000001.sam > ~/SRR000001.bam

ERROR message - sequence and quality are inconsistent

[haojam]

paired-end reads were obtained from GAII (Illumina Genome Analyzer Pipeline v.1.4.0)