Identifying SNPs using mpileup
A quick way to call SNPs is to use mpileup which is one of the tools provided by samtools. However if you want control over your results, consider the GATK pipeline found Processing files and Running GATK|here. The steps to process the file are almost the same compared to the GATK pipeline, however mpileup allows you to provide the different groups (Read Groups) as separate bam files as well.Below is a recommended pipeline when working with mpileup.
First create the bowtie index files for alignment.
module load bowtie2 module load samtools bowtie2-build PTC_Human.fasta PTC_Human samtools faidx PTC_Human.fasta
The following steps have to executes for all samples
- align using Bowtie2
- convert to bam
- sort the bam
- remove duplicates
- create a bam index for dedup bam files
bowtie2 -x PTC_Human -U Cohen.fastq -S Cohen.sam samtools view -bS Cohen.sam > Cohen.bam samtools sort Cohen.bam Cohen_sorted samtools index Cohen_sorted.bam
For further curation
- realign indels within each samples
Once they are all processed you can create read groups for all the samples and then merge then into one file. This is not necessary for mpileup but it is for GATK. So if you are planning to use both methods it is not a bad idea.
- create read groups
- merge all samples to one bam file
- sort and index this merged bam file
Running mpileup assuming no read groups and not merged
- run mpileup
- convert bcf to vcf
- filter using vcfutils
- index vcf file to visualize on IGV
samtools mpileup -uf PTC_Human.fasta \ Cohen.bam \ Linder.bam \ Rikhi.bam \ Sherman.bam > PTC_human.bcf
Run bcf tools to call the snps:
- b means save output as bcf
- v means output potential variant sites only
- c means call snps
- g call genotypes at variant sites.
bcftools view -bvcg PTC_human.bcf > PTC_human.raw.bcf &