---

title: “Variant calling tutorial” output: html_notebook

#script for Variant calling tutorial
#Calling variants in reads mapped by bowtie
#Load SAMtools
#Prepare your directories
#Index the FASTA reference file
#Convert mapped reads from SAM to BAM, sort, and index
#Call genome variants
#Calling variants in reads mapped by BWA or Bowtie2
#Filtering VCF files with grep
#Comparing the results of different mappers using bedtools
module load samtools
samtools
which samtools
mkdir samtools_bowtie
cp bowtie/SRR030257.sam samtools_bowtie
cp bowtie/NC_012967.1.fasta samtools_bowtie
#Index the FASTA reference file
#First, you need to index the reference file. (This isn't indexing it for read mapping. It's indexing it so that SAMtools can quickly jump to a certain base in the reference.)Then run this command to index the reference file.

samtools faidx samtools_bowtie/NC_012967.1.fasta
#Take a look at the new *.fai file that was created by this command.
less samtools_bowtie/NC_012967.1.fasta.fai

#Convert mapped reads from SAM to BAM, sort, and index
#SAM is a text file, so it is slow to access information about how any given read was mapped. SAMtools and many of the commands that we will run later work on BAM files (essentially GZIP compressed binary forms of the text SAM files).
samtools view -b -S -o samtools_bowtie/SRR030257.bam samtools_bowtie/SRR030257.sam
#Sort and index the BAM file.
samtools sort samtools_bowtie/SRR030257.bam samtools_bowtie/SRR030257.sorted
samtools index samtools_bowtie/SRR030257.sorted.bam
#Call genome variants
#we use the mpileup command from samtools to compile information about the bases mapped to each reference position.Output BCF file. This is a binary form of the text Variant Call Format (VCF).

samtools mpileup -u -f samtools_bowtie/NC_012967.1.fasta samtools_bowtie/SRR030257.sorted.bam > samtools_bowtie/SRR030257.bcf
#Convert BCF to human-readable VCF:
  bcftools view -v -c -g samtools_bowtie/SRR030257.bcf > samtools_bowtie/SRR030257.vcf
  
  #Comparing the results of different mappers using bedtools
  mkdir comparison
cp samtools_bowtie/SRR030257.vcf comparison/bowtie.vcf
cp samtools_bwa/SRR030257.vcf comparison/bwa.vcf
cp samtools_bowtie2/SRR030257.vcf comparison/bowtie2.vcf
cd comparison
 #load module
 module load bedtools
#Finding common mutations.
bedtools intersect -a bowtie.vcf -b bwa.vcf > common_bowtie_bwa.vcf
#Finding mutations that are unique for each mapper.
bedtools subtract -a bowtie.vcf -b common_bowtie_bwa.vcf > unique_bowtie.vcf
bedtools subtract -a bwa.vcf -b common_bowtie_bwa.vcf > unique_bwa.vcf