1.1 Input Data

Specific information pertaining to the input data for your project is shown in Table 1 below. Check first to ensure that the designations in the Sample Name column match with the associated fastq file(s) in column File Name(s). It is critically important that these designations are correct. If sequencing included a paired-end read, both fastq files should be present in the list of file names used. If data were your data uses a Unique Molecular Index (UMI) to mark duplicate sequencing fragments, then an additional column will describe the source of UMI information.

The deliverable directory contains numerous files consisting of alignment maps, variant calls, and quality control metrics. These files contain the core information from your project. Content description and explanations of these files are given in the appropriate sections below. Many of these files contain metadata which contain additional information which we do not explicitly pull into this report. Information such as the sequencing flowcell, instrument ID, and cluster are embedded in the names of the sequence alignments.

1.2 Fastq Format

A FASTQ file normally uses four lines per sequence.
Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description.
Line 2 is the raw sequence letters.
Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.
Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.

For Illumina FASTQ files the sequence identifier contains information that can be used to identify metadata about the run.
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read>:<is filtered>:<control number>:<sample number>

2.1 Read Trimming

During sequencing library construction, short oligonucleotides are ligated to the ends DNA fragments to be sequenced, so that they can be combined with primers for PCR amplification. While sequencing the fragments, if the read length is greater than that of the target DNA, the adapter sequence next to the unknown DNA sequence of interest is also sequenced, sometimes only partially. To recover the true target DNA sequence, it is required to identify the adapter sequence and remove it. Moreover, to improve further the quality of the sequence reads, skewer also trims the 3’ end of the fragment until a Phred quality of 20 is reached. This is equivalent to a base call accuracy of 99%. More information regarding Phred scores can be obtained from Ewing et al. (1998) and in Ewing and Green (1998).

The Bioinformatics Resource Center performs a quality control step to remove any sequencing adapters and low quality bases. The trimming software skewer (Jiang et al. 2014) was used to pre-process raw fastq files. Skewer implements an efficient dynamic programming algorithm designed to remove adapters and/or primers from sequence reads. Reads that are trimmed too short to be of usable length are discarded from downstream analysis along with the corresponding read pair, in the case of paired-end sequencing.

2.2 Sequence Quality Statistics

After quality filtering, BRC computes QC statistics on your read data including: read length and GC content distribution.

The graphs of summarized data should be fairly uniform for all samples in each QC metric. If one or more metrics in a sample [library] are considerably different, or if they are all highly variable, it may be indicative of a problem that warrants additional attention.

3.1 Read Mapping

The trimmed single- or paired-end reads are aligned against the selected reference genome sequence using BWA-MEM (Li 2013). BWA is a software package for mapping low-divergent sequences against a large reference genome. The BWA-MEM algorithm works by seeding alignments with maximal exact matches (MEMs) and then extending seeds with the affine-gap Smith-Waterman algorithm (SW) through local sequence alignment. BWA may produce multiple primary alignments for different part of a query sequence. This is a crucial feature for long sequences. However, some downstream tools do not work with split alignments.

BWA also finds all hits containing one more mismatches. If the best hit is highly repetitive, then BWA finds all equally best hits only. Base quality is NOT considered in evaluating hits. In the paired-end mode, BWA pairs all hits that it finds. It further performs Smith-Waterman alignment for unmapped reads to rescue reads with a high error rate, and for high-quality anomalous pairs to fix potential alignment errors. The final output from sequence alignment will be a Sequence Alignment Map (SAM)) file. BRC provides the binary (BAM) equivalent of these files, which has a smaller file size but contains the same information as SAM. The binary BAM file can be converted back into a human readable text format by using the SAMTOOLS view command.

3.2 Regions of Interest

The input reads were mapped to the full reference genonme. However, many experiments use capture or amplicon enrichment methods to selectivly target regions of the geneome that are most relavent to the project. For most projects these regions represent a the exons plus some padding bases containing the surrounding introninc sequence. We report metrics for only the regions of interest specifed. The metrics will be highly skewed if we were to report on regions of the genome for which there were no data. It is also important to know that we call variants from only the regions of interst. We provide the a gzipped file with the merged regions of interest for this project.

