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:
#Theta calculation
Running estimation of theta is very simple now that we've gotten everything set up. We need to modify one line in the Thetas_Calculator.sh script:
TAXON=test
If you would like to change the window size and step amount, you can modify the following variables:
WIN=100
STEP=50
You can also turn off sliding windows with:
SLIDING_WINDOW=false
Great! Now we can run the script:
sbatch -p bigmemm scripts/Thetas_Calculator.sh scripts/Thetas_Calculator.conf
We can visualize this the same way as before with the site frequency spectrum. The file we are interested in is called test_Diversity.thetas.gz.pestPG:
$scp name@server.address.edu:~/rilab/aw-tutorial/results/test_Diversity.thetas.gz.pestPG ./
Once you upload the file into the web app, you should get two graphs that look like this:
You can view these graphs interactively by using the controls on the side. You can switch between estimators of theta and neutrality statistics with the drop down menus and zoom in on regions of interest with the Base start and end position number fields. You can also add lowess curves to the graphs and include gene annotations by uploading a GFF3 file that corresponds to your species.
In order to carry out admixture analysis, we have to compile the ngsAdmix program.
cd ngsPopGen
g++ ngsadmix32.cpp -O3 -lpthread -lz -o NGSadmix
Then we just have to change one variable in the Admixutre.conf file.
TAXON=test
The variable K specifies the ending point for the number of ancestral populations to run admixture with (starting ith 2). For example, if K=5, admixture will be run from K={2,3,4,5}. Now we can submit it like this:
sbatch -p bigmemm scripts/Admixture.sh scripts/Admixture.conf
Once this has finished running, you can copy over the results (*.qopt) and view them in the graphing application.
$scp name@server.edu:~/rilab/aw-tutorial/test.3.qopt ./
The results should look something like this:
