Team III Genome Assembly Group

Group Members: Deepali Kundnani, Aparna Maddala, Swetha Singu, Yiqiong Xiao, Ruize Yang

Lectures

• Team 3 Genome Assembly: Background and Strategy
• Team 3 Genome Assembly: Results

Introduction

Genome assembly is the process of combining the reads generated by sequencing technology into contiguous sequences that ultimately span a large portion of a genome. This can either be accomplished by examining how the reads overlap or examining their similarity to a reference genome. A standard genome assembly pipeline includes quality control, trimming, post-trimming quality control, sequence assembly, and assembly validation [1]. The objective of this project was to assemble bacterial genomes from 50 fastq files provided by CDC PulseNET. All of the samples were generated by paired-end sequencing on the Illumina platform.

Methods

Pipeline for Testing

Figure 1: Proposed assembly pipeline. fastp and MultiQC are used for quality control/trimming, sequence assembling is done in SPAdes, and assembly validation is conducted in Quast and BUSCO. The yellow bubbles represent the output files generated at each step.

Quality Control/Trimming

Tools used:

• fastp (https://github.com/OpenGene/fastp)
• MultiQC (https://multiqc.info/)

Tool Selection

For quality control, we compared two tools, fastp and FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We first proved that the two programs generated identical information when run on identical fastq files, after which we compared the information displayed in the reports for both. While FastQC creates highly informative per-base sequence quality graphs, it runs significantly slower than fastp. We were able to extract biologically significant information from the charts generated from fastp, so we decided to use fastp for quality control.

Afterwards, we compared fastp's trimming features with those of Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Our data did not contain adapters, so we did not need an additional tool like CutAdapt (https://cutadapt.readthedocs.io/en/stable/) to remove them. We showed that most, if not all, of Trimmomatic's trimming features can be replicated in fastp. Furthermore, fastp contains a feature specifically for paired-end data where it can use a high-read to correct low confidence bases in its mate. Fastp has the added advantage of combining both quality control and trimming into a single step, increasing the speed and usability of our pipeline. Since we had 50 input files in our pipeline, we used MultiQC to consolidate the 50 separate quality control reports generated by fastp into a single report.

Pre-Trim Quality Control

Figure 2: Pre-trim quality control curves of our 50 samples. Plot generated using fastp and MultiQC. A) Read 1 quality curves; B) Read 2 quality curves.

Read 1 had high quality reads, so each of the 50 read 1 samples passed the mean quality score threshold of 28. However, read 2 was of lower quality, with only 25 samples passing the same mean quality threshold. This was mostly due to low quality regions at the 5' end of read 2. Because of this notable difference in quality between the two reads, we decided to trim read 1 and read 2 using different parameters. According to a comparison of trimming parameters on the assembly quality of Cyanoderma ruficeps genomes, trimming does not have a clear effect on assembly quality but speeds up most assembly implementations [2]. Therefore, in order to obtain high quality data to pass to the assemblers, we decided to completely trim the low-quality regions at the 3' and 5' ends of each read while retaining the high-quality middle regions.

Trimming Parameters

The following arguments were supplied to fastp in order to trim our data: -f 5 -F 30 -t 10 -e 28 -c -5 3 -M 27.

• -f 5 - globally trims 5 bases from 5' end of mate 1
• -F 30 - globally trims 30 bases from 5' end of mate 2
• -t 10 - globally trims 10 bases from 3' end of both mates
• -e 28 - discards reads with an average quality score under 28
• -c - turns on paired-end base correction, which slightly increased the quality of the 3' end of mate 2
• -5 3 - turns on sliding window trimming from 5' end with a window size of 3
• -M 27 - sets a quality threshold of 27 for sliding window

Assembly

Parameters considered when selecting assembly tool

[namedrop the two comparative studies we referenced here]

Tools considered

MaSuRCA

Version: 3.3.5

MaSuRCA, or the Maryland Super-Read Celera Assembler, uses a combination of the OLC and de Bruijn Graph algorithms to assemble data. MaSuRCA normally creates longer contigs and offers higher N50 values. It works best on untrimmed data, which is possibly responsible for its slower computational time [2].

Results

Post-Trimming Quality Control

