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
call SNP
After that, I got only 268 SNP, similar work for indel but only 78 were found. ![Frown](https://www.seqanswers.com/core/images/smilies/frown.png)
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!
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
Code:
$samtools mpileup -f Pacbio_ref.fastq illumina_paired_map_pacbio.sorted.bam > mpileup_output
Code:
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
![Frown](https://www.seqanswers.com/core/images/smilies/frown.png)
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!
Comment