3.3 Mapping statistics

While informative, the sequence-only metrics reported above are alone insufficient to assess the usability of sequencing data. Data summaries derived from read mapping to the a reference provide additional QC information. For example, sequencing depth should be sufficient to perform basic downstream analyses. As part of this report, mapping statistics for each fastq file were also computed to help gauge various QC metrics including comparisons at the (1) sample-level (2) chromosome-level and (3) aggregate gene-level.

Sample-level information provides the broadest level of mapping QC and is shown in Table 3. It includes the total number of alignments and percentages distributed across primary and secondary alignments. This report is best used to gauge how well reads from each sample mapped (aligned) to the specific reference assembly chosen. Poor primary mapping (unique alignments) may indicate an inappropriate or poorly assembled reference genome. It may also suggest significant repeated regions, which would be reflected in an increase in secondary mapping (non-unique alignments).

The BRC recommends that you check visually the results of the sample-level mappings to ensure the reads align to the expected regions. The Integrative Genomics Viewer (IGV) application is a software that can be downloaded with an academic email address from IGV. It is platform independent and only requires an up-to-date Java installation.

Indexed BAM files of all alignments are available in the delivery directory. To view the coverages, first start IGV and load your genome of interest (IGV has most genomes available for download from within the IGV browser). Second, select Load from File… from the File menu. A dialog box will open. Navigate to and select the BAM file corresponding to the specific sampleName you are interested in viewing.

3.4 Post Alignment Processing

After sequencing reads are aligned to the reference, they are then sorted by position using samtools (Li et al. 2009). Position sorted SAM/BAM files can be indexed for random access so that the entire file does not have to be read into memory when displaying sequence alignments for a given locus. Libraries with high sequencing depth may sample the same fragment of DNA multiple times leading to duplication of reads. Furthermore, the PCR amplification step in the library preparation may preferentially amplify certain fragments due to their length and/or sequence context. BRC marks these duplicate reads using the Picard markduplicates function. Picard uses the position of the fragment ends to determine if alignments are likely to come from the same original fragment of input DNA.

Some projects use of a Unique Molecular Index (UMI). UMIs tag each molecule of input DNA with a unique barcode of random nucleotides so that the UMI barcode can be used to distinguish if alignments with the same positional ends actually come from a two separate input DNA molecules or if the alignments are duplicate reads from a single input fragment. The use of UMIs is important if you wish to have accurate depth of coverage for projects with extremely high sequencing depth or projects that utilize non-random fragment selection (i.e. amplicon based enrichment or enzymatic shearing) during library construction. BRC uses Je, a tool that uses the Picard API, to process the UMI information if applicable (Girardot et al. 2016).

The "On-Target" percentage and "Mean Depth" count only the sequence reads that overlap with the regions of interest specified for this report (see above). Therefore these values are highly dependent on what regions of interest are defined for this project. For example, if the regions of interest represent a subset of genes from an exome test you can expect the On-Target percentage to be a reflection of how much of your total sequenced data corresponds to your subsetted genes. Likewise, the mean coverage represents the mean coverage depth for only the regions defined in the regions of interest.

4.1 Variant Call Format (VCF)

A variant call format (VCF) file is the most common way to communicate a list of variations from the reference genome. VCF is a flexible text based format that is often stored in a gzip compressed format (*.vcf.gz extension). A detailed description of the nuances of this file format can be found at here. Briefly, the format contains meta-information lines, a header line, and then data lines. Each data line is tab-delimited with 8 fixed fields per record, described in the table below. The INFO field is the fixed field that contains a large amount of metadata about the variant. The data contained in the INFO field is variable and will be specified in the VCF header section.

Following the 8 fixed fields, there are optional fields describing the genotype information for each sample (the same types of genotype data must be present for all samples). A FORMAT field is given specifying the order of the colon-separated genotype data for each sample field. See the VCF documentation for more information about genotype fields.

4.2 Small Variant Statistics

The Single nucleotide variant (SNV) constitutes the largest class of variation from a reference genome. The most frequent variations in a population are known as polymorphisms. Thus, this type of variant is sometimes referred to as a single nucleotide polymorphism (SNP). A typical human genome has almost 4 million SNPs. Only a small fraction of these occur the coding region of the genome, commonly known as the exome.

