Hello everybody,
Could anyone tell me how to get SNPs using SAMtools???. I am learning how to do it but can´t get the SNPs, I only get the Indels. I am working with a mutant I created from an E coli genome (just for practice). The mutant was created with 5 mismatches, 5 insertions and 5 deletions (using Biopieces).
The lenght and coverage is here:
read_fasta -i <fasta ref genome>.fna | shred_seq -s 100 -c 50 | add_ident -k SEQ_NAME | write_fasta -x -o <rrreads>.fna
The alignment was made with BWA . Then I used SAMtools protocol (as described in some of the posts) like this :
samtools mpileup -uf <refgenome>.fasta result.sort.bam > results.bcf
and then ,
bcftools view -vcg results.bcf - > test.vcf
If anyone can help me .. Thanks!!! :)
•
written
by
Pilar E. •
1 answer
First simple thing to try: samtools mpileup -Buf. Sometimes, the BAQ calculations will wrongly believe real SNPs to be undiagnosed indels, and will drop their quality in the pileup step such that they can't be called as SNPs
Since you know where your missing SNPs are, look at them in the pileup text file. If the SNP is detected, but the quality score is terrible, despite being fine in the .bams, that's what happened.
Log in to answer this question.
-
Mpileup Command Line
written by Giffredo •
-
Consensus Fasta Sequence Of A Given Genomic Region For Individuals From 1000G
written by Pierre
-
Mpileup Problem
written by stolarek.ir
-
Extract Fasta Sequences From File Which Are Not Present In Other Fasta File
written by namratapatel183 •
-
What Better Way To Get Consensus Sequence Of Mapping Against Reference Genome Using Samtools Mpileu…
written by thiagomafra
-
Error Using Samtools To View Content In .Sam File
written by soosus •
-
Mpileup Output And Quality Scores
written by NextGenSeek •
-
Emule "Samtools View" Using Pysam
written by Geparada
-
Samtools To Analyze Bwa Output
written by Dawnoflife •
-
Error Running Blast2Sam.Pl
written by Rm
might be the alignment step, too. did you enable mismatches for the alingments? Because, if just by chance only indels were allowed but 0 mismatches, then you cannot call SNPs. Can you provide the bwa parameters?
Thanks for answering !Here is the BWA :
read_fasta -i reference.fasta | mutate_seq -n 5 | indel_seq -i 5 -d 5 | write_fasta -xo mutant.fna read_fasta -i mutant.fna | shred_seq -s 100 -c 50 | add_ident -k SEQ_NAME | write_fasta -xo reads.fna mkdir BWA (for alignment files) bwa index -p BWA/Ecoli reference.fasta bwa aln BWA/Ecoli reads.fna > reads.sai bwa samse -f result.sam BWA/Ecoli reads.sai reads.fna
Do you see anthing wrong??? Thanks for helping Jeremy! :)
Pilar,
I am having trouble with this same situation.
I have a VCF file with only indels. Did you figure out why this occured?
Jim