Table of Contents¶

Monday, Day 1:¶

11am: Meet in Valley Hall (room TBA), setup, and introductions.

• Logging into the cloud (XSEDE Jetstream)

Noon - 1:15: Lunch

1:15 - 2pm: Group meeting with all of DIBSI

2pm - 4pm:

• Working on Command Line
• Command Line Blast

Homework: Skim the Critical Assessment of Metagenome Interpretation (CAMI) Paper

Tuesday, Day 2:¶

9am - Noon:

• Introduction to dataset we will be using
• Evaluating your short-read data set quality Noon - 1:15: Lunch

1:15 - 4pm:

• Lecture about assembly from Titus Brown
• Trying quality trimming and assembly with your own data!

Wednesday, Day 3:¶

9am - Noon:

• Assembling your short read data set with MEGAHIT
• Quickly searching and comparing your samples with sourmash

Noon - 1:15: Lunch

1:15 - 4pm:

• Mapping short reads to the assembly time permitting
• Binning genomes out of your metagenome

Thursday, Day 4:¶

9am - Noon:

• Annotating your dataset with Prokka
• Quantifying abundance across samples with Salmon

Noon - 1:15: Lunch

1:15 - 4pm:

• Anvi’o
• A brief discussion of workflows & repeatability

Friday, Optional Day 5: TBD (plotting, visualization, jupyter notebooks...)¶

What next?

Resources:

SEQ Answers

Biostars

Data Carpentry

DIB Summer Institute

Undone¶

See the complete table of contents

Technical information¶

The github repository for this workshop is public at https://github.com/ngs-docs/2017-ucsc-metagenomics

1. Learning goals¶

For you:

• get a first (or second) look at tools;
• gain some experience in the basic command line;
• get 80% of way to a complete analysis of some data;
• introduction to philosophy and perspective of data analysis in science;

2. Safe space and code of conduct¶

This is intended to be a safe and friendly place for learning!

Please see the Software Carpentry workshop Code of Conduct: http://software-carpentry.org/conduct.html

In particular, please ask questions, because I guarantee you that your question will help others!

3. Instructor introductions¶

Harriet Alexander - postdoc at UC Davis.

Phil Brooks - postdoc at UC Davis.

Titus Brown - prof at UC Davis in the School of Vet Med.

Shannon Joslin - grad student at UC Davis.

Taylor Reiter - grad student at UC Davis.

4. Jetstream and cloud computing - why?!¶

• simplifies software installation;
• can be used for bigger analyses quite easily;
• good for “burst” capacity (just got a data set!)
• accessible everywhere;

5. Sticky notes and how they work... + Minute Cards¶

Basic rules:

• no sticky note - “working on it”
• green sticky note - “all is well”
• red sticky note - “need help!”

Place the sticky notes where we can see them from the back of the room – e.g. on the back of your laptop.

Next: Booting a Jetstream Computer Instance for your use!

Request to log in to the Jetstream Portal¶

Click the login link in the upper right.

Use “XSEDE”¶

Choose “XSEDE” as your account provider (it should be the default) and click on “Continue”.

Fill in the username and password and click “Sign in”¶

Fill in the username, which is ‘tx160085’ for the ANGUS workshop, and then enter the password (which we will tell you in class).

Select Projects and “Create New Project”¶

Now, this is something you only need to once if you have your own account - but if you’re using a shared account like tx160085, you will need a way to keep your computers separate from everyone else’s.

We’ll do this with Projects, which give you a bit of a workspace in which to keep things that belong to “you”.

Click on “Projects” up along the top.

Name the project for yourself, click “create”¶

Enter your name into the Project Name, and something simple like “ANGUS” into the description. Then click ‘create’.

Select the newly created project¶

Click on your new project!

Within the project, select “new”¶

Now, select ‘New’ and then “Instance” from the dropdown menu to start up a new machine.

Find the “Ubuntu 16.04” image, click on it¶

Enter “Ubuntu 16.04” into the search bar - make sure it’s from June 21st, 2017.

Name it something simple and select ‘m1.medium’¶

Change the name after what we’re doing - “workshop tutorial”, for example, but it doesn’t matter – and select ‘m1.medium’.

Wait for it to become active¶

It will now be booting up! This will take 2-10 minutes, depending. Just wait! Don’t reload or anything.

Click on your new instance to get more information!¶

Now, you can either click “Open Web Shell”, or, if you know how to use secure shell (ssh), you can ssh in as user ‘tx160085’ on the IP address of the machine - see circled information below. Click Here to access the tutorial for setting up a ssh connection from your local computer to your Jetstream instance.

Miscellany¶

There’s a possibility that you’ll be confronted with this when you log in to jetstream:

A refresh of the page should get you past it. Please try not to actually move any instances to a new project; it’s probably someone else’s and it could confuse them :)

Suspend your instance¶

You can save your workspace so you can return to your instance at a later time without losing any of your files or information stored in memory, similiar to putting your physical computer to sleep. At the Instance Details screen, select the “Suspend” button.

This will open up a dialogue window. Select the “Yes, suspend this instance” button.

It may take Jetstream a few minutes to process, so wait until the progress bar says “Suspended.”

Shutting down your instance¶

You can shut down your workspace so you can return to your instance another day without losing any of your files, similiar to shutting down your physical computer. You will retain your files, but you will lose any information stored in memory, such as your history on the command line. At the Instance Details screen, select the “Stop” button.

This will open up a dialogue window. Select the “Yes, stop this instance” button.

It may take Jetstream a few minutes to process, so wait until the progress bar says “Shutoff.”

Deleting your instance¶

To completely remove your instance, you can select the “delete” buttom from the instance details page.

This will open up a dialogue window. Select the “Yes, delete this instance” button.

It may take Jetstream a few minutes to process your request. The instance should disappear from the project when it has been successfully deleted.

Learning Objectives¶

• Describe what the shell is and how it is used.
• Summarize reasons why learning the shell is beneficial.

What is the shell?¶

The shell is a program that presents a command line interface which allows you to control your computer using commands entered with a keyboard instead of controlling graphical user interfaces (GUIs) with a mouse/keyboard combination.

There are many reasons to learn about the shell.

• For most bioinformatics tools, you have to use the shell. There is no graphical interface. If you want to work in metagenomics or genomics you’re going to need to use the shell.
• The shell gives you power. The command line gives you the power to do your work more efficiently and more quickly. When you need to do things tens to hundreds of times, knowing how to use the shell is transformative.
• To use remote computers or cloud computing, you need to use the shell.

Unix is user-friendly. It’s just very selective about who its friends are.

How to access the shell¶

The shell is already available on Mac and Linux. For Windows, you’ll have to download a separate program.

Mac¶

On Mac the shell is available through TerminalApplications -> Utilities -> TerminalGo ahead and drag the Terminal application to your Dock for easy access.

Windows¶

For Windows, we’re going to be using gitbash.Download and install gitbash on your computer. Open up the program.

Starting with the shell¶

We will spend most of our time learning about the basics of the shell by manipulating some experimental data.

Now we’re going to download the data for the tutorial. For this you’ll need internet access, because you’re going to get it off the web.

Open the shell.

Enter the command:

git clone https://github.com/edamame-course/edamame-data.git

This command will grab all of the data needed for this workshop from the internet. (We’re not going to talk about git right now, but it’s a tool for doing version control.)

Now let’s go in to that directory

cd
cd edamame-data

The command cd stands for ‘change directory’

In this directory, there should be some things we just downloaded. Let’s check. Type:

ls

ls stands for ‘list’ and it lists the contents of a directory.

There’s a few directories there, but not too many. Let’s go look in the data directory.

cd shell
ls

In there, all mixed up together are files and directories/folders. If we want to know which is which, we can type:

ls -F

Anything with a "/" after it is a directory. Things with a "" after them are programs. If there’s nothing there it’s a file.

You can also use the command ls -l to see whether items in a directory are files or directories. ls -l gives a lot more information too, such as the size of the file

So, we can see that we have several files, directories and a program. Great!

Arguments¶

Most programs take additional arguments that control their exact behavior. For example, -F and -l are arguments to ls. The ls program, like many programs, take a lot of arguments. But how do we know what the options are to particular commands?

Most commonly used shell programs have a manual. You can access the manual using the man program. Try entering:

man ls

This will open the manual page for ls. Use the space key to go forward and b to go backwards. When you are done reading, just hit q to quit.

Programs that are run from the shell can get extremely complicated. To see an example, open up the manual page for the find program. No one can possibly learn all of these arguments, of course. So you will probably find yourself referring back to the manual page frequently.

The Unix directory file structure (a.k.a. where am I?)¶

As you’ve already just seen, you can move around in different directories or folders at the command line. Why would you want to do this, rather than just navigating around the normal way.

When you’re working with bioinformatics programs, you’re working with your data and it’s key to be able to have that data in the right place and make sure the program has access to the data. Many of the problems people run in to with command line bioinformatics programs is not having the data in the place the program expects it to be.

Moving around the file system¶

Let’s practice moving around a bit.

We’re going to work in that shell directory we just downloaded.

First let’s navigate there using the regular way by clicking on the different folders.