Another common class of genomic variation is the small insertion, deletion, or InDel (insertion and deletion). These variants can have implications for significant changes to the transcripts and/or proteins if they occur in a coding region or splice site. If the length of an InDel is not a multiple of three nucleotides, the reading frame of the the transcript will be altered with a frame shift. The statistics for small InDel mutations are shown in the table below.

4.3 Large Variant Detection

A significant portion of genetic differences between genomes comes in the form of gross changes to the number and/or position of genomic regions. Copy number variations (CNVs) are a class of genomic variant that leads to a gain or loss of a relatively large fragment of the genome. Many CNVs are the result of structural variations. In general, a structural variation can be a deletion, duplication, insertion, inversion, or translocation event that changes the chromosome on a relatively large scale. Large genomic variants can be more difficult to detect with short read technologies. Many methods have been developed to try to address this problem, several of which use a matched control sample (which may not be available for every project).

4.4 Variant Annotation

Following variant detection, BRC annotates the resulting list of mutations with the variant annotation tool SNPeff (Cingolani et al. 2013). SnpEff uses the known gene annotations of the reference genome to predict coding effects such as: synonymous or non-synonymous amino acid changes, gains or losses of start/stop codons, and frame shifts due to InDels. Annotations also include introns, upstream/downstream untranslated regions, splice sites, or intergenic regions. The tool classifies the predicted mutation according to its predicted impact (HIGH, MODERATE or LOW) on gene function. The information from the SNPeff annotation can be found in the INFO field for each variant in the VCF file.

5.1 Did it work?

The identification of mutations in which you are most interested will undoubtedly dictate a measure of “success” or “failure” regarding your experiment. For example, if you had no expectations a priori for any specific mutations and simply wanted a candidate list of variants found in your sample, no matter the outcome, you would likely deem the experiment a success.

Alternatively, if you were interested in testing for specific mutations among a known group of candidate genes, but the analysis showed no evidence of any mutations in your candidate list, you may hastily - and perhaps erroneously - conclude that the experiment was a “failure”.

To evaluate how confident you can be in the results, you must go beyond the variant list itself, and evaluate additional aspects of the analyses and perform some basic checks. Make sure that you review the QC metrics and check that you had sufficient coverage in the regions of interest. There are limitations to calling variants when insufficient data is available for the entire region of interest. It is also possible that the experiment methodology was not appropriate to detect the specific type of mutations needed. For example, detecting structural variants using short read sequencing data is complicated and does not always allow for the detection of complex variations from the reference.

5.2 Where do I go from here?

In nearly all cases, creating a list of variants is not the final step of the analysis. The American College of Medical Genetics and Genomics (ACMG) has published recommendations on how to best interpret the pathogenicity of variants (Richards et al. 2015). At a minimum, further biological insight into an experiment may be gained by examination of genes identified with suspected deleterious mutations, particularly if the prior knowledge of genes or gene pathways involved may highlight specific variants that deserve further scrutiny. You may want to consider the expected inheritance pattern for the phenotype. For a simple recessive phenotype with a Mendelian inheritance pattern, you would expect to find a homozygous mutation or multiple deleterious heterozygous mutations in the responsible gene. Complex traits may be expected to have contributions from several loci. You may consider using Phenotype Ontology terms to narrow your candidate gene search to a list of genes which are known to be associated with the observed phenotype. You may wish to explore additional methodologies like Gene Set Enrichment Analysis (GESA) (Subramanian et al. 2005). Additional insight may also be gained by establishing connections with other experimental data such as RNAseq.

Please feel free to contact our office should you have any questions regarding your data.

6.1 Bioinformatics Resource Center

Bioinformatics Resource Center
Genetics & Biotechnology Center, Room 2130
425 Henry Mall, Madison WI 53706
Email: brc@biotech.wisc.edu
Website: https://bioinformatics.biotech.wisc.edu

6.2 Software

Platform: x86_64-pc-linux-gnu (64-bit)

6.3 References