call indel and SNP by varscan

[doublealice]

Hi all,

I met problem when I work on a 7M bacterium genome. I mapped 140x of paired-end reads (151x2) from illumina to a Pacbio reference genome by BWA mem. Then used samtool to get a mpileup file. Finally using varscan to call indel and SNP. Below is commands I used:

map:
$bwa mem -M -t 6 Pacbio_ref reads_1.fq reads_2.fq | samtools view -buS - | samtools sort - illumina_paired_map_pacbio.sorted
$samtools index illumina_paired_map_pacbio.sorted.bam
[/CODE]

get mpileup

$samtools mpileup -f Pacbio_ref.fastq illumina_paired_map_pacbio.sorted.bam > mpileup_output

call SNP

java -Xmx6g -Djava.io.tmpdir=temp -jar VarScan.v2.3.9.jar  \
mpileup2snp mpileup_output \
--min-coverage 20  \
--min-reads2 2 \
--min-avg-.01 \
--p-value 0.05  \
--output-vcf 1 > 05snp.vcf

After that, I got only 268 SNP, similar work for indel but only 78 were found.

I don't think those numbers should be so low. Would you please help to check if there is anything wrong on any step? I think maybe the parameters in BWA were not set correctly. Please help!

Thanks!

[Jane M]

Don't know if it can solve your problem but you can try to add -B option to samtools mpileup, as discussed here for example:

http://massgenomics.org/2012/03/5-things-to-know-about-samtools-mpileup.html

http://seqanswers.com/forums/showthread.php?t=18462&highlight=varscan2