First we did something like go to the folder of our username. Then we opened 'edamame-data' then 'shell'

Let’s draw out how that went.

Harriet draw image

Now let’s draw some of the other files and folders we could have clicked on.

This is called a hierarchical file system structure, like an upside down tree with root (/) at the base that looks like this.

Exercise

Now we’re going to try a hunt. Move around in the ‘hidden’ directory and try to find the file 'youfoundit.txt'

Examining the contents of other directories¶

By default, the ls commands lists the contents of the working directory (i.e. the directory you are in). You can always find the directory you are in using the pwd command. However, you can also give ls the names of other directories to view. Navigate to the home directory if you are not already there.

Type:

cd

Then enter the command:

ls edamame-data

This will list the contents of the edamame-data directory without you having to navigate there.

The cd command works in a similar way. Try entering:

cd
cd edamame-data/shell/hidden

and you will jump directly to hidden without having to go through the intermediate directory.

Exercise

Try finding the 'anotherfile.txt' file without changing directories.

cd e<tab>

ls C01<tab><tab>

Saving time with shortcuts, wild cards, and tab completion¶

cd
cd edamame-data
cd shell

ls ~

ls ..

ls ../../

~/edamame-data/shell/MiSeq

ls *fastq

ls /usr/bin/*.sh

ls *R*fastq

Command History¶

You can easily access previous commands. Hit the up arrow. Hit it again. You can step backwards through your command history. The down arrow takes your forwards in the command history.

^-C will cancel the command you are writing, and give you a fresh prompt.

^-R will do a reverse-search through your command history. This is very useful.

You can also review your recent commands with the history command. Just enter:

history

to see a numbered list of recent commands, including this just issues history command. You can reuse one of these commands directly by referring to the number of that command.

If your history looked like this:

259  ls *
260  ls /usr/bin/*.sh
261  ls *R*fastq

then you could repeat command #260 by simply entering:

!260

(that’s an exclamation mark).

Short Exercise

1. Find the line number in your history for the last exercise (listing files in /bin) and reissue that command.

Examining Files¶

We now know how to switch directories, run programs, and look at the contents of directories, but how do we look at the contents of files?

The easiest way to examine a file is to just print out all of the contents using the program cat. Enter the following command:

cat C01D01F_sub.fastq

This prints out the contents of the C01D01F_sub.fastq file.

Short Exercises

1. Print out the contents of the ~/edamame-data/shell/MiSeq/Centralia_mapping_files/Collapsed_Centralia_full_map.txt file. What does this file contain?
2. Without changing directories, (you should still be in edamame-data), use one short command to print the contents of all of the files in the /home/username/edamame-data/shell/MiSeq directory.

Make sure we’re in the right place for the next set of the lessons. We want to be in the shell directory. Check if you’re there with pwd and if not navigate there. One way to do that would be

cd ~/edamame-data/shell/MiSeq

cat is a terrific program, but when the file is really big, it can be annoying to use. The program, less, is useful for this case. Enter the following command:

less C01D01F_sub.fastq

less opens the file, and lets you navigate through it. The commands are identical to the man program.

Some commands in less

| key | action | | ——- | ———- | | “space” | to go forward | | “b” | to go backwarsd | | “g” | to go to the beginning | | “G” | to go to the end | | “q” | to quit |

less also gives you a way of searching through files. Just hit the “/” key to begin a search. Enter the name of the word you would like to search for and hit enter. It will jump to the next location where that word is found. Try searching the dictionary.txt file for the word “cat”. If you hit “/” then “enter”, less will just repeat the previous search. less searches from the current location and works its way forward. If you are at the end of the file and search for the word “cat”, less will not find it. You need to go to the beginning of the file and search.

For instance, let’s search for the sequence 1101:14341 in our file. You can see that we go right to that sequence and can see what it looks like.

Remember, the man program actually uses less internally and therefore uses the same commands, so you can search documentation using “/” as well!

There’s another way that we can look at files, and in this case, just look at part of them. This can be particularly useful if we just want to see the beginning or end of the file, or see how it’s formatted.

The commands are head and tail and they just let you look at the beginning and end of a file respectively.

head C01D01F_sub.fastq
tail C01D01F_sub.fastq

The -n option to either of these commands can be used to print the first or last n lines of a file. To print the first/last line of the file use:

head -n 1 C01D01F_sub.fastq
tail -n 1 C01D01F_sub.fastq

Searching files¶

We showed a little how to search within a file using less. We can also search within files without even opening them, using grep. Grep is a command-line utility for searching plain-text data sets for lines matching a string or regular expression. Let’s give it a try!

Let’s search for that sequence 22029:7208 in the C01D01F_sub.fastq file.

grep 22029:7208 C01D01F_sub.fastq

We get back the whole line that had '22029:7208' in it. What if we wanted all four lines, the whole part of that FASTQ sequence, back instead.

grep -A 3 22029:7208 C01D01F_sub.fastq

The -A flag stands for “after match” so it’s returning the line that matches plus the three after it. The -B flag returns that number of lines before the match.

Exercise

1. Search for the sequence 'TTATCCGGATTTATTGGGTTTAAAGGGT' in the C01D01F_sub.fastq file and in the output have the sequence name and the sequence. e.g.@M00967:43:000000000-A3JHG:1:2114:11799:28499 1:N:0:188 TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGATGCAG
2. Search for that sequence in all the FASTQ files.

Redirection¶

We’re excited we have all these sequences that we care about that we just got from the FASTQ files. That is a really important motif that is going to help us answer our important question. But all those sequences just went whizzing by with grep. How can we capture them?

We can do that with something called “redirection”. The idea is that we’re redirecting the output to the terminal (all the stuff that went whizzing by) to something else. In this case, we want to print it to a file, so that we can look at it later.

The redirection command for putting something in a file is >

Let’s try it out and put all the sequences that contain 'TTATCCGGATTTATTGGGTTTAAAGGGT' from all the files in to another file called 'good-data2.txt'

grep -B 2 TTATCCGGATTTATTGGGTTTAAAGGGT *.fastq > good-data2.txt

The prompt should sit there a little bit, and then it should look like nothing happened. But type ls. You should have a new file called good-data2.txt. Take a look at it and see if it has what you think it should.

There’s one more useful redirection command that we’re going to show, and that’s called the pipe command, and it is |. It’s probably not a key on your keyboard you use very much. What | does is take the output that scrolling by on the terminal and then can run it through another command. When it was all whizzing by before, we wished we could just slow it down and look at it, like we can with less. Well it turns out that we can! We pipe the grep command through less

grep TTATCCGGATTTATTGGGTTTAAAGGGT *.fastq | less

Now we can use the arrows to scroll up and down and use q to get out.

We can also do something tricky and use the command wc. wc stands for word count. It counts the number of lines or characters. So, we can use it to count the number of lines we’re getting back from our grep command. And that will magically tell us how many sequences we’re finding. We’re

grep TTATCCGGATTTATTGGGTTTAAAGGGT *.fastq | wc

That tells us the number of lines, words and characters in the file. If we just want the number of lines, we can use the -l flag for lines.

grep TTATCCGGATTTATTGGGTTTAAAGGGT *.fastq | wc -l

Redirecting is not super intuitive, but it’s really powerful for stringing together these different commands, so you can do whatever you need to do.

The philosophy behind these command line programs is that none of them really do anything all that impressive. BUT when you start chaining them together, you can do some really powerful things really efficiently. If you want to be proficient at using the shell, you must learn to become proficient with the pipe and redirection operators: |, >, >>.

Creating, moving, copying, and removing¶

Now we can move around in the file structure, look at files, search files, redirect. But what if we want to do normal things like copy files or move them around or get rid of them. Sure we could do most of these things without the command line, but what fun would that be?! Besides it’s often faster to do it at the command line, or you’ll be on a remote server like Amazon where you won’t have another option.

The stability.files file is one that tells us what sample name goes with what sequences. This is a really important file, so we want to make a copy so we don’t lose it.

Lets copy the file using the cp command. The cp command backs up the file. Navigate to the data directory and enter:

cp good-data2.txt good-data2.backup.txt

Now good-data2.backup.txt has been created as a copy of good-data2.txt.

Let’s make a backup directory where we can put this file.

The mkdir command is used to make a directory. Just enter mkdir followed by a space, then the directory name.

mkdir backup

We can now move our backed up file in to this directory. We can move files around using the command mv. Enter this command:

mv good-data2.backup.txt backup/

This moves good-data2.backup.txt into the directory backup/ or the full path would be ~/edamame-data/shell/MiSeq/backup

The mv command is also how you rename files. Since this file is so important, let’s rename it:

mv backup/good-data2.backup.txt backup/good-data2.backup_IMPORTANT

Now the file name has been changed to good-data2.backup_IMPORTANT. Let’s delete the backup file now:

rm backup/good-data2.backup_IMPORTANT

The rm file removes the file. Be careful with this command. It doesn’t just nicely put the files in the Trash. They’re really gone.

Short Exercise

Do the following:

1. Rename the good-data2.backup_IMPORTANT file to good-data2.backup.txt.
2. Create a directory in the MiSeq directory called new
3. Then, copy the good-data2.backup.txt file into new

By default, rm, will NOT delete directories. You can tell rm to delete a directory using the -r option. Let’s delete that new directory we just made. Enter the following command: rm -r new

Running programs¶

Commands like ls, rm, echo, and cd are just ordinary programs on the computer. A program is just a file that you can execute. The program which tells you the location of a particular program. For example:

which ls

Will return “/bin/ls”. Thus, we can see that ls is a program that sits inside of the /bin directory. Now enter:

which find

You will see that find is a program that sits inside of the /usr/bin directory.

So ... when we enter a program name, like ls, and hit enter, how does the shell know where to look for that program? How does it know to run /bin/ls when we enter ls. The answer is that when we enter a program name and hit enter, there are a few standard places that the shell automatically looks. If it can’t find the program in any of those places, it will print an error saying “command not found”. Enter the command:

echo $PATH

This will print out the value of the PATH environment variable. More on environment variables later. Notice that a list of directories, separated by colon characters, is listed. These are the places the shell looks for programs to run. If your program is not in this list, then an error is printed. The shell ONLY checks in the places listed in the PATH environment variable.

Navigate to the shell directory and list the contents. You will notice that there is a program (executable file) called hello.sh in this directory. Now, try to run the program by entering:

hello.sh

You should get an error saying that hello.sh cannot be found. That is because the directory /home/username/edamame-data/shell is not in the PATH. You can run the hello.sh program by entering:

./hello.sh

Remember that . is a shortcut for the current working directory. This tells the shell to run the hello.sh program which is located right here. So, you can run any program by entering the path to that program. You can run hello.sh equally well by specifying:

/home/username/edamame-data/shell/hello.sh

Or by entering:

~/edamame-data/shell/hello.sh

When there are no / characters, the shell assumes you want to look in one of the default places for the program.

Where can I learn more about the shell?¶

• Software Carpentry tutorial - The Unix shell
• The shell handout - Command Reference
• explainshell.com
• http://tldp.org/HOWTO/Bash-Prog-Intro-HOWTO.html
• man bash
• Google - if you don’t know how to do something, try Googling it. Other people have probably had the same question.
• Learn by doing. There’s no real other way to learn this than by trying it out. Write your next paper in nano (really emacs or vi), open pdfs from the command line, automate something you don’t really need to automate.

Bonus:¶

backtick, xargs: Example find all files with certain text

alias -> rm -i

variables -> use a path example

.bashrc

du

ln

ssh and scp

Regular Expressions

Permissions

Chaining commands together

Updating the software on the machine¶

Copy and paste the following commands:

sudo apt-get update && sudo apt-get -y install python ncbi-blast+

(make sure to hit enter after the paste – sometimes the last line doesn’t paste completely.)

This updates the software list and installs the Python programming language and NCBI BLAST+.

Running BLAST¶

First! We need some data. Let’s grab the mouse and zebrafish RefSeq protein data sets from NCBI, and put them in our home directory. If you’ve just logged in, you should be there already, but if you’re unsure, just run cd and hit enter. Now, we’ll use curl to download the files:

curl -O ftp://ftp.ncbi.nih.gov/refseq/M_musculus/mRNA_Prot/mouse.1.protein.faa.gz
curl -O ftp://ftp.ncbi.nih.gov/refseq/M_musculus/mRNA_Prot/mouse.2.protein.faa.gz
curl -O ftp://ftp.ncbi.nih.gov/refseq/M_musculus/mRNA_Prot/mouse.3.protein.faa.gz

curl -O ftp://ftp.ncbi.nih.gov/refseq/D_rerio/mRNA_Prot/zebrafish.1.protein.faa.gz

If you look at the files in the current directory, you should see four files, along with a directory called lost+found which is for system information:

ls -l

should show you:

total 29056
drwxr-xr-x 2 titus titus     4096 May  5 08:26 Desktop
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Documents
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Downloads
-rw-rw-r-- 1 titus titus  3610407 Jun 14 18:11 mouse.1.protein.faa.gz
-rw-rw-r-- 1 titus titus  6080985 Jun 14 18:11 mouse.2.protein.faa.gz
-rw-rw-r-- 1 titus titus  7520591 Jun 14 18:11 mouse.3.protein.faa.gz
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Music
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Pictures
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Public
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Templates
drwxr-xr-x 2 titus titus     4096 Jun 14 18:03 Videos
-rw-rw-r-- 1 titus titus 12500932 Jun 14 18:11 zebrafish.1.protein.faa.gz

All four of the files are FASTA protein files (that’s what the .faa suggests) that are compressed with gzip (that’s what the .gz means).

Uncompress them:

gunzip *.faa.gz

and let’s look at the first few sequences in the file:

head mouse.1.protein.faa 

These are protein sequences in FASTA format. FASTA format is something many of you have probably seen in one form or another – it’s pretty ubiquitous. It’s a text file, containing records; each record starts with a line beginning with a ‘>’, and then contains one or more lines of sequence text.

Let’s take those first two sequences and save them to a file. We’ll do this using output redirection with ‘>’, which says “take all the output and put it into this file here.”

head -11 mouse.1.protein.faa > mm-first.fa

So now, for example, you can do cat mm-first.fa to see the contents of that file (or less mm-first.fa).

Now let’s BLAST these two sequences against the entire zebrafish protein data set. First, we need to tell BLAST that the zebrafish sequences are (a) a database, and (b) a protein database. That’s done by calling ‘makeblastdb’:

makeblastdb -in zebrafish.1.protein.faa -dbtype prot

Next, we call BLAST to do the search:

blastp -query mm-first.fa -db zebrafish.1.protein.faa

This should run pretty quickly, but you’re going to get a lot of output!! To save it to a file instead of watching it go past on the screen, ask BLAST to save the output to a file that we’ll name mm-first.x.zebrafish.txt:

blastp -query mm-first.fa -db zebrafish.1.protein.faa -out mm-first.x.zebrafish.txt

and then you can ‘page’ through this file at your leisure by typing:

less mm-first.x.zebrafish.txt

(Type spacebar to move down, and ‘q’ to get out of paging mode.)

Let’s do some more sequences (this one will take a little longer to run):

head -500 mouse.1.protein.faa > mm-second.fa
blastp -query mm-second.fa -db zebrafish.1.protein.faa -out mm-second.x.zebrafish.txt

will compare the first 83 sequences. You can look at the output file with:

less mm-second.x.zebrafish.txt

(and again, type ‘q’ to get out of paging mode.)

To get an output format that reads well into downstream applications, it is helpful to add the flag -outfmt 6

blastp -query mm-second.fa -db zebrafish.1.protein.faa -out mm-second.x.zebrafish.tbl.txt -outfmt 6

To see the results:

head mm-second.x.zebrafish.tbl.txt | less -S

less -S means it is scrollable in screen from left to right. type ‘q’ to exit less

To find out how to customize blast outputs, it is helpful to look at the BLAST® Command Line Applications User Manual.

Notes:

• you can execute multiple commands at a time;

• You might see a warning -

  Selenocysteine (U) at position 310 replaced by X

  what does this mean?

• why did it take longer to BLAST mm-second.fa than mm-first.fa?

Things to mention and discuss:

• blastp options and -help.
• command line options, more generally - why so many?
• automation rocks!

Reminder: shut down your instance!

Other topics to discuss:

• when you shut down, you lose all your data
• what computer(s) is this all happening on?

Installing software for the workshop¶

First, let’s install software for short read quality assessment, trimming and python virtual environments:

sudo apt-get -y update && \
sudo apt-get -y install trimmomatic python-pip \
   samtools zlib1g-dev ncurses-dev python-dev
sudo apt-get install -y python3.5-dev python3.5-venv make \
    libc6-dev g++ zlib1g-dev

 wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
 unzip fastqc_v0.11.5.zip
 cd FastQC
 chmod +x fastqc
 cd

Now, create a python 3.5 virtual environment and install software within:

python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer

pip install -U https://github.com/dib-lab/sourmash/archive/master.zip

Running Jupyter Notebook¶

Let’s also run a Jupyter Notebook. First, configure it a teensy bit more securely, and also have it run in the background:

jupyter notebook --generate-config

cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8000

EOF

Now, run!

jupyter notebook &

On the Jetstream webshell, you can get the Web page address for the jupyter notebook you just launched. Press enter, then execute

echo http://$(hostname):8000/

Copy the output from echo and paste it into the url bar of a new tab in your web browser.

Note, the password is ‘davis’.

ssh -N -f -L localhost:8000:localhost:8000 username@remotehost

Data source¶

We’re going to be using a subset of data from Hu et al., 2016. This paper from the Banfield lab samples some relatively low diversity environments and finds a bunch of nearly complete genomes.

(See DATA.html for a list of the data sets we’re using in this tutorial.)

1. Copying in some data to work with.¶

We’ve loaded subsets of the data onto an Amazon location for you, to make everything faster for today’s work. We’re going to put the files on your computer locally under the directory ~/data:

mkdir ~/data

Next, let’s grab part of the data set:

cd ~/data
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz

Let’s make sure we downloaded all of our data using md5sum.:

md5sum SRR1976948_1.fastq.gz SRR1976948_2.fastq.gz

You should see this:

37bc70919a21fccb134ff2fefcda03ce  SRR1976948_1.fastq.gz
29919864e4650e633cc409688b9748e2  SRR1976948_2.fastq.gz

Now if you type:

ls -l

you should see something like:

total 346936
-rw-rw-r-- 1 ubuntu ubuntu 169620631 Oct 11 23:37 SRR1976948_1.fastq.gz
-rw-rw-r-- 1 ubuntu ubuntu 185636992 Oct 11 23:38 SRR1976948_2.fastq.gz

These are 1m read subsets of the original data, taken from the beginning of the file.

One problem with these files is that they are writeable - by default, UNIX makes things writeable by the file owner. This poses an issue with creating typos or errors in raw data. Let’s fix that before we go on any further:

chmod u-w *

We’ll talk about what these files are below.

1. Copying data into a working location¶

First, make a working directory; this will be a place where you can futz around with a copy of the data without messing up your primary data:

mkdir ~/work
cd ~/work

Now, make a “virtual copy” of the data in your working directory by linking it in –

ln -fs ~/data/* .

These are FASTQ files – let’s take a look at them:

less SRR1976948_1.fastq.gz

(use the spacebar to scroll down, and type ‘q’ to exit ‘less’)

Question:

• where does the filename come from?
• why are there 1 and 2 in the file names?

Links:

• FASTQ Format

2. FastQC¶

We’re going to use FastQC to summarize the data. We already installed ‘fastqc’ on our computer for you.

Now, run FastQC on two files:

fastqc SRR1976948_1.fastq.gz
fastqc SRR1976948_2.fastq.gz

Now type ‘ls’:

ls -d *fastqc.zip*

to list the files, and you should see:

SRR1976948_1_fastqc.zip
SRR1976948_2_fastqc.zip

Inside each of the fatqc directories you will find reports from the fastqc. You can download these files using your Jupyter Notebook console, if you like; or you can look at these copies of them:

• SRR1976948_1_fastqc/fastqc_report.html
• SRR1976948_2_fastqc/fastqc_report.html

Questions:

• What should you pay attention to in the FastQC report?
• Which is “better”, file 1 or file 2? And why?

Links:

• FastQC
• FastQC tutorial video

There are several caveats about FastQC - the main one is that it only calculates certain statistics (like duplicated sequences) for subsets of the data (e.g. duplicate sequences are only analyzed for the first

3. Trimmomatic¶

Now we’re going to do some trimming! We’ll be using Trimmomatic, which (as with fastqc) we’ve already installed via apt-get.

The first thing we’ll need are the adapters to trim off:

curl -O -L http://dib-training.ucdavis.edu.s3.amazonaws.com/mRNAseq-semi-2015-03-04/TruSeq2-PE.fa

Now, to run Trimmomatic:

TrimmomaticPE SRR1976948_1.fastq.gz \
              SRR1976948_2.fastq.gz \
     SRR1976948_1.qc.fq.gz s1_se \
     SRR1976948_2.qc.fq.gz s2_se \
     ILLUMINACLIP:TruSeq2-PE.fa:2:40:15 \
     LEADING:2 TRAILING:2 \
     SLIDINGWINDOW:4:2 \
     MINLEN:25

You should see output that looks like this:

...
Input Read Pairs: 1000000 Both Surviving: 885734 (88.57%) Forward Only Surviving: 114262 (11.43%) Reverse Only Surviving: 4 (0.00%) Dropped: 0 (0.00%)
TrimmomaticPE: Completed successfully

Questions:

• How do you figure out what the parameters mean?
• How do you figure out what parameters to use?
• What adapters do you use?
• What version of Trimmomatic are we using here? (And FastQC?)
• Do you think parameters are different for RNAseq and genomic data sets?
• What’s with these annoyingly long and complicated filenames?
• why are we running R1 and R2 together?

For a discussion of optimal trimming strategies, see MacManes, 2014 – it’s about RNAseq but similar arguments should apply to metagenome assembly.

Links:

• Trimmomatic

4. FastQC again¶

Run FastQC again on the trimmed files:

fastqc SRR1976948_1.qc.fq.gz
fastqc SRR1976948_2.qc.fq.gz

And now view my copies of these files:

• SRR1976948_1.qc_fastqc/fastqc_report.html
• SRR1976948_2.qc_fastqc/fastqc_report.html

Let’s take a look at the output files:

less SRR1976948_1.qc.fq.gz

(again, use spacebar to scroll, ‘q’ to exit less).

5. MultiQC¶

MultiQC aggregates results across many samples into a single report for easy comparison.

Install MultiQC within the py3 environment:

pip install multiqc

Now, run Mulitqc on both the untrimmed and trimmed files within the work directory:

multiqc .

And now you should see output that looks like this:

[INFO   ]         multiqc : This is MultiQC v1.0
[INFO   ]         multiqc : Template    : default
[INFO   ]         multiqc : Searching '.'
Searching 15 files..  [####################################]  100%
[INFO   ]          fastqc : Found 4 reports
[INFO   ]         multiqc : Compressing plot data
[INFO   ]         multiqc : Report      : multiqc_report.html
[INFO   ]         multiqc : Data        : multiqc_data
[INFO   ]         multiqc : MultiQC complete

Now we can view the output file using Jupyter Notebook.

Questions:

• is the quality trimmed data “better” than before?
• Does it matter that you still have adapters!?

Optional: K-mer Spectral Error Trimming

Next: Run the MEGAHIT assembler

While the assembly runs...¶

• Discuss CAMI paper.
• What does, and doesn’t, assemble?
• How good is assembly anyway?

After the assembly is finished¶

Let’s first take a look at the assembly:

    less combined/final.contigs.fa

Now we can run a few stats on our assembly. To do this we will use QUAST:

    cd ~/
    git clone https://github.com/ablab/quast.git -b release_4.5
    export PYTHONPATH=$(pwd)/quast/libs/

Now, run QUAST on the assembly:

    cd ~/assembly
    ~/quast/quast.py combined/final.contigs.fa -o combined-report
    cat combined-report/report.txt

What does this say about our assembly? What do the stats not tell us? For thought and discussion: a few nice slides from Jessica Blanton on assessing assembly quality.

At this point we can do a bunch of things:

• annotate the assembly with Prokka;
• evaluate the assembly’s inclusion of k-mers and reads;
• set up a BLAST database so that we can search it for genes of interest;
• quantify the abundance of the contigs or genes in the assembly, using the original read data set Salmon;
• Identify contigs with similar tetra-mer abundance and coverage with Binning;

Installing Prokka¶

Download and extract the latest version of prokka:

cd ~/
git clone https://github.com/tseemann/prokka.git

We also will need some dependencies such as bioperl:

sudo apt-get -y install bioperl libdatetime-perl libxml-simple-perl libdigest-md5-perl

This may take a little while.

and we need an XML package from perl

sudo bash
export PERL_MM_USE_DEFAULT=1
export PERL_EXTUTILS_AUTOINSTALL="--defaultdeps"
perl -MCPAN -e 'install "XML::Simple"'
exit

Now, you should be able to add Prokka to your $PATH and set up the index for the sequence database:

export PATH=$PATH:$HOME/prokka/bin
prokka --setupdb

To make sure the database loaded directly:

prokka --listdb

You should see something like:

tx160085@js-157-212:~$ prokka --listdb
[17:04:15] Looking for databases in: /home/tx160085/prokka/bin/../db
[17:04:15] * Kingdoms: Archaea Bacteria Mitochondria Viruses
[17:04:15] * Genera: Enterococcus Escherichia Staphylococcus
[17:04:15] * HMMs: HAMAP
[17:04:15] * CMs: Bacteria Viruses

Prokka uses a core set of the Uniprot-DB Kingdom sets against which it blasts your samples. It is possible to search in a more specific dataset, e.g. the genus Enterococcus, by adding a few flags to the command.

–usegenus –genus Enterococcus

Question: What do you think you would do for adding to the default databases?

Prokka should be good to go now– you can check to make sure that all is well by typing prokka. This should print the help screen with all available options. You can find out more about Prokka databases here.

Running Prokka¶

Make a new directory for the annotation:

cd ~/
mkdir annotation
cd annotation

Link the metagenome assembly file into this directory:

ln -fs ~/mapping/subset_assembly.fa .

Now it is time to run Prokka! There are tons of different ways to specialize the running of Prokka. We are going to keep it simple for now, though. It will take a little bit to run.

prokka subset_assembly.fa --outdir prokka_annotation --prefix metagG --metagenome --kingdom Bacteria

Question: Look at the results of the prokka analysis as it prepares your output file. What types of categories are you seeing flash by on the screen?

Don’t worry, the program tends to pause here:

Running: cat prokka_annotation\/sprot\.faa | parallel --gnu --plain -j 6 --block 242000
--recstart '>' --pipe blastp -query --db /home/tx160085/prokka/bin/../db/kingdom/Bacteria/sprot
-evalue 1e-06 -num_threads 1 -num_descriptions 1 -num_alignments 1 -seg no > prokka_annotation\/sprot\.blast 2> /dev/null

This will generate a new folder called prokka_annotation in which will be a series of files, which are detailed here.

In particular, we will be using the *.ffn file to assess the relative read coverage within our metagenomes across the predicted genomic regions.

Question: Take a moment and look inside the output files.:

cd ~/annotation/prokka_annotation
less -S *.fsa

less reminders:

*Press space_bar to page down *Press q to exit the less commands

Questions?¶

• What can I annotate with prokka?
• Alternatives?
• How do I submit my annotated files to Genbank? EBI??
• Why is it called Prokka?

Install Kraken¶

cd ~
git clone https://github.com/DerrickWood/kraken.git
cd ~/kraken
mkdir ~/kraken/bin
./install_kraken.sh ~/kraken/bin
export PATH=$PATH:$HOME/kraken/bin

Install Kraken Mini DB¶

mkdir ~/KRAKEN
cd ~/KRAKEN
wget http://ccb.jhu.edu/software/kraken/dl/minikraken.tgz
tar -xvf minikraken.tgz

Running Kraken¶

cd ~/annotation
mkdir kraken_annotation

kraken --db ~/KRAKEN/minikraken_20141208/ --threads 2 --fasta-input subset_assembly.fa --output kraken_annotation/subset_assembly.kraken


kraken-translate --db ~/KRAKEN/minikraken_20141208/  kraken_annotation/subset_assembly.kraken > kraken_annotation/subset_assembly.kraken.labels

Kraken has now provided a taxonomic assignment to all of the clusters.

To generate a summary table:

cd ~/annotation
kraken-report --db ~/KRAKEN/minikraken_20141208 kraken_annotation/subset_assembly.kraken > kraken_annotation/subset_assembly.kraken.report

The top of the file lists all the unclassified sequences, to look at the file and skip over these, do the following:

grep -v ^U ~/annotation/kraken_annotation/subset_assembly.kraken.report | head -n20

The output of kraken-report is tab-delimited, with one line per taxon. The fields of the output, from left-to-right, are as follows:

1. Percentage of reads covered by the clade rooted at this taxon
2. Number of reads covered by the clade rooted at this taxon
3. Number of reads assigned directly to this taxon
4. A rank code, indicating (U)nclassified, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. All other ranks are simply ‘-‘.
5. NCBI taxonomy ID
6. indented scientific name

Example output:

tx160085@js-157-212:~/annotation/kraken_annotation$ grep -v ^U subset_assembly.kraken.report | head -n20
 89.60  8311    8311    U       0       unclassified
 10.40  965     0       -       1       root
 10.40  965     3       -       131567    cellular organisms
 10.37  962     43      D       2           Bacteria
  6.51  604     0       P       200918        Thermotogae
  6.51  604     0       C       188708          Thermotogae
  6.51  604     0       O       2419              Thermotogales
  6.51  604     8       F       188709              Thermotogaceae
  5.11  474     0       G       28236                 Petrotoga
  5.11  474     0       S       69499                   Petrotoga mobilis
  5.11  474     474     -       403833                    Petrotoga mobilis SJ95
  1.22  113     0       G       1184396               Mesotoga
  1.22  113     0       S       1184387                 Mesotoga prima
  1.22  113     113     -       660470                    Mesotoga prima MesG1.Ag.4.2
  0.04  4       0       G       651456                Kosmotoga
  0.04  4       0       S       651457                  Kosmotoga olearia
  0.04  4       4       -       521045                    Kosmotoga olearia TBF 19.5.1
  0.02  2       1       G       2335                  Thermotoga
  0.01  1       0       S       177758                  Thermotoga lettingae
  0.01  1       1       -       416591                    Thermotoga lettingae TMO

Why use Kraken?

For a simulated metagenome of 100 bp reads in its fastest mode of operation, , Kraken processed over 4 million reads per minute on a single core, over 900 times faster than Megablast and over 11 times faster than the abundance estimation program MetaPhlAn. Kraken’s accuracy is comparable with Megablast, with slightly lower sensitivity and very high precision. Citation

However, kraken is only as sensitive as the provided database, so for unusual samples, a custom database needs to be constructed . The accuracy is very sensitive to the quantity of samples in the database.

Install Prodigal¶

cd ~
wget https://github.com/hyattpd/Prodigal/releases/download/v2.6.3/prodigal.linux
tar -xvf v2.6.3.tar.gz
chmod 775 ~/prodigal.linux

Running Prodigal¶

Using prodigal with the same set of data, we can get a list of predicted genes.

cd ~/annotation
mkdir prodigal_annotation
~/prodigal.linux -p meta -a prodigal_annotation/subset_assembly.faa -d prodigal_annotation/subset_assembly.fna -f gff -o prodigal_annotation/subset_assembly.gff -i subset_assembly.fa

References¶

• http://www.vicbioinformatics.com/software.prokka.shtml
• https://www.ncbi.nlm.nih.gov/pubmed/24642063
• https://github.com/tseemann/prokka/blob/master/README.md
• http://denbi-metagenomics-workshop.readthedocs.io/en/latest/classification/kraken.html

Objectives¶

• Compare reads to assemblies
• Create and search a custom database
• Compare datasets and build a tree
• Determine what’s in a metagenome (Taxonomic classification)

At the beginning¶

Create / log into an m1.medium Jetstream instance, and run these two commands:

cd
curl -O https://s3-us-west-1.amazonaws.com/spacegraphcats.ucdavis.edu/microbe-genbank-sbt-k31-2017.05.09.tar.gz
tar xzf microbe-genbank-sbt-k31-2017.05.09.tar.gz

– they take a long time :).

K-mers!¶

K-mers are a fairly simple concept that turn out to be tremendously powerful.

A “k-mer” is a word of DNA that is k long:

ATTG - a 4-mer
ATGGAC - a 6-mer

Typically we extract k-mers from genomic assemblies or read data sets by running a k-length window across all of the reads and sequences – e.g. given a sequence of length 16, you could extract 11 k-mers of length six from it like so:

AGGATGAGACAGATAG

becomes the following set of 6-mers:

AGGATG
 GGATGA
  GATGAG
   ATGAGA
    TGAGAC
     GAGACA
      AGACAG
       GACAGA
        ACAGAT
         CAGATA
          AGATAG

k-mers are most useful when they’re long, because then they’re specific. That is, if you have a 31-mer taken from a human genome, it’s pretty unlikely that another genome has that exact 31-mer in it. (You can calculate the probability if you assume genomes are random: there are 4³¹ possible 31-mers, and 4³¹ = 4,611,686,018,427,387,904. So, you know, a lot.)

The important concept here is that long k-mers are species specific. We’ll go into a bit more detail later.

K-mers and assembly graphs¶

We’ve already run into k-mers before, as it turns out - when we were doing genome assembly. One of the three major ways that genome assembly works is by taking reads, breaking them into k-mers, and then “walking” from one k-mer to the next to bridge between reads. To see how this works, let’s take the 16-base sequence above, and add another overlapping sequence:

AGGATGAGACAGATAG
    TGAGACAGATAGGATTGC

One way to assemble these together is to break them down into k-mers –

becomes the following set of 6-mers:

AGGATG
 GGATGA
  GATGAG
   ATGAGA
    TGAGAC
     GAGACA
      AGACAG
       GACAGA
        ACAGAT
         CAGATA
          AGATAG -> off the end of the first sequence
           GATAGG <- beginning of the second sequence
            ATAGGA
             TAGGAT
              AGGATT
               GGATTG
                GATTGC

and if you walk from one 6-mer to the next based on 5-mer overlap, you get the assembled sequence:

AGGATGAGACAGATAGGATTGC

Graphs of many k-mers together are called De Bruijn graphs, and assemblers like MEGAHIT and SOAPdenovo are De Bruijn graph assemblers - they use k-mers underneath.

Why k-mers, though? Why not just work with the full read sequences?¶

Computers love k-mers because there’s no ambiguity in matching them. You either have an exact match, or you don’t. And computers love that sort of thing!

Basically, it’s really easy for a computer to tell if two reads share a k-mer, and it’s pretty easy for a computer to store all the k-mers that it sees in a pile of reads or in a genome.

Long k-mers are species specific¶

So, we’ve said long k-mers (say, k=31 or longer) are pretty species specific. Is that really true?

Yes! Check out this figure from the MetaPalette paper:

here, the Koslicki and Falush show that k-mer similarity works to group microbes by genus, at k=40. If you go longer (say k=50) then you get only very little similarity between different species.

Using k-mers to compare samples against each other¶

So, one thing you can do is use k-mers to compare genomes to genomes, or read data sets to read data sets: data sets that have a lot of similarity probably are similar or even the same genome.

One metric you can use for this comparisons is the Jaccard distance, which is calculated by asking how many k-mers are shared between two samples vs how many k-mers in total are in the combined samples.

only k-mers in both samples
----------------------------
all k-mers in either or both samples

A Jaccard distance of 1 means the samples are identical; a Jaccard distance of 0 means the samples are completely different.

This is a great measure and it can be used to search databases and cluster unknown genomes and all sorts of other things! The only real problem with it is that there are a lot of k-mers in a genome – a 5 Mbp genome (like E. coli) has 5 m k-mers!

About a year ago, Ondov et al. (2016) showed that MinHash approaches could be used to estimate Jaccard distance using only a small fraction (1 in 10,000 or so) of all the k-mers.

The basic idea behind MinHash is that you pick a small subset of k-mers to look at, and you use those as a proxy for all the k-mers. The trick is that you pick the k-mers randomly but consistently: so if a chosen k-mer is present in two data sets of interest, it will be picked in both. This is done using a clever trick that we can try to explain to you in class - but either way, trust us, it works!

We have implemented a MinHash approach in our sourmash software, which can do some nice things with samples. We’ll show you some of these things next!

Installing sourmash¶

To install sourmash, run:

sudo apt-get -y update && \
sudo apt-get install -y python3.5-dev python3.5-venv make \
    libc6-dev g++ zlib1g-dev

this installs Python 3.5.

Now, create a local software install and populate it with Jupyter and other dependencies:

python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer

pip install -U https://github.com/dib-lab/sourmash/archive/master.zip

Generate a signature for Illumina reads¶

Download some reads and a reference genome:

mkdir ~/data
cd ~/data
wget http://public.ged.msu.edu.s3.amazonaws.com/ecoli_ref-5m-trim.pe.fq.gz
wget https://s3.amazonaws.com/public.ged.msu.edu/ecoliMG1655.fa.gz

Compute a scaled MinHash signature from our reads:

mkdir ~/sourmash
cd ~/sourmash

sourmash compute --scaled 10000 ~/data/ecoli_ref*pe*.fq.gz -o ecoli-reads.sig -k 31

Compare reads to assemblies¶

Use case: how much of the read content is contained in the reference genome?

Build a signature for an E. coli genome:

sourmash compute --scaled 10000 -k 31 ~/data/ecoliMG1655.fa.gz -o ecoli-genome.sig

and now evaluate containment, that is, what fraction of the read content is contained in the genome:

sourmash search -k 31 ecoli-reads.sig ecoli-genome.sig --containment

and you should see:

# running sourmash subcommand: search
loaded query: /home/ubuntu/data/ecoli_ref-5m... (k=31, DNA)
loaded 1 signatures from ecoli-genome.sig
1 matches:
similarity   match
----------   -----
 46.6%       /home/ubuntu/data/ecoliMG1655.fa.gz

Why are only 50% or so of our k-mers from the reads in the genome!? Any ideas?

Try the reverse - why is it bigger?

sourmash search -k 31 ecoli-genome.sig ecoli-reads.sig --containment

(...but 99% of our k-mers from the genome are in the reads!?)

This is basically because of sequencing error! Illumina data contains a lot of errors, and the assembler doesn’t include them in the assembly!

Make and search a database quickly.¶

Suppose that we have a collection of signatures (made with sourmash compute as above) and we want to search it with our newly assembled genome (or the reads, even!). How would we do that?

Let’s grab a sample collection of 50 E. coli genomes and unpack it –

mkdir ecoli_many_sigs
cd ecoli_many_sigs

curl -O -L https://github.com/dib-lab/sourmash/raw/master/data/eschericia-sigs.tar.gz

tar xzf eschericia-sigs.tar.gz
rm eschericia-sigs.tar.gz

cd ../

This will produce 50 files named ecoli-N.sig in the ecoli_many_sigs –

ls ecoli_many_sigs

Let’s turn this into an easily-searchable database with sourmash index –

sourmash index -k 31 ecolidb ecoli_many_sigs/*.sig

What does the database look like and how does the search work?

sourmash search (and gather) compare the query and the database (sequence bloom tree) and report matches based containment.

One point to make with this is that the search can quickly narrow down which signatures match your query, without losing any matches. It’s a clever example of how computer scientists can actually make life better :).

And now we can search!

sourmash search ecoli-genome.sig ecolidb.sbt.json -n 20

You should see output like this:

# running sourmash subcommand: search
select query k=31 automatically.
loaded query: /home/tx160085/data/ecoliMG165... (k=31, DNA)
loaded SBT ecolidb.sbt.json
Searching SBT ecolidb.sbt.json
49 matches; showing first 20:
similarity   match
----------   -----
 75.4%       NZ_JMGW01000001.1 Escherichia coli 1-176-05_S4_C2 e117605...
 72.2%       NZ_GG774190.1 Escherichia coli MS 196-1 Scfld2538, whole ...
 71.4%       NZ_JMGU01000001.1 Escherichia coli 2-011-08_S3_C2 e201108...
 70.1%       NZ_JHRU01000001.1 Escherichia coli strain 100854 100854_1...
 69.0%       NZ_JH659569.1 Escherichia coli M919 supercont2.1, whole g...
 64.9%       NZ_JNLZ01000001.1 Escherichia coli 3-105-05_S1_C1 e310505...
 63.0%       NZ_MOJK01000001.1 Escherichia coli strain 469 Cleandata-B...
 62.9%       NZ_MOGK01000001.1 Escherichia coli strain 676 BN4_676_1_(...
 62.0%       NZ_JHDG01000001.1 Escherichia coli 1-176-05_S3_C1 e117605...
 59.9%       NZ_MIWF01000001.1 Escherichia coli strain AF7759-1 contig...
 52.7%       NZ_KE700241.1 Escherichia coli HVH 147 (4-5893887) acYxy-...
 51.7%       NZ_APWY01000001.1 Escherichia coli 178200 gec178200.conti...
 49.3%       NZ_LVOV01000001.1 Escherichia coli strain swine72 swine72...
 49.3%       NZ_MIWP01000001.1 Escherichia coli strain K6412 contig_00...
 49.0%       NZ_LQWB01000001.1 Escherichia coli strain GN03624 GCID_EC...
 48.9%       NZ_JHGJ01000001.1 Escherichia coli O45:H2 str. 2009C-4780...
 48.1%       NZ_CP011331.1 Escherichia coli O104:H4 str. C227-11, comp...
 47.7%       NZ_JHNB01000001.1 Escherichia coli O103:H25 str. 2010C-45...
 47.7%       NZ_JHRE01000001.1 Escherichia coli strain 302014 302014_1...
 47.6%       NZ_JHHE01000001.1 Escherichia coli O103:H2 str. 2009C-327...

identifying what genome is in the signature.

Compare many signatures and build a tree.¶

Adjust plotting (this is a bug in sourmash :) –

echo 'backend : Agg' > matplotlibrc

Compare all the signatures:

sourmash compare ecoli_many_sigs/* -o ecoli_cmp

and then plot:

sourmash plot --labels ecoli_cmp

which will produce a file ecoli_cmp.matrix.png and ecoli_cmp.dendro.png which you can then download via FileZilla and view on your local computer.

Now open jupyter notebook and visualize:

from IPython.display import Image
Image("ecoli_cmp.matrix.png")

Here’s a PNG version:

What’s in my metagenome?¶

At the beginning, we downloaded and unpacked a GenBank index of all the microbial genomes – you can see a basic description here, CTB’s blog post – this one contains sketches of all 100k Genbank microbes. (See available sourmash databases for more information.)

After this database is unpacked, it produces a file genbank-k31.sbt.json and a whole bunch of hidden files in the directory .sbt.genbank-k31.

Next, run the ‘gather’ command to see what’s in your ecoli genome –

at the command line run:

sourmash gather -k 31 ecoli-reads.sig ../genbank-k31.sbt.json -o ecoli-reads-gather-reults.csv

and you should get:

# running sourmash subcommand: sbt_gather
loaded query: /home/tx160085/data/ecoli_ref-... (k=31, DNA)
loaded SBT ../genbank-k31.sbt.json

overlap     p_query p_match
---------   ------- --------
4.9 Mbp      46.6%  100.0%      APIN01000001.1 Escherichia coli str. ...
3.4 Mbp      32.4%    1.5%      AQAX01000001.1 Escherichia coli P0304...
2.5 Mbp      23.6%    0.7%      JHDK01000001.1 Escherichia coli 1-110...
160.0 kbp     1.5%    0.6%      AMPE01000001.1 Citrobacter sp. L17 co...
0.5 Mbp       4.8%    0.4%      BBVS01000001.1 Escherichia albertii D...
3.3 Mbp      31.2%    0.4%      APYP01000001.1 Escherichia coli P0299...
3.6 Mbp      34.5%    0.4%      AQCO01000001.1 Escherichia coli MP021...
4.3 Mbp      40.9%    0.4%      JONB01000001.1 Escherichia coli 3-020...
2.8 Mbp      26.0%    0.5%      AZPL01000001.1 Shigella flexneri 2001...
2.8 Mbp      26.0%    0.2%      MIIY01000001.1 Shigella sp. FC1737 NO...
2.2 Mbp      20.7%    0.2%      AFET01000004.1 Escherichia coli AA86 ...
found less than 10.0 kbp in common. => exiting

found 11 matches total;
the recovered matches hit 49.2% of the query

You can use this on metagenomes (assembled and unassembled) as well; you’ve just got to make the signature files.

To see this in action, here is gather running on a signature generated from some sequences that assemble (but don’t align to known genomes) from the Shakya et al. 2013 mock metagenome paper.

wget https://github.com/dib-lab/sourmash/raw/master/doc/_static/shakya-unaligned-contigs.sig
sourmash gather -k 31 shakya-unaligned-contigs.sig ../genbank-k31.sbt.json

This should yield:

# running sourmash subcommand: gather
loaded query: mqc500.QC.AMBIGUOUS.99.unalign... (k=31, DNA)
loaded SBT genbank-k31.sbt.json

overlap     p_query p_match
---------   ------- --------
1.4 Mbp      11.0%   58.0%      JANA01000001.1 Fusobacterium sp. OBRC1 c
1.0 Mbp       7.7%   25.9%      CP001957.1 Haloferax volcanii DS2 plasmi
0.9 Mbp       7.5%   11.8%      BA000019.2 Nostoc sp. PCC 7120 DNA, comp
0.7 Mbp       5.9%   23.0%      FOVK01000036.1 Proteiniclasticum ruminis
0.7 Mbp       5.3%   17.6%      AE017285.1 Desulfovibrio vulgaris subsp.
0.6 Mbp       4.9%   11.1%      CP001252.1 Shewanella baltica OS223, com
0.6 Mbp       4.8%   27.3%      AP008226.1 Thermus thermophilus HB8 geno
0.6 Mbp       4.4%   11.2%      CP000031.2 Ruegeria pomeroyi DSS-3, comp
480.0 kbp     3.8%    7.6%      CP000875.1 Herpetosiphon aurantiacus DSM
410.0 kbp     3.3%   10.5%      CH959317.1 Sulfitobacter sp. NAS-14.1 sc
1.4 Mbp      10.9%   11.8%      LN831027.1 Fusobacterium nucleatum subsp
0.5 Mbp       4.1%    5.3%      CP000753.1 Shewanella baltica OS185, com
420.0 kbp     3.3%    7.7%      FNDZ01000023.1 Proteiniclasticum ruminis
150.0 kbp     1.2%    4.5%      CP015081.1 Deinococcus radiodurans R1 ch
150.0 kbp     1.2%    8.2%      CP000969.1 Thermotoga sp. RQ2, complete
290.0 kbp     2.3%    4.1%      CH959311.1 Sulfitobacter sp. EE-36 scf_1
1.2 Mbp       9.4%    5.0%      CP013328.1 Fusobacterium nucleatum subsp
110.0 kbp     0.9%    3.5%      FREL01000833.1 Enterococcus faecalis iso
0.6 Mbp       5.0%    2.8%      CP000527.1 Desulfovibrio vulgaris DP4, c
340.0 kbp     2.7%    3.3%      KQ235732.1 Fusobacterium nucleatum subsp
70.0 kbp      0.6%    1.2%      CP000850.1 Salinispora arenicola CNS-205
60.0 kbp      0.5%    0.7%      CP000270.1 Burkholderia xenovorans LB400
50.0 kbp      0.4%    2.6%      CP001080.1 Sulfurihydrogenibium sp. YO3A
50.0 kbp      0.4%    3.2%      L77117.1 Methanocaldococcus jannaschii D
found less than 40.0 kbp in common. => exiting

found 24 matches total;
the recovered matches hit 73.4% of the query

In our recent preprint using this, we showed that

It is straightforward to build your own databases for use with search and gather; this is of interest if you have dozens or hundreds of sequencing data sets in your group. Ping us if you want us to write that up.

Final thoughts on sourmash¶

There are many tools like Kraken and Kaiju that can do taxonomic classification of individual reads from metagenomes; these seem to perform well (albeit with high false positive rates) in situations where you don’t necessarily have the genome sequences that are in the metagenome. Sourmash, by contrast, can estimate which known genomes are actually present, so that you can extract them and map/align to them. It seems to have a very low false positive rate and is quite sensitive to strains.

Above, we’ve shown you a few things that you can use sourmash for. Here is a (non-exclusive) list of other uses that we’ve been thinking about –

• detect contamination in sequencing data;
• index and search private sequencing collections;
• search all of SRA for overlaps in metagenomes;

Chat with Luiz, Phil, or Titus if you are interested in these use cases!

Installing binners¶

MaxBin

 cd
 curl  https://downloads.jbei.org/data/microbial_communities/MaxBin/getfile.php?MaxBin-2.2.2.tar.gz > MaxBin-2.2.2.tar.gz
 tar xzvf MaxBin-2.2.2.tar.gz
 cd MaxBin-2.2.2
 cd src
 make
 cd ..
 ./autobuild_auxiliary
 export PATH=$PATH:~/MaxBin-2.2.2

MetaBAT

cd 
curl -L https://bitbucket.org/berkeleylab/metabat/downloads/metabat-static-binary-linux-x64_v0.32.4.tar.gz > metabatv0.32.4.tar.gz
tar xvf metabatv0.32.4.tar.gz

Time to finally run the Binners! Note: MaxBin can take a lot of time to run and bin your metagenome. As this is a workshop, we are doing three things that sacrifice quality for speed.

1. We are only using 2 of the 6 datasets that were generated for the this project. Most binning software, rely upon many samples to accurately bin data. And, we have subsampled the data to make it faster to proess.
2. We did quick and dirty contig estimation using samtools. If you use MaxBin, I would recommend allowing them to generate the abudance tables using Bowtie2. But, again, this will add to the run time.
3. We are limiting the number of iterations that are performed through the MaxBin expectation-maximization algorithm (5 iterations instead of 50+). This will likely limit the quality of the bins we get out. So, users beware and read the user’s manual before proceeding with your own data analysis.

Get a rough count of the total number of reads mapping¶

Q. Why is the read counts per contigs not an accurate number?

cd ~/mapping

for i in *sorted.bam
  do
    samtools idxstats $i > $i.idxstats
    cut -f1,3 $i.idxstats > $i.counts
  done

Binning 1 - MaxBin¶

–

Maxbin uses read coverage & tetranucleotide frequencies for each contig, and marker gene counts for each bin

First, we will get a list of the count files that we have to pass to MaxBin

mkdir ~/mapping/binning
cd ~/mapping/binning
ls ../*counts > abundance.list

Now, on to the actual binning

run_MaxBin.pl -contig ../subset_assembly.fa -abund_list abundance.list -max_iteration 5 -out mbin

This will generate a series of files. Take a look at the files generated. In particular you should see a series of *.fasta files preceeded by numbers. These are the different genome bins predicted by MaxBin.

Take a look at the mbin.summary file. What is shown?

Now, we are going to generate a concatenated file that contains all of our genome bins put together. We will change the fasta header name to include the bin number so that we can tell them apart later.

for file in mbin.*.fasta
  do 
    num=${file//[!0-9]/}
    sed -e "/^>/ s/$/ ${num}/" mbin.$num.fasta >> maxbin_binned.concat.fasta
  done

And finally make an annotation file for visualization

echo label > maxbin_annotation.list
grep ">" maxbin_binned.concat.fasta |cut -f2 -d ' '>> maxbin_annotation.list

Binning 2 - MetaBAT¶

–

MetaBAT uses read coverage, coverage variance, & tetranucleotide frequencies for each contig. This is done with a custom script

cd ~/mapping
~/metabat/jgi_summarize_bam_contig_depths --outputDepth depth_var.txt *sorted.bam

Setup another “binning” directory for this tool’s results

mkdir ~/mapping/binning_metabat
cd ~/mapping/binning_metabat

Run MetaBAT script

Note that we are outputting info to a logfile

~/metabat/metabat -i ../subset_assembly.fa -a ../depth_var.txt --verysensitive -o metabat -v > log.txt

Make the .fasta file of all binned sequences

for file in metabat.*.fa
  do
    num=${file//[!0-9]/} 
   sed -e "/^>/ s/$/ ${num}/" metabat.$num.fa >> metabat_binned.concat.fasta
done 

Make an annotation file of the bin numbers for annotation in VizBin

echo label > metabat_annotation.list
grep ">" metabat_binned.concat.fasta |cut -f2 -d ' '>> metabat_annotation.list

Visualizing the bins¶

Now that we have our binned data from both MetaBAT and MaxBin there are several different things we can do. One thing we might want to do is check the quality of the binning– a useful tool for this is CheckM. Today, for the sake of time, we will visualize the bins that we just generated using VizBin.

First, install VizBin::

cd
sudo apt-get install libatlas3-base libopenblas-base default-jre
curl -L https://github.com/claczny/VizBin/blob/master/VizBin-dist.jar?raw=true > VizBin-dist.jar

VizBin can run in OSX or Linux but is very hard to install on Windows. To simplify things we are going to run VizBin in the desktop emulator through JetStream (which is ... a bit clunky). So, go back to the Jetstream and open up the web desktop simulator.

Open the terminal through the desktop simulator and open VizBin:

java -jar VizBin-dist.jar

This should prompt VizBin to open in another window. First we will look at the output of the MaxBin assembly. Click the choose button to open file browser to navigate to the binning folder (~/mapping/binning). There you will find the concatenated binned fasta file (maxbin_binned.concat.fasta). Upload this file and hit run.

What do you see? Read up a bit on VizBin to see how the visualization is generated.

Now, upload the maxbin_annotation.list file as an annotation file to VizBin. The annotation file contains the bin id for each of the contigs in the assembly that were binned.

Now, do the same for MetaBat!

Compare the results of the two binning methods-

• How many bins were found?
• How distinct are the bins?

Installing Salmon¶

Download and extract the latest version of Salmon and add it to your PATH:

cd
. ~/py3/bin/activate
pip install pandas
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.7.2/Salmon-0.7.2_linux_x86_64.tar.gz
tar -xvzf Salmon-0.7.2_linux_x86_64.tar.gz
cd Salmon-0.7.2_linux_x86_64
export PATH=$PATH:$HOME/Salmon-0.7.2_linux_x86_64/bin

Running Salmon¶

Make a new directory for the quantification of data with Salmon:

mkdir ~/quant
cd ~/quant

Grab the nucleotide (*ffn) predicted protein regions from Prokka and link them here. Also grab the trimmed sequence data (*fq)

ln -fs ~/annotation/prokka_annotation/metagG.ffn .
ln -fs ~/annotation/prokka_annotation/metagG.gff .
ln -fs ~/data/*.abundtrim.subset.pe.fq.gz .

Create the salmon index:

salmon index -t metagG.ffn -i transcript_index --type quasi -k 31

Salmon requires that paired reads be separated into two files. We can split the reads using the split-paired-reads.py from the khmer package:

pip install khmer

for file in *.abundtrim.subset.pe.fq.gz
do
  tail=.fq.gz
  BASE=${file/$tail/}
  split-paired-reads.py $BASE$tail -1 ${file/$tail/}.1.fq -2 ${file/$tail/}.2.fq
done

Now, we can quantify our reads against this reference:

for file in *.pe.1.fq
do
tail1=.abundtrim.subset.pe.1.fq
tail2=.abundtrim.subset.pe.2.fq
BASE=${file/$tail1/}
salmon quant -i transcript_index --libType IU \
      -1 $BASE$tail1 -2 $BASE$tail2 -o $BASE.quant;
 done

(Note that –libType must come before the read files!)

This will create a bunch of directories named after the fastq files that we just pushed through. Take a look at what files there are within one of these directories:

find SRR1976948.quant -type f

Working with count data¶

Now, the quant.sf files actually contain the relevant information about expression – take a look:

head -10 SRR1976948.quant/quant.sf

The first column contains the transcript names, and the fourth column is what we will want down the road - the normalized counts (TPM). However, they’re not in a convenient location / format for use; let’s fix that.

Download the gather-counts.py script:

curl -L -O https://raw.githubusercontent.com/ngs-docs/2016-metagenomics-sio/master/gather-counts.py

and run it:

python2 ./gather-counts.py

This will give you a bunch of .counts files, which are processed from the quant.sf files and named for the directory from which they emanate.

Now, we can create one file:

for file in *counts
do
   name=${file%%.*}
   sed -e "s/count/$name/g" $file > tmp
   mv tmp $file
done
paste *counts |cut -f 1,2,4 > Combined-counts.tab

Plotting the results¶

In Jupyter Notebook, open a new Python3 notebook and name it. Then, into the first cell enter:

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

In another cell (to make sure we are in the correct directory) enter:

cd ~/quant

Now, we can read our data into a pandas dataframe. This is similar to dataframes in R or an Excel spreadsheet. Paste the following into a new cell.:

count_df=pd.read_table('Combined-counts.tab', index_col='transcript')
count_df

And, finally we can plot it!:

fig, ax = plt.subplots(1) #set up a figure and axis handle
count_df.plot(kind='scatter', x='SRR1976948', y='SRR1977249', ax=ax) #plot scatter plot

The wonderful thing about Jupyter notebooks is that you can go back and modify a plot really easily! Let’s try a few modifications with the above plot.

References¶

• http://salmon.readthedocs.io/en/latest/salmon.html
• http://biorxiv.org/content/early/2016/08/30/021592

Downloading data¶

Now, go to a new directory and grab the data:

mkdir ~/mapping
cd ~/mapping

curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249.abundtrim.subset.pe.fq.gz

And extract the files:

for file in *fq.gz
  do
  gunzip $file
done

We will also need the assembly; rather than rebuilding it, you can download a copy that we saved for you:

curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz
gunzip subset_assembly.fa

Mapping the reads¶

First, we will need to to index the megahit assembly:

bwa index subset_assembly.fa

to The reads are in paired-end/interleaved format, so you’ll need to add the -p flag to indicate to bwa that these are paired end data:

Map the reads:

for i in *fq
  do
  bwa mem -p subset_assembly.fa $i > ${i}.aln.sam
done

Converting to BAM to visualize¶

First, index the assembly for samtools:

samtools faidx subset_assembly.fa

Then, convert both SAM files to BAM files:

for i in *.sam
do
   samtools import subset_assembly.fa $i $i.bam
   samtools sort $i.bam $i.bam.sorted
   samtools index $i.bam.sorted.bam
done

Visualizing the read mapping¶

Find a contig name to visualize:

grep -v ^@ SRR1976948.abundtrim.subset.pe.fq.aln.sam | \
    cut -f 3 | sort | uniq -c | sort -n | tail

Pick one e.g. k99_13588.

Now execute:

samtools tview SRR1976948.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400

(use arrow keys to scroll, ‘q’ to quit; a key for what you are looking at: `pileup format<https://en.wikipedia.org/wiki/Pileup_format>`__.)

Look at it in both mappings:

samtools tview SRR1977249.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400

Why is the mapping so good??

Note: no strain variation :).

Grab some untrimmed data:

curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz

Now align this untrimmed data:

gunzip -c SRR1976948_1.fastq.gz | head -800000 > SRR1976948.1
gunzip -c SRR1976948_2.fastq.gz | head -800000 > SRR1976948.2

bwa aln subset_assembly.fa SRR1976948.1 > SRR1976948_1.untrimmed.sai
bwa aln subset_assembly.fa SRR1976948.2 > SRR1976948_2.untrimmed.sai

bwa sampe subset_assembly.fa SRR1976948_1.untrimmed.sai SRR1976948_2.untrimmed.sai SRR1976948.1 SRR1976948.2 > SRR1976948.untrimmed.sam

i=SRR1976948.untrimmed.sam
samtools import subset_assembly.fa $i $i.bam
samtools sort $i.bam $i.bam.sorted
samtools index $i.bam.sorted.bam

And now look:

samtools tview SRR1976948.untrimmed.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:500

You can also use ‘Tablet’ to view the downloaded BAM file - see the Tablet paper.

Assemble the slice¶

(Re)install megahit:

cd
git clone https://github.com/voutcn/megahit.git
cd megahit
make

Go back to the slice directory and extract paired ends:

cd /mnt/slice
extract-paired-ends.py slice.fq

Assemble!

~/megahit/megahit --12 slice.fq.pe -o slice

The contigs will be in slice/final.contigs.fa.

Installing anvi’o (and a few other programs)¶

The first thing we need to do is install anvi’o. To install anvi’o we will be be using Anaconda.

cd ~
wget https://repo.continuum.io/archive/Anaconda3-4.4.0-Linux-x86_64.sh
bash Anaconda3-4.4.0-Linux-x86_64.sh

Now, follow the prompts for the Anaconda installation. To finish the install it will ask you if you would like to add anaconda to the $PATH in you .bashrc. You should say yes. Now, you just need to source your .bashrc to make sure you can use conda.

source .bashrc

Anaconda should now be installed. We will now use anaconda (conda) to install anvi’o (and all its dependencies) as well as source an environment in which to to run conda.

Now, install anvi’o using conda, create an environment in which to run it, and source the environment:

conda create -n anvio232 -c bioconda -c conda-forge gsl anvio
source activate anvio232

Anvi’o should now be installed. But, let’s double check that it worked. They have a nice little test case to check that everything is working well as follows:

anvi-self-test --suite mini

This prompt will start anvi’o processing and ultimately it will generate an interactive window with the anvi’o environment. This is accessible through port 8080 (typically, though it might create go to a different port that will be specified) at your IP address.

To find your IP address you can look back at the Jetstream instance page and copy the IP address.

Now, open a new tab in your browser* and paste in the following:

*Note: This only works in Google Chrome.

[Your IP Address]:8080

This should open up the anvi’o interface which is interactive and pretty good looking.

Now, we just need to install a few other programs, namely, samtools and Bowtie2, which we will use for mapping and looking at our mapped data.

wget https://downloads.sourceforge.net/project/bowtie-bio/bowtie2/2.3.2/bowtie2-2.3.2-linux-x86_64.zip
unzip bowtie2-2.3.2-linux-x86_64.zip

echo 'export PATH=$PATH:~/bowtie2-2.3.2' >> ~/.bashrc
source ~/.bashrc
sudo apt-get -y install samtools

Alright, now onto a complete re-analysis of our data with the anvi’o pipeline.

Getting it into Anvi’o format¶

Anvi’o takes in 1) an assembly and 2) the raw read data. We have both of those already created, so let’s go ahead and download those data (trimmed reads and asssemblies):

mkdir ~/anvio-work
cd ~/anvio-work

curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gz