Table des matières
Here we will use the BWA aligner to map short reads to a reference genome, and then call variants (differences between the reads and the reference).
Download data
Goal: get the sequence data!
Run:
mkdir ~/data cd ~/data curl -O http://dib-training.ucdavis.edu.s3.amazonaws.com/2017-ucdavis-igg201b/SRR2584857.fq.gz
Map data
Goal: execute a basic mapping
Run the following commands to install bwa:
cd curl -L https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2/download > bwa-0.7.15.tar.bz2
tar xjvf bwa-0.7.15.tar.bz2 cd bwa-0.7.15 make
sudo cp bwa /usr/local/bin echo 'export PATH=$PATH:/usr/local/bin' >> ~/.bashrc source ~/.bashrc
2 Make & change into a working directory:
3 Copy and gunzip the reference:
wget https://github.com/ctb/2017-ucdavis-igg201b/raw/master/lab1/ecoli-rel606.fa.gz gunzip ecoli-rel606.fa.gz
4 Prepare it for mapping:
bwa index ecoli-rel606.fa
5 Map! (This will take about 2 minutes.)
bwa mem -t 6 ecoli-rel606.fa ~/data/SRR2584857.fq.gz > SRR2584857.sam
6 Observe!
Visualize mapping
Goal: make it possible to go look at a specific bit of the genome.
Install samtools:
sudo apt-get -y install samtools
Convert the SAM file into a BAM file that can be sorted and indexed:
samtools view -hSbo SRR2584857.bam SRR2584857.sam
Sort the BAM file by position in genome:
samtools sort SRR2584857.bam SRR2584857.sorted
Index the BAM file so that we can randomly access it quickly:
samtools index SRR2584857.sorted.bam
Visualize with tview:
samtools tview SRR2584857.sorted.bam ecoli-rel606.fa
tview commands of relevance:
left and right arrows scroll
q to quit
CTRL-h and CTRL-l do “big” scrolls
g ecoli:3931002 will take you to a specific location.
Call variants!
Goal: find places where the reads are systematically different from the genome.
Now we can call variants using samtools mpileup:
samtools mpileup -uD -f ecoli-rel606.fa SRR2584857.sorted.bam | \
bcftools view -bvcg - > variants.raw.bcf
This will take a few minutes… and output a file that is not human readable! But we can quickly convert it into the ‘variant call format’ that is human readable:
bcftools view variants.raw.bcf > variants.vcf
Look at the VCF file
The output VCF file contains a list of all the variants that samtools thinks are there. What’s in it?
The official VCF specification is a great read…if you’re suffering from insomnia. Let’s skip this and just take a quick look at the file.
Look at the non-commented lines along with the header:
The first five columns: CHROM POS ID REF ALT. It’s a little easier to see if you run
grep -v ^## variants.vcf | less -S
Use your left and right arrows to scroll, and ‘q’ to quit.
Examine one of the variants with tview:
samtools tview SRR2584857.sorted.bam ecoli-rel606.fa -p ecoli:920514
‘q’ to quit, left arrow to scroll a bit left.
Well, at least that variant looks real… Look at the VCF file with bedtools.¶
bedtools docs
Download and build bedtools:
cd ~/ curl -O -L https://github.com/arq5x/bedtools2/releases/download/v2.26.0/bedtools-2.26.0.tar.gz tar -xzf bedtools-2.26.0.tar.gz
cd bedtools2 make sudo make install
Go back to work:
Download a GFF3 file with annotations for E. coli:
wget https://github.com/ctb/2017-ucdavis-igg201b/raw/master/lab3/ecoli-rel606.gff.gz .
Run bedtools intersect:
bedtools intersect -a ecoli-rel606.gff.gz -b variants.vcf -wa -u
Documentation for bedtools intersect
Extract reads with samtools.
Execute:
samtools view SRR2584857.sorted.bam 'ecoli:920514-920514' > out.sam wc -l out.sam
and this will give you the coverage of the relevant position. Discussion points / extra things to cover¶
What are the drawbacks to mapping-based variant calling? What are the positives? Where do reference genomes come from?
Source : https://angus.readthedocs.io/en/2017/variant-calling.html