...
| Tip |
|---|
Unlike other times in the class where we are concerned about being good TACC citizens and not hurting other people by the programs we run, assembly programs are exceptionally memory intensive and attempting to run on the head node may result in the program returning a memory error rather than useable results. When it comes time to assemble your own reference genome, remember to give each sample its own compute node rather than having multiple samples split a single node. |
Assembling even small bacterial genomes can be incredibly time intensive (as well as memory intensive as highlighted above). Fortunately for this class, we can make use of the plasmid spades option to assemble and even smaller plasmid genome that is ~2000 bp long in only a few minutes. I analyzing this data on an idev node and then submitting the other data analysis for the bacterial genomes as a job to run overnight.
Data
Download the paired end fastq files which have had their adapters trimmed from the $BI/gva_course/Assembly/ directory.
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
cp $BI/gva_course/Assembly/*.fastq.gz $SCRATCH/GVA_SPAdes_tutorial
cd $SCRATCH/GVA_SPAdes_tutorial |
SPAdes Assembly
Now let's use SPAdes to assemble the reads. As always its a good idea to get a look at what kind of options the program accepts using the -h option. SPAdes is actually written in python and the base script name is "spades.py". There are additional scripts that change many of the default options such as metaspades.py, plasmidspades.py, and rnaspades.py or these options can be set from the main spades.py script with the flags --meta, --plasmid, --rna respectively. For this tutorial lets use plasmidspades.py
| Expand | ||
|---|---|---|
| ||
The first option in the basic option is: -o<output_dir>directory to store all the resulting files (required) And we will need to supply the read files to the program. In this case we are looking for the following options: -1<filename>file with forward paired-end reads -2<filename>file with reverse paired-end reads |
Once you have figured out what options you need to use see if you can come up with a command to run on the single end and have the output go into a new directory called single_end using all 48 threads that are available (-t 48). The following command is expected to take less than 2 minutes.
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
plasmidspades.py -t 48 -o plasmid -1 SH1_1P.fastq.gz -2 SH1_2P.fastq.gz |
Set 2: SPAdes example Data
This assembly may take more than 2 hours to complete and like the next data set it may be better to run this as a job overnight and analyze the data in tomorrow's class. If you choose to attempt to run this data as an interactive session remember that if your idev session runs out of time, the program will stop working and you will have to start over. If you want to submit it as a job, take the commands you want to run (ie the commands that start with "spades.py" and put them in a file named 'commands' using nano. You can include both commands for this data set and the next data set in the same commands file. Finally use the 'launcher_creator.py' command to generate a ".slurm" file and check with instructor that everything is set up correctly so the job will run overnight.
Data
SPAdes provides two example data sets and benchmarks for how long the pipeline might take to run http://spades.bioinf.spbau.ru/release3.11.1/manual.html#sec1.3. In this section we will download the data for the standard E. coli data and compare their results with ours.
| Code Block | ||||
|---|---|---|---|---|
| ||||
mkdir $SCRATCH/GVA_SPAdes_tutorial # you likely already did this when you ran the selftest cd $SCRATCH/GVA_SPAdes_tutorial wget http://spades.bioinf.spbau.ru/spades_test_datasets/ecoli_mc/s_6_1.fastq.gz wget http://spades.bioinf.spbau.ru/spades_test_datasets/ecoli_mc/s_6_2.fastq.gz |
SPAdes Assembly
Now let's use SPAdes to assemble the reads. As always its a good idea to get a look at what kind of options the program accepts using the -h option. SPAdes is actually written in python and the base script name is "spades.py". There are additional scripts that change many of the default options such as metaspades.py, plasmidspades.py, and rnaspades.py or these options can be set from the main spades.py script with the flags --meta, --plasmid, --rna respectively.
...
Program self tests are typically safe to run on the head node, but the rest of the Tutorial assumes that you are on an idev node. If you are not sure please ask for help.
Data
| Code Block | ||
|---|---|---|
| ||
mkdir $SCRATCH/GVA_SPAdes_tutorial # you likely already did this when you ran the selftest cp $BI/ngs_course/velvet/data/*/* $SCRATCH/GVA_SPAdes_tutorial cd $SCRATCH/GVA_SPAdes_tutorial |
...
Often your read pairs will be "separate" with the corresponding paired reads at the same index in two different files (each with exactly the same number of reads).
SPAdes Assembly
Now let's use SPAdes to assemble the reads. As always its a good idea to get a look at what kind of options the program accepts using the -h option. SPAdes is actually written in python and the base script name is "spades.py". There are additional scripts that change many of the default options such as metaspades.py, plasmidspades.py, and rnaspades.py or these options can be set from the main spades.py script with the flags --meta, --plasmid, --rna respectively.
...
Explore each output directory that was created for each set of reads you interrogated. The actionable information is in the contigs.fa fasta file. The contig file is a fasta file ordered based on the length of the individual contig in decreasing order. The names of each individual contig lists the number of the contig (largest contig being named NODE_1 next largest being named NODE_2 and so on) followed by the length of the contig, and the coverage (labeled as cov on the line). Generally, the lower number of total contigs and the larger the length of each are regarded as better assemblies, but the number of chromosomes present in the organism is an important factor as well.
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
# Count the total number of contigs: grep -c "^>" single_end/contigs.fafasta # Determine the length of the 5 largest contigs: grep "^>" single_end/contigs.fafasta | head -n 5 # Determine the length of the 20 smallest contigs: grep "^>" single_end/contigs.fafasta | tail -n 20 # Determine the length of the 100th through 110th contigs: grep "^>" single_end/contigs.fafasta | head -n 110 | tail -n 10 |
If you ran multiple different combinations of reads for the simulated data how did the insert size effect the number of contigs? the length of the largest contigs? Why might larger insert sizes not help things very much?
| Expand | ||
|---|---|---|
| ||
With better read pairs that link more distant locations in the genome, there are fewer contigs, and contigs are are longer, giving us a more complete picture of linkage across the genomeThe length of repetitive elements in the genome plays a large role in how easily it can be assembled as large repeats need even larger insert sizes to be spanned by single read pairs. |
The complete E. coli genome is about 4.6 Mb. Why weren't we able to assemble it, even with the "perfect" simulated data?
| Expand | ||||
|---|---|---|---|---|
| Sometimes errors in reads lead to dead-ends in the graphs that are trimmed when they should not be.||||
There are 7 nearly identical ribosomal RNA operons in E. coli spaced throughout the chromosome. Since each is >3000 bases, contigs cannot be connected across them using this data. |
What comes next when working with your own data?
...