Introduction
We analyzed a single lambda phage evolution experiment as both an introduction to variant visualization in the first part of Wednesday's class, but breseq is designed to be able to with "small" genomes, not just "tiny" genomes. It also has more advanced tools for visualizing mutations across multiple genomes and more complete identification of all mutations.
Objectives
- Run breseq on 7 different E. coli clones evolved for 40,000 generations
- Use the packaged gdtools commands to generate a comparison table of all mutations from all 7 samples
- Learn shortcuts for compressing, transfering, and extracting multiple folders/files with single commands
- Visualize the data off of TACC
Table of Contents
Setting up the run
The data files for this example are in the following path: $BI/ngs_course/ecoli_clones/data/ go ahead and copy them to a new folder in your $SCRATCH directory called BDIB_breseq_multi-sample:
if everything worked correctly, you should see the following files. We've provided a bit more context to what those files actually are:
File Name | Description | Sample |
|---|---|---|
| Paired-end Illumina 36-bp reads | 0K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 2K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 5K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 10K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 15K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 20K generation evolved E. coli strain |
| Paired-end Illumina 36-bp reads | 40K generation evolved E. coli strain |
| Reference Genome | E. coli B str. REL606 |
breseq will take a little longer to run on these sequences, so this give us an opportunity to run several commands at the same time making use of the multiple cores on a single processor. You'll want each command (line) in the commands file to look something like this:
breseq -j 6 -r NC_012967.1.gbk -o output_<XX>K SRR030252_1.fastq SRR030252_2.fastq &> <XX>K.log.txt
Notice we've added some additional options:
| part | puprose |
|---|---|
| &> <XX>00K.log.txt | Redirect both the standard output and the standard error streams to a file called <XX>00k.log.txt. It is important that you replace the <XX> to send it to different files, but KEEP the &> as those are telling the command line to send the streams to that file. |
| -o output_<xx>00k | all of those output directories should be put in the specified directory, instead of the current directory. If we don't include this (and change the <XX>), then we will end up writing the output from all of the runs on top of one other. The program will undoubtedly get confused, possibly crash, and generally it will be a mess. |
| & | run the preceding command in the background. This is required so all the commands will run at once |
Executing the run:
chmod +x commands # makes the commands file executable module unload samtools ./commands # executes the commands file. the './' portion tells the command line to execute a file from the current directory as it is not in your $PATH variable.
Since we started running all of the breseq runs in the background, and are redirecting all of the output to log files, we need to make sure that things are actually running as we intended. The following 2 commands are good ways to check on this:
Even with just using the R1 data, this will take several minutes to complete (~40 minutes). While it is running, you should move onto reading through the next section learning how you will create comparison tables to easily visualize the data across the different samples. If you have all that figured out you may want to start ne of the other tutorials while the above commands finish. Just remember to check back periodically as the 'completion' notifications of background jobs can be easy to miss.
Using gdtools to compare samples.
Assuming everything is running as it should (you should ask if you are not sure), we need to start exploring the gdtools commands. gdtools is a series of commands designed to function on "genome diff" (.gd) files (which are the main command line type of breseq command line analysis) that come packaged with breseq. As we have seen with numerous other programs throughout the course typing gdtools by itself will return a list of options. In the case of gdtools, it returns the list of subcommands that are available. Each has its own use, but for the purpose of this tutorial we are interested in comparing the output of multiple samples to one another.
In our command we sent the output of each sample to its own folder using the -o output_<XX>K portion of the command line. Among the multiple directories that breseq created within each output directory, you will find one named output which contains a file named output.gd. These are the files that we need to feed into the gdtools compare command (along with the NC_012967.1.gbk reference file) to generate the comparison table we want. In this course we have mentioned (though not emphasized) that using consitant nomenclature is important which it is, but now we come to a point where the consistancy has created a problem. Because all of the files are named "output.gd" we can neither move all of the commands into a single folder (they would overwrite one another as they have the same name) or the resulting table would display sample IDs as "output" "output" "output" ... not the most helpful information to say the least.
To circumvent that we will need to rename each file to something more meaningful (say the name of the directory that the output/output.gd file was found in (aka the sample name)). While this command can be run on TACC we'll wait until we transfer the files back to our local machine and do it there.
Visualizing output
Can you figure out how to archive all of the output directories and copy only those files (and not all of the very large intermediate files) back to your machine? - without deleting any files?
For 7 samples this is certainly doable, what happens when you have 1,000 samples? If only this was a big enough problem for enough researchers that someone stepped up and wrote a small script to do this in an automated way... Oh wait there is the BioITeam, and within the bin directory you will find a script named export-breseq.sh ... Normally this script assumes that all samples output were directed to a single output folder before being named based on the specific samples. In other words, rather than using -o output_<XX>K as the output directory, we would use -o breseq_output/output_<XX>K so all of the output folders are in a single subdirectory. The export-breseq.sh command will still work, but the output will end up 1 directory higher than where you started ... in this case, it will end up in the $SCRATCH directory.
cds cd BDIB_breseq_multi-sample export-breseq.sh mv ../05_Output* .
Now lets transfer the 05_Output_Export.tar.gz file back to our local computer
Now you will see 7 .tar.gz files in the directory, and you can extract each of them 1 at a time by using the tar -xvzf <FILE> command. Again this seems acceptable for 7 samples but not when dealing with 100s of samples. The 1 liner below will extract .tar.gz files sequentially:
for f in *.tar.gz;do tar xvzf $f;done
This command will extract all files ending in .tar.gz regardless of their origin. As such the following 2 commands might be of use to you:
echo "for f in *.tar.gz;do tar xvzf $f;done" > detar chmod +x detar
Now you have a file named detar that executes the command without having to copy paste it. The only thing left is to move it to somewhere in your path. This will work on either your personal computer or TACC the same way.
Now you can click through each individual sample's output to see the different mutations in any given sample, but the whole point of this tutorial has been to look at multiple samples at once.
Generating compare tables.
As mentioned above, we need to change the names of the output.gd files to something more informative. You can manually do this by going into each sample directory, each output directory, and using the cp command to rename the output.gd file to something more informative.
cd output_00K_output cp output.gd 00K_output.gd# ETC
Once again there is an easier way to do this all at once while generating a new directory to hold all of the outputs:
# navigate to the directory containing all of the extracted output directories: the ls command should return the following folders: # output_00K_output output_05K_output output_15K_output output_40K_output output_02K_output output_10K_output output_20K_output mkdir gds for f in *_output; do n=$(echo $f|sed 's/_output//'); cp $f/output.gd gds/$n.gd;done cd gds
Now that we have nicely named .gd files (the ls command can confirm this), its time to use the gdtools compare command. As your local computer doesn't have breseq on it we will need to move these .gd files back to TACC so we can access the gdtools commands. This is done on purpose to show you how to transfer files from your local computer to TACC (something you will do almost as frequently as the inverse).
Now working again on TACC. Earlier we determined the formatting for the gdtools compare command. Either of the following commands will work. Compare the 2 and decide which you think is better:
gdtools compare -r ../NC_012967.1.gbk output_00K.gd output_02K.gd output_05K.gd output_10K.gd output_15K.gd output_20K.gd output_40K.gd # OR more simply... gdtools compare -r ../NC_012967.1.gbk *.gd
Now transfer the new annotated.html file back to your local computer.
Open the annotated.html and examine the results, does this make sense with what you would expected from an evolution experiment?
