-
Notifications
You must be signed in to change notification settings - Fork 12
Tutorial
Welcome! This is a short guide to population genetics analysis using angsd-wrapper. We will be using a test data set containing sequence from Oryza sativa and Oryza glumaepetula. First, we need to clone the angsd-wrapper repository. You need to have git installed to do this. Alternatively, you can download a zip file from the releases page or use the download zip button on the home page of the repository.
$ git clone https://github.com/arundurvasula/angsd-wrapper.git aw-tutorial
$ cd aw-tutorial/
angsd-wrapper comes with its own version of ANGSD to prevent compatibility breaking changes in ANGSD from affecting angsd-wrapper and comes with a few other programs. In order to compile these programs, you must run make:
$ cd angsd; make
$ cd ../ngsPopGen/; make
$ cd ../ngsF/; make
Great! Now we're ready to get the example data.
$ cd ../data
$ wget http://de.iplantcollaborative.org/dl/d/3A541C91-A66A-4651-949D-4E65028C4A2F/iplant.zip
$ unzip -j iplant.zip
The example data should now be in the data folder. Now we need to change the Common_Variables.conf file to point to the correct paths.
$ nano scripts/Common_Variables.conf
Modify this line to point to your cloned directory:
PROJECT_DIR=${HOME}/rilab/aw-tutorial
And these two lines to point to the reference and ancestral sequences.
# Reference genome sequence
REF_SEQ=${DATA_DIR}/reference.Oryza_sativa.IRGSP-1.0.23.dna.genome_chr.fa.gz
# Ancestral sequence
ANC_SEQ=${DATA_DIR}/ancestral.merid_japonica_chr.fa.gz
Great! Now we're ready to calculate a site frequency spectrum on the data. The script to do this is called scripts/Site_Frequency_Spectrum.sh and the variables are stored in scripts/Site_Frequency_Spectrum.conf
In the configuration file, we need to point the TAXON variable to our test data:
TAXON=test
TAXON_LIST=${DATA_DIR}/${TAXON}_samples.txt
TAXON_INBREEDING=${DATA_DIR}/${TAXON}_F.txt
Now we can run the site frequency spectrum calculation by running the script:
$ sbatch -p bigmemm scripts/Site_Frequency_Spectrum.sh scripts/Site_Frequency_Spectrum.conf
You should get a few files in your results directory. The one that contains the site frequency spectrum to graph is called test_DerivedSFS. The other files can be used as inputs to other commands. In order to graph the site frequency spectrum, we need to open up the shiny web app. If the server you are connected to has a web browser, you can do it directly on there. If it doesn't, you need to copy the results onto your local computer and run the web app there.
If this is the case, you need to clone angsd-wrapper on your local computer:
$ git clone https://github.com/arundurvasula/angsd-wrapper.git aw-tutorial
$ cd aw-tutorial
And copy your results file (I use scp):
$ scp name@server.address.edu:~/rilab/aw-tutorial/results/test_DerivedSFS ./
You need to make sure you have the dependencies installed:
> install.packages("shiny")
> install.packages("ape")
> install.packages("Hmisc")
> install.packages("lattice")
> source("http://bioconductor.org/biocLite.R")
> biocLite("genomeIntervals")
> install.packages("data.table")
Once they are installed, you can run the shiny app:
$ R
> library(shiny)
> runApp("scripts/shiny")
A webpage should open up in your browser. To graph the site frequency spectrum, click the SFS tab and choose your file. Once you upload the file you should get a graph of your site frequency spectrum. You can save this graph by right clicking and saving the image. If you are following the tutorial, you should get a funky looking graph like this:
New tutorial outline (comments welcome)
- clone angsd-wrapper
- wget data (http://de.iplantcollaborative.org/dl/d/3A541C91-A66A-4651-949D-4E65028C4A2F/iplant.zip)
- set up common.conf
- run SFS
- run theta
- visualize both with shiny
old:
Welcome! This page will walk you through basic population genetic analysis using ANGSD and angsd-wrapper.
There are 3 types of data needed for analysis in ANGSD:
- A reference sequence (
.fasta) - An ancestral reference sequence (
.fasta) - Individual sequences mapped to the reference sequence (
.bam)
In order to simplify the input, ANGSD requires that the file paths to the bam files be placed in one file. For example, the data/BKN_samples.txt included in the repository has a list of paths to bam files. Thus, to get started with ANGSD analysis, our data folder will contain the following files:
data/{TAXON}_samples.txt
data/refseq.fasta
data/ancseq.fasta
NOTE: The {TAXON} part of the samples file is important. angsd-wrapper will look for files named exactly like that. See the wiki page for more information.
Now we are ready to do our first analysis, a site frequency spectrum. In ANGSD, this takes 2 steps that have been simplified in angsd-wrapper:
- create a binary file containing a prior for the SFS
- obtain a maximum likelihood estimate of the SFS
Using angsd-wrapper, we can accomplish this by modifying the included configuration file (or creating a new one) to include the path to the ancestral sequence, the reference sequence, and the TAXON variable (which will point to the list of bams).
Once you have modified the configuration file to include this information, you can begin the analysis. If you are just running the analysis on your local computer, you can begin it by running (from the angsd directory):
$ bash scripts/ANGSD_SFS.sh scripts/sfs_example.conf
If you are running it on a cluster, you will likely need to create a submission script. In slurm, you can do something like the following:
sfs_submit.sh:
#!/bin/bash
#SBATCH -D /home/adurvasu/angsd-wrapper
#SBATCH -J Slurm-SFS
#SBATCH -o /home/adurvasu/angsd-wrapper/results/out-%j.txt
#SBATCH -e /home/adurvasu/angsd-wrapper/results/error-%j.txt
bash scripts/ANGSD_SFS.sh scripts/sfs_example.conf
Then submit the script to your cluster's queue:
$ sbatch -p serial scripts/sfs_submit.sh