Figure 3: Post-trim quality control curves of our 50 samples. Plot generated using fastp and MultiQC. A) Read 1 quality curves; B) Read 2 quality curves.

After trimming, the quality of our read 1 data improved slightly. However, the largest difference can be observed in read 2 data, with all 50 of the samples passing the mean quality threshold of a Phred score of 28. The 5' and 3' ends of some of the read 2 samples are of low quality, but more aggressive trimming would most probably also eliminate high-confidence regions.

Figure 4: Percentage of reads in each sample that passed the filters, generated using fastp and MultiQC.

A vast majority of reads in each sample passed our filter, with at most 11% of reads falling below our quality and length thresholds and at most []% of bases trimmed from each sample.

Inter-Assembler Assessments

For quality assessment of Genome Assembly, we make use of the widely used software QUAST(QUaslity Assessment Tool for Genome Assembly) with BUSCO(Benchmarking Universal Single-Copy Orthologs) to complement the analysis. With Quast, we mainly on N50, L50, number of Contigs, N's per 100kb. BUSCO provides us with single copy, duplicated, missing and fragmented Benchmarking orthologs. Good assemblies are more complete with higher N50 and lower L50. Assemblies with poor error correction lead to higher duplicated BUSCOs and ones with data loss lead to missing or fragmented BUSCOs. Lower completeness scores with many fragmented BUSCOs can also be suggestive of misassemblies that Quast might not be able to detect without alignment to a reference genome and that is where BUSCO provides a better understanding fo the assembly quality.

Figure x: Comparison of 1.SPAdes 2.SKESA 3. Abyss with KmerGenie 4. Masurca on a subset of 10 samples selected based on trimming quality.

We rule out Abyss based on presence of N's per 100Kbp. In the previous analysis of different kmers on Abyss, the genome assemblies showed N's if the the kmers were too low or too high. Its possible that the kmers suggested by KmerGenie are lower than optimal.

It is also evident that MaSuRCA has a lower Completeness score and missing BUSCOs suggesting the loss of valuable information from the Assembly. This again can be due high kmer selection which was seen to be around 99 for most of the samples. MaSuRCA can be a good tool for reads with higher sequencing depth or coverage and for organisms cultured prior to sequencing. It is safe to assume the worst-case scenario of working with non-cultured organisms having lower sequencing depth, we rule out MaSuRCA from the list of Assemblers for final Validation. With MaSuRCA we also run the risk of losing smaller genomes of plasmids that can play a role in horizontal gene transfer.

Since the comparisons were non-conclusive between Spades and SKESA by just using 10 libraries, we assemble genomes from all 50 libraries provided to us for finding the best assembly given this set of samples.

Figure x: Comparison of SPAdes vs SKESA on all given 50 samples

As seen from metrics from QUAST and BUSCO, Spades clearly outperforms SKESA in terms of giving us better N50 and lower L50 values for genomes assembled from the same libraries. With SKESA we also saw a varied distribution of these scores which indicates that SKESA, even though it has a sophisticated way of handling assemblies for each sample, might not be that effective. And when we look at the BUSCO completeness scores with fragmented/missing BUSCOs we see a similar trend.

Figure x: Correlation of Quast and BUSCO scores for each samples

Unlike MaSuRCA, while looking at each sample's N50 with BUSCO score, we do see a correlation of SKESA's lower N50's with lower completeness scores. These comparisons help us conclude that SPAdes is the best tool to be used for Genome Assembly regardless of the quality of the libraries it is provided and hence is the best assembler for our pipeline

References

1. Dominguez Del Angel V, Hjerde E, Sterck L et al. Ten steps to get started in Genome Assembly and Annotation [version 1; peer review: 2 approved]. F1000Research 2018, 7(ELIXIR):148 (https://doi.org/10.12688/f1000research.13598.1)
2. Yang, S.-F.; Lu, C.-W.; Yao, C.-T.; Hung, C.-M. To Trim or Not to Trim: Effects of Read Trimming on the De Novo Genome Assembly of a Widespread East Asian Passerine, the Rufous-Capped Babbler (Cyanoderma ruficeps Blyth). Genes 2019, 10, 737.