--- id: dunovo name: Calling very rare variants description: >- In this tutorial we are going to perform discovery of low frequency variants from duplex sequencing data. As an example we use the ABL1 dataset published by Schmitt and colleagues (SRA accession SRR1799908). title_default: dunovo steps: - title: Finding rare variants content: >- Most popular variant callers focus on the common case of sequencing a diploid individual to find heterozygous and homozygous variants. This is a well-studied problem with its own challenges, but at least you can expect your variants to be present in either 100%, 50%, or 0% of your sample DNA. If you observe a variant present in 99%, 56%, or 2% of the reads at a site, you can probably assume the allele is actually present at 100%, 50%, or 0%, respectively, in your sample.

But in this tutorial, we’re looking for rare variants. So our true frequency might actually be 13%, 1%, or even 0.4%. The challenge then becomes distinguishing these situations from sequencing errors. Next-generation sequencers produce noise at this level, making it challenging to make this distinction in data produced with standard resequencing methods. backdrop: true - title: Duplex sequencing content: >- Duplex sequencing is a method that addresses the problem of distinguishing sequencing signal from noise. It can increase sequencing accuracy by over four orders of magnitude. Duplex sequencing uses randomly generated oligomers to uniquely tag each fragment in a sample after random shearing. The tagged fragments are then PCR amplified prior to sequencing, so that many reads can be obtained from each original molecule. The tags in each read can then be used to identify which original fragment the read came from. Identifying multiple reads from each fragment allows building a consensus of the original sequence of the fragment, eliminating errors.

The key to duplex sequencing, as opposed to other types of consensus-based methods (review here), is that both ends of the original fragment are tagged such that its strands can be distinguished. Knowing which strand each read comes from allows us to recognize errors even in the first round of PCR. backdrop: true - title: Duplex sequencing content: >- Processing the raw reads into consensus sequences consists of four main steps: Du Novo is a tool which can carry out these steps. Unlike most other such tools, it can do so without the use of a reference sequence, and it can correct for errors in the tags which can contribute to data loss. backdrop: true - title: The value of single-strand consensus sequences content: >- The DCSs have the ultimate accuracy, yet the SSCSs can also be very useful when ampliconic DNA is used as an input to a duplex experiment. Let us illustrate the utility of SSCSs with the following example. Suppose one is interested in quantifying variants in a virus that has a very low titer in body fluids. Since the duplex procedure requires a substantial amount of starting DNA (between 0.2 and 3 micrograms) the virus needs to be enriched. This can be done, for example, with a PCR designed to amplify the entire genome of the virus. Yet the problem is that during the amplification heterologous strands will almost certainly realign to some extent forming heteroduplex molecule. backdrop: true - title: Generating consensus sequences content: >- The starting point of the analysis is sequencing reads (in FASTQ format) produced from a duplex sequencing library. backdrop: true - title: Getting data in and assessing quality content: >- We uploaded the Schmitt et al. 2015 data directly from SRA as shown in this screencast. backdrop: true - title: Importing the raw data element: '#history-options-button' content: >- In order to import the data from another history click on the gear icon at the top. Click on “Import from File” at the bottom of the menu and insert this link in the box under Archived History URL: https://usegalaxy.org/history/export_archive?id=7ac09d1db287dbba placement: left - title: Getting data in and assessing quality element: '#history-options-button' content: >- This created two datasets in our galaxy history: one for forward reads and one for reverse. placement: left - title: Getting data in and assessing quality content: >- We then need to evaluate the quality of the data by running FastQC on both datasets (forward and reverse). You can read about using tool FastQC here. backdrop: true - title: Assessing quality element: '#tool-search-query' content: Search for 'FastQC' tool. placement: right textinsert: FastQC - title: Assessing quality element: '#tool-search' content: Click on the 'FastQC' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fdevteam%2Ffastqc%2Ffastqc%2F0.72"] .tool-old-link - title: Assessing quality element: '#tool-search' content: Run the tool with the default parameters. position: right - title: Assessing quality element: '.history-right-panel .list-items > *:first' content: >- One can see that these data are of excellent quality and no additional processing is required before we can start the actual analysis. position: left - title: Processing reads into duplex consensus sequences with Du Novo content: >- Now we are ready to collapse the raw reads into duplex consensus sequences. backdrop: true - title: Sorting reads into families content: >- The Du Novo: Make families tool will separate the 12bp tags from each read pair and concatenate them into a 24bp barcode. Then it will use the barcodes to sort the reads into families that all descend from the same original fragment. backdrop: true - title: Sorting reads into families element: '#tool-search-query' content: Search for 'Du Novo Make families' tool. placement: right textinsert: Du Novo Make families - title: Sorting reads into families element: '#tool-search' content: Click on the 'Du Novo Make families' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fdunovo%2Fmake_families%2F2.15"] .tool-old-link - title: Sorting reads into families element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Correcting barcodes content: >- The grouping reads based on barcode relies on exact barcode matches. Any PCR or sequencing error in the barcode sequence will prevent the affected reads from being joined with their other family members.
Du Novo includes a tool which can correct most of these errors and recover the affected reads. This can increase the final yield of duplex consensus reads by up to 11%fix link backdrop: true - title: Correcting barcodes element: '#tool-search-query' content: Search for 'Du Novo Correct barcodes' tool. placement: right textinsert: Du Novo Correct barcodes - title: Correcting barcodes element: '#tool-search' content: Click on the 'Du Novo Correct barcodes' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fdunovo%2Fcorrect_barcodes%2F2.15"] .tool-old-link - title: Correcting barcodes element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Aligning families content: >- After grouping reads that came from the same original fragment, we need to align them with each other. This next tool will perform a multiple sequence alignment on each family. backdrop: true - title: Aligning families element: '#tool-search-query' content: Search for 'Du Novo Align families' tool. placement: right textinsert: Du Novo Align families - title: Aligning families element: '#tool-search' content: Click on the 'Du Novo Align families' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fdunovo%2Falign_families%2F2.15"] .tool-old-link - title: Aligning families element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Making consensus sequences content: >- Now, we need to collapse the aligned reads into consensus sequences. This next tool will process each group of aligned reads that came from the same single-stranded family into a consensus. Then it will align the consensus sequences from the two strands of each original molecule, and call a consensus between them.

