Kalain Hündin

Cette page est en lecture seule. Vous pouvez afficher le texte source, mais ne pourrez pas le modifier. Contactez votre administrateur si vous pensez qu'il s'agit d'une erreur.

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