Introductory RAD-seq Activity


Today, let’s go through a brief experimental design activity to help demonstrate some of the features of RAD-seq and how they can be harnessed to explore numerous questions in ecology, evolutionary biology, and genetics. To perform this activity, we will use Harvard’s Odyssey compute cluster.

Background & Resources

To introduce the topic of RAD-seq, we are relying on the excellent review article by Andrews et al. 2016. Nature Reviews Genetics. We are also using Catchen et al. 2013. Molecular Ecology as an empirical example of RAD-seq for discussion. Slides of key figures from these two resources are available here.

Setting up Our Analysis Environment

Given we cannot run jobs on the login nodes, we are going to start an interactive job that will allow us to go through this lesson. You can execute the following command to do this. It may take a few minutes to clear the queue, but eventually your prompt should return.

srun -p test --pty -n 2 --mem 4000 -t 0-06:00 --account=oeb275r /bin/bash

We are going to use the very useful environments feature of Anaconda Python to manage the software dependencies we have for this activity. Anaconda is already installed on Odyssey and can be loaded as follows:

module load Anaconda3/5.0.1-fasrc02

Now we can create an analysis environment that will always be available for us to use for this activity, which will run using Python version 3.7.

conda create -n radseq python=3.7

This will take a few minutes to run and along the way, we will answer Yes (y) to allow installation. When it is done, we can activate our newly created environment.

source activate radseq
# on some systems:
# conda activate radseq

We must also install some other important software in our environment, namely Biopython. For whatever reason, installing biopython through conda is not allowing later software to run properly for me, so I have commands for installing via conda (preferred) and using pip instead, which is Python’s default package manager.

# preferred
 conda install -c conda-forge biopython 
# only if conda does not work:
# pip install --user biopython

We must also install bedtools, which is required for some of our software. This is available on Odyssey already, but let’s install it into our radseq environment to demonstrate how to install using conda.

conda install -c bioconda bedtools

Now let’s create a directory for our activity called radseq_activity and move into it so we can do our work there.

mkdir radseq_activity
cd radseq_activity

We are also going to be using some software that I wrote for performing in silico RAD-seq called RADis. I have, unfortunately, not packaged this into easily-installable software, but we can still make use of it fairly easily. Let’s download the GitHub repository for this purpose.

git clone

One more thing we will need to do is add the software in this repository to our software $PATH, which is where the computer looks for software to run by default. Note that this is a temporary addition that will not persist beyond this session or if this directory is moved. There are ways to permanently add these software to your $PATH that I can show if anyone is interested.

cd RADis
export PATH=$PATH:`pwd`
cd ..

In Silico RAD-seq Activity

In order to perform our in silico RAD-seq, we need a sequence for our digital restriction enzyme digestion. For time and ease, we will use chromosome 22 from the human genome, which is about 50 Mb in length. We can download this from the UCSC Genome Browser database.

gunzip chr22.fa.gz

We need to create a restriction enzymes file to run the in silico restriction digest. This software performs a double digest, so two enzymes need to be used. We are going to work in 4 groups, each with a different combination of enzymes as follows.

# group 1
sau3AI	GATC	0

# group 2
sau3AI	GATC	0

# group 3
mspI	CCGG	1

# group 4
mspI	CCGG	1

Each group should take their restriction enzymes and paste them into a new text file called enzymes.txt. This can be done using the terminal-based text editor nano. Note that there are TAB characters between each of the columns.

Once we have our reference genome (human chromosome 22) and our enzymes, we can run an in silico digestion with each enzyme using the following command. Sit back and be patient, as I’m not a professional software engineer and my program is slower than I would like (should take about 10 minutes at most). --input chr22.fa --output chr22 --enzymes enzymes.txt

Once the command ends, you will see a couple new files in the working directory when we run ls.

# group 1 output
# chr22_sbfI_CCTGCAGG_8.bed
# chr22_sau3AI_GATC_4.bed

Now that we have the individual restriction digest files for single enzymes, we can run the double digest. As part of this command, we need to also set the lower and upper limit on the fragment lengths we will select for. Given there are 4-5 people in a group, let’s have each person select a different range of fragment sizes.

# person 1
200 - 400
# person 2
300 - 500
# person 3
400 - 600
# person 4
500 - 700
# person 5
600 - 800

Now we can run our double digest as follows. Be sure to use the output from the longer enzyme as the -r input. This will run much faster than the other command.

# example for group 1, person 1 -g chr22.fa -r chr22_sbfI_CCTGCAGG_8.bed -c chr22_sau3AI_GATC_4.bed -l 200 -u 400 -o chr22_sbfI_sau3AI_200_400_output

You will see two outputs of this command.

# example for group 1, person 1
# chr22_sbfI_sau3AI_200-400_output.bed.gz
# chr22_sbfI_sau3AI_200-400_output.fasta.gz

Both are gzipped, but we can easily unzip them.

gunzip *.gz

And look at their contents. One is a BED file.

head chr22_sbfI_sau3AI_200_400_output.bed
# chr22	10520281	10520562	sbfI-sau3AI	281	+
# chr22	10527904	10528214	sbfI-sau3AI	310	-
# chr22	10528218	10528432	sbfI-sau3AI	214	+
# chr22	11054755	11054977	sbfI-sau3AI	222	-
# chr22	11265186	11265460	sbfI-sau3AI	274	+
# chr22	11438712	11439070	sbfI-sau3AI	358	+
# chr22	11553161	11553482	sbfI-sau3AI	321	+
# chr22	11716993	11717223	sbfI-sau3AI	230	+
# chr22	11953435	11953797	sbfI-sau3AI	362	-
# chr22	11956673	11957032	sbfI-sau3AI	359	-

The other is a FASTA file.

head chr22_sbfI_sau3AI_200_400_output.fasta
# >chr22:10520281-10520562(+)
# ggtggggcttcactggggacccatcaccttccacccaggagcctgtctgcctcccacaaccatccatggcacccaggctgctggcatcaaggggcacttgcaggccagtgccaagccaccctcgtaccccctcatcttcccctcccatgctcctgctcctcagtgtccaaagtccagaaggggctgaggtggcaggggactgacatgtcagcactgcttccaatgtgtgcacacctggctgggcagtgacagcaccctgctgggtcccaaccccactctga
# >chr22:10527904-10528214(-)
# >chr22:10528218-10528432(+)
# >chr22:11054755-11054977(-)
# >chr22:11265186-11265460(+)

Now the critical piece of information we need for experimental design considerations is the number of fragments, as marker number/density enables different types of questions to be addressed. We can easily determine this by counting the number of lines in our BED file.

wc -l chr22_sbfI_sau3AI_200_400_output.bed
# 952

Or the number of sequences in our FASTA file.

grep -c ">" chr22_sbfI_sau3AI_200_400_output.fasta
# 952

Activity Questions

Within your group, work together to determine how fragment number varies across different fragment size selections. What pattern do you see and why might this be important?

As a class, what do you notice about the number of fragments sampled based on the restriction enzymes and the fragment size windows used? What aspect of the restriction enzyme appears to impact fragment number?

Based on your research interests, think of an interesting question (or questions) you might want to ask and how genome-wide data could be useful. Assuming you are restricted to using RAD-seq only, what experimental design considerations must you consider to answer your research question?