Normally, the tool only produces the final double-stranded consensus sequences. But we will make use of the single-stranded consensus sequences later, so we’ll tell it to keep those as well. backdrop: true - title: Making consensus sequences element: '#tool-search-query' content: Search for 'Du Novo Make consensus reads' tool. placement: right textinsert: Du Novo Make consensus reads - title: Making consensus sequences element: '#tool-search' content: Click on the 'Du Novo Make consensus reads' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fdunovo%2Fdunovo%2F2.15"] .tool-old-link - title: Making consensus sequences element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Filtering consensuses content: >- You may have realized that when calling a “consensus” between two sequences, if the two disagree on a base, there’s no way to know which is correct. So in these situations, Du Novo uses the IUPAC ambiguity letter for the two different bases (e.g. W = A or T). Also, when calling single-stranded consensus sequences, if there aren’t enough high-quality bases to call a position (during the previous steps we set this threshold to 70%), it gives an N
This information could be useful for some analyses, but not for our variant calling. The tool Sequence Content Trimmer will help with filtering these out. With the settings below, it will move along the read, tracking the frequency of ambiguous (non-ACGT) bases in a 10bp window. If it sees more than 2 ambiguous bases in a window, it will remove the rest of the read, starting with the first offending base in the window. We’ll also tell it to remove entirely any read pair containing a read that got trimmed to less than 50bp. backdrop: true - title: Filtering the consensus sequences element: '#tool-search-query' content: Search for 'Sequence Content Trimmer' tool. placement: right textinsert: Sequence Content Trimmer - title: Filtering the consensus sequences element: '#tool-search' content: Click on the 'Sequence Content Trimmer' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fsequence_content_trimmer%2Fsequence_content_trimmer%2F0.1"] .tool-old-link - title: Filtering the consensus sequences element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Calling variants with duplex consensus sequences content: >- At this point we have trimmed DCSs. We can now proceed to call variants. This involves aligning the variants against the reference genome, then counting variants.
We’re not specifically interested in the reference sequence, since all we care about is sequence content of the consensus reads. But we’ll be using the reference sequence to figure out where all the reads come from. This lets us stack them on top of each other, with equivalent bases lined up in columns. Then we can step through each column, count how many times we see each base, and and compile a list of variants. backdrop: true - title: Align against the genome with BWA-MEM element: '#tool-search-query' content: Search for 'Map with BWA-MEM' tool. placement: right textinsert: Map with BWA-MEM - title: Align against the genome with BWA-MEM element: '#tool-search' content: Click on the 'Map with BWA-MEM' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fdevteam%2Fbwa%2Fbwa_mem%2F0.7.17.1"] .tool-old-link - title: Align against the genome with BWA-MEM element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Left Aligning indels content: >- To normalize the positional distribution of indels we use the tool BamLeftAlign utility from the FreeBayes package.
This is necessary to avoid erroneous polymorphisms flanking regions with indels (e.g., in low complexity loci) backdrop: true - title: Left-align indels element: '#tool-search-query' content: Search for 'BamLeftAlign' tool. placement: right textinsert: BamLeftAlign - title: Left-align indels element: '#tool-search' content: Click on the 'BamLeftAlign' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fdevteam%2Ffreebayes%2Fbamleftalign%2F1.1.0.46-0"] .tool-old-link - title: Left-align indels element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Calling the variants content: >- Now we’ll use our aligned consensus reads to find variants.
Normally, in a diploid resequencing experiment, you would call variants relative to the reference. So, you’d report sites which are different from the reference (and whether they’re hetero- or homozygous). In our case, we’re interested in rare variants. So what we’ll report is the sites where there is more than one allele, and what the frequency is of the less-common allele (the minor allele). This has the potential to include every small sequencing error (even though we’re using duplex, there still are errors). So to reduce the noise, we’ll set a lower threshold at 1% minor allele frequency (MAF). backdrop: true - title: Finding variants in the alignment content: >- To identify sites containing variants we use the tool Naive Variant Caller (NVC) tool. This reads the alignment and counts the number of bases of each type at each site. backdrop: true - title: Count the variants element: '#tool-search-query' content: Search for 'Naive Variant Caller (NVC)' tool. placement: right textinsert: Naive Variant Caller (NVC) - title: Count the variants element: '#tool-search' content: Click on the 'Naive Variant Caller (NVC)' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fblankenberg%2Fnaive_variant_caller%2Fnaive_variant_caller%2F0.0.4"] .tool-old-link - title: Count the variants element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Count the variants element: '#history-options-button' content: >- The tool Naive Variant Caller (NVC) generates a VCF file that can be viewed at genome browsers such as IGV. Yet one rarely finds variants by looking at genome browsers. We’ll want to use tools to search for variants that fit our criteria. placement: left - title: Finding minor alleles content: >- Now we’ll want to parse the VCF produced by the NVC, determine what the major and minor allele is at each site, and calculate their frequencies. The Variant Annotator from the NGS: Variant Analysis section can do this. backdrop: true - title: Read the variants file element: '#tool-search-query' content: Search for 'Variant Annotator' tool. placement: right textinsert: Variant Annotator - title: Read the variants file element: '#tool-search' content: Click on the 'Variant Annotator' tool to open it. placement: right postclick: - >- a[href$="/tool_runner?tool_id=toolshed.g2.bx.psu.edu%2Frepos%2Fnick%2Fallele_counts%2Fallele_counts_1%2F1.2"] .tool-old-link - title: Read the variants file element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Filtering out the noise content: >- Now we have a file containing the base counts for every site covered by at least 10 reads. We’d like to filter through this data to find sites with a reasonable chance of being a real variant, not sequencing error.
The tool Variant Annotator produces a simple tab-delimited file, with one site per line. We can use the tool Filter tool to process this kind of file. We’ll use the filter c16 >= 0.01 to remove lines where the value in column 16 is less than 0.01. Column 16 contains the minor allele frequency, so this will remove all sites with a MAF less than 1%. backdrop: true - title: Filter the raw variants list element: '#tool-search-query' content: Search for 'Filter' tool. placement: right textinsert: Filter - title: Filter the raw variants list element: '#tool-search' content: Click on the 'Filter' tool to open it. placement: right postclick: - 'a[href$="/tool_runner?tool_id=Filter1"] .tool-old-link' - title: Filter the raw variants list element: '#tool-search' content: |- Run the tool with the following parameters: position: right - title: Filter the raw variants list element: '#history-options-button' content: >- The polymorphism we are interested in (and the one reported by Schmitt et al. 2015) is at the position 130,872,141 and has a frequency of 1.3%. The other site (position 130,880,141) is a known common variant rs2227985, which is heterozygous in this sample. placement: left - title: Conclusion content: >- You should now understand duplex sequencing, rare variants, and be able to process the former to find the latter. backdrop: true - title: Key points content: >- backdrop: true