Hello, my name is Susan Hsiao, and I am an assistant professor and director of bioinformatics in the Laboratory of Personalized Genomic Medicine at Columbia University Medical Center. Today I will be speaking on basic principles in bioinformatics for the identification and interpretation of genomic alterations. I have the following disclosures. I've received honoraria from Illumina, Loxo Oncology, Medscape, and research funding from Bristol-Myers Squibb. So let's start out with this slide. This slide is showing the infrastructure requirements for clinical next-generation sequencing. As you can see, there is data flowing and interfacing between several systems throughout the entire process, talking to systems for test ordering, accessioning, billing, and reporting into the EMR. This talk is really only going to focus on the aspects that relate to generating, interpreting next-generation sequencing data, so really just the middle part of the blue box. But I think it's helpful to kind of keep this in mind as you think about this. It's really bioinformatics analysis is just one small part of the infrastructure that needs to be in place or set up when you're setting up next-generation sequencing testing. There are a large number of tools and databases available for analysis and interpretation. This picture here is from a review in 2016, and the number of tools out there that are constantly being developed, this number just continually increases. My talk today is not going to focus on the specific methodology of any one of these tools, and it's not going to involve programming or biostatistics. The talk today is really a basic talk on principles and approaches to NGS data analysis. And my goal is that by the end of this talk, you'll all feel comfortable with understanding the principles of, the format of, and how NGS data is generated and analyzed. Generally speaking, you can break the bioinformatics process down into primary, secondary, and tertiary analyses. Primary analysis, so in this yellow box here, is the process for sequence generation. This is the process where the signal from the sequencing instrument is converted into reads. From there we go on to secondary analysis, which is shown in the blue box here. In secondary analysis, the sequence is processed, the reads are aligned to a reference genome, and variants are called. Finally, tertiary analysis, which is shown in the red box, this is the process where using various databases and algorithms, variants that we identified are annotated and interpreted. So we're going to go into more detail about each of these processes in the next few slides. Sequence generation, as we heard from Dr. Martis' talk, can be, for example, by the detection of a flow of ions or the detection of a fluorescent signal. This data is converted into nucleotide sequences and stored in FASTQ format. FASTQ format is a text format. If you have access to or you download a FASTQ file, you can actually just open this up in a text editor such as WordPad or Notepad. What you'd see is that it consists of four lines, as is shown in this example here. The first line begins with the at symbol and is followed by a sequence identifier. The second line is the actual nucleotide sequence. The third line begins with the plus symbol and is optionally followed by a sequence identifier. And finally, the fourth line is the corresponding quality value for each nucleotide. This quality value is in ASCII format, so you can look it up, actually, if you're curious what score each symbol corresponds to. And so the second and the fourth lines are going to be the same length because it's the nucleotide sequence and each corresponding quality value. So I mentioned quality value, what is this quality score, what does it mean? Quality scores are typically expressed as a FRED score or Q. This is a property that is logarithmically related to the base calling error probability. So a FRED score of 20, for example, gives a probability of 1 in 100 for an incorrect base call and an accuracy of 99%. A FRED score of 30 gives a probability of 1 in 1000 for an incorrect base call and an accuracy of 99.9%. Typically for most clinical sequencing, we want data with at least a Q score of 20 or more often at least 30. So how do we evaluate the quality of the sequence that was generated? There's several ways. I'm showing here one popular tool that's utilized by many laboratories, which is FastQC. This has a really nice way of visualizing the quality of your sequence data. So an example is shown here. The x-axis is the position in the read and the y-axis is the quality score. The image at the left shows good quality data with high Q scores. You see they're all in this green region. They're all well above a Q score of 30. The image at the right shows poor quality data. As we move along the position of the read, the quality scores go down. There's a wide range of them and you can see how it drops down to single digits even. Here's another visualization of the quality of the data from FastQC. This shows the average quality per read. On the left, we see a nice tight distribution of read quality and the average is about 36 or 37, so great quality score. The right side shows poor quality reads. If you look at the distribution, you see the distribution is wide and the average quality is lower. Sequencers often come in with built-in tools to visualize quality data for the sequence that is generated. This is an example of the sequence analysis viewer for Illumina instruments. This viewer provides information about the cluster density and the flow cell, the Q score bi-cycle, and the Q score distribution. This is all very helpful information when you're trying to troubleshoot or run and look at if there are any issues with the quality of the sequence generation. Here's an example of quality control data from the Torrent Suite software for ion-torrent instruments. This is similar, it shows information about the ion sphere particle density and what percent were these reads were clonal or polyclonal and read length. Okay, so we've covered primary analysis, we've talked about sequence generation, the data format that the reads are in, and how we assess the quality control of the sequence generation process. Let's move on to talk about secondary analysis. So mapping and alignment. Many people use these two terms interchangeably, but they're actually somewhat different. Mapping is a process to identify the origin of the sequence. Alignment is a process to identify differences between sequences. This picture shows how various short reads are mapping and aligning to a reference sequence. Generally speaking, for most applications, you're doing both mapping and alignment. Most aligners will do both. But in theory, you could choose to just perform mapping without the need to identify differences between sequences. Maybe if you were doing transcriptome analysis, if you only want to, you know, identify transcripts and count them, you don't really necessarily need to know the difference between sequences. Or in a metagenomics type of application, perhaps you're only looking for viral or bacterial sequences, and you're just mapping them to a reference, you know, not necessarily looking for differences. So how is alignment data stored? A very commonly used format for storing alignments is the SAM file format. There's also the binary version of this, which is called BAM format. This is an alignment format for storing read alignments against reference sequences. It consists of a header section. This is optional, begins with an at symbol, and also has an alignment section with 11 mandatory fields. I'm going to go into a little more detail about some of these fields and what they mean in the next slide. So this is an example of what a SAM file could look like. So at the top, labeled chord is just the numbering of my reference sequence. The next line, ref, is the reference sequence. And the third line, R001, is the read that we're trying to align. So this is my example, and then next is how we would describe this alignment in a SAM file. So again, the at symbol is the headers, so the first two lines are headers with just some information, and then the third line is our actual alignment section. So the read here, the read name here, like I said, is R001, so that's the first field. This also has information about the reference sequence name, ref, and position of align read is the fourth field. So if you look at our example above, you can see that our read, R001, starts aligning at position 7 of the reference read. Next we have mapping quality, which is expressed in a FRED scale, so FRED score of 30 is pretty good. The next field is the cigar string. At first glance, this looks like some random numbers and letters, but this is actually describing the alignment that we see. Let's break this down. First part of it starts with 8m. If you look in the example above, at R001, the first 8 bases match to the reference sequence. Next part of the cigar string is 2i, and it's saying that there are two insertions, and again, when we look at our example here, we see an insertion of an A and a G relative to the reference sequence. That's followed by 4m, so there's four more bases that are matching to the reference sequence, So how do we assess mapping quality and what are some factors that can affect our ability to map and align? Well, as I mentioned before, quality is expressed on a FRED scale, so that's one thing to look at if you're having any issues with your mapping quality. Other factors to keep in mind are that mapping is affected by repetitive regions and also can be affected by pseudogenes. If the reads themselves are made up of bases with low quality, that will also affect your mapping. And finally, larger indels and structural variants may require specific considerations, specific settings, or specific software to detect. The next part of sequence processing is variant calling. Again, there's a wide variety of variant caller options out there. Most variant callers use an algorithm, either based on a probabilistic model or some heuristics, to predict the likelihood of a variant versus artifact at each locus. And this is based on the quality scores and the allele counts of the reads at that locus. An example here on the right, you can see a homozygous variant in individual 1 and a heterozygous variant in individual 2. We need to have sufficient reads of good quality here for this to be distinguished from an artifact. Additionally, many variant callers include filters such as strand bias filter, filters for homopolymer or low complexity regions. And these also can help reduce some of these artifacts. Variant calls are typically stored in variant call format, or VCF for short. This is also a text-based file. This has several header rows with information about each of these fields. And then in the body, we have the list of variants and we can find information such as a chromosome, physician, any identifiers, reference and alternate allele coverage and quality scores. So a nice compact way of storing our variant calls. Some important considerations to keep in mind for variant calling. Limited nucleic acid input amount. So if you have small, limited samples, these can lead to large numbers of PCR duplicates, especially as seen in amplicon-based assays. Occasionally, we will need higher sensitivity to detect low variant nucleic fractions. This could be for if you're looking for somatic variants or you're looking for variants that might be mosaic. You can include unique molecular identifiers. This utilizes base call consensus and UMI counting to reduce DNA polymerase or sequencing errors and improve accuracy of variants at very low variant allelic fractions. Final consideration is making sure with your variant calling software that you can handle control samples. For example, if you're running a trio for an exome or if you're running paired normal tissue if you're looking at a cancer sample. So how do we look at the quality control metrics for the sequencing process that we just talked about? So let's start with the alignment part. There's a couple of parameters that you can kind of look at and keep in mind. The first one is aligned reads. Next is percentage of on-target reads, and third is coverage. We'll talk in a little bit more detail about each of these concepts in the next few slides. Read duplicates. So what are they exactly? They arise from PCR amplification of adapter ligated DNA fragments. Duplication rates vary based on the amount of starting material used for the library preparation or when there's varied fragment size. That's because smaller fragments can get preferentially amplified. So shown in this image here, you can see some duplicate reads. So these are sort of like the third, fourth, and fifth reads. They all start at the same position. This can be removed during the alignment process. This depends on your aligner and what setting you use. Not every aligner does this. But if you really want to be looking at unique reads, that's a consideration you might want to set up. Next let's talk about on-target reads. So how do you determine if the reads you obtained are in the area you're interested? Well, generally speaking, you want to define your region of interest, or ROI. Frequently this is done using a BED file. So a BED file is a tab-delimited text file. It has three required fields, nine additional optional fields. The three required bed fields are Chrome, or C-H-R-O-M, this is the name of the chromosome, and ChromeStart and ChromeEnd. So that's the starting and ending position of the feature of the chromosome, whatever feature you're looking at, for example, exons. And this will allow you to look at your on-target and off-target reads, which can be very helpful to determine if you were able to successfully enrich for your regions of interest during your library preparation. Now let's talk about depth, or coverage. So this is the number of independent times each base in a targeted region has been sequenced. This is typically expressed as an average X coverage. For example, 20X is an average of 20 sequence reads per base. This can be calculated as shown here, the coverage is the read length times the number of reads over the haploid genome length. The average coverage that you get depends on the assay and to what depth you want to sequence it, and your validation for whatever sensitivity and specificity you wanted. A typical whole genome might have a coverage of, say, maybe 30 to 60X, a whole exome might have a coverage of 80 to 150X, and a targeted amplicon hotspot panel might have a coverage of 1,000 or 5,000X. So that was how we assess some of the alignment part of things. Now let's talk about quality metrics for the barren calls. One metric that a lot of people look at is the transition to transversion ratio, which is simply computed as the number of transition SMPs divided by the number of transversion SMPs. So as you can see in this image here, the number of possible transversions is twice as many as the number of possible transitions. This would lead to a transition-transversion ratio of 0.5 if mutations occurred at equal rates. They don't actually do. Studies have found that actual transition-transversion ratio differs by genomic regions. So it's approximately 2 for whole genome and approximately 3 for whole exome data. Another quality metric that some people use is the heterozygous-homozygous ratio, which is simply the number of heterozygous variants over the number of homozygous variants. It is 2.0 for variants that are in Hardy-Weinberg equilibrium. But again, studies have shown that this is actually population dependent. So it is about 2.0 for folks of African descent, 1.7 for Americans, 1.6 for Europeans, and 1.4 for Asians. So if you want to use this, you know, it might be worthwhile to keep in mind the race of, you know, your patient population and, you know, as this will vary. Other variant-level quality metrics to be aware of include strand bias. This is a disagreement in the genotype inferred from the forward and reverse strands. So sometimes you'll see a variant called only in the forward strand and not the reverse strand. So this is something to be aware of. Usually this will result in a low quality score and your variant caller might already have a filter to remove these calls. Another thing to keep in mind is allele balance. So allele balance is the percentage of reads that support the alternative allele. Generally speaking, this should be about 50% for heterozygous variants. And of course, if you're looking at somatic variants, that's going to vary a little bit. Finally, another metric to keep in mind is SNP density. So this is the frequency of variants observed within a region. Too high of a value without any other explanation suggests noise in the data. So, you know, I mean, maybe you're sequencing a tumor that has MSI. There you might expect a large number of variants, but generally speaking, you know, if you're sequencing, you know, normal or constitutional samples, you really wouldn't expect to see something like, you know, 10 or 20 mutations per megabase. So there you might want to do a closer look at your data and the variant calls. So let's take a look at an example of some quality metrics. So in this first example, this is some metrics from a large hybrid capture-based targeted cancer panel. It's pretty large. It's several megabases in size and it covers selected exons and introns. So total reads 93 million. The proportion of reads that align and are that are on target looks pretty good. Average coverage 950x looks good for a tumor case and the percent of the region of interest that was covered was 98%. So all looks pretty good. And it's all suggesting no real issues with the library or the library preparation. Let's take a look at some of these other metrics listed here, which will give us an idea about the variant calls. Transition-transversion ratio 2.2. And this was about average for what is typically observed in this assay. So not too worried about that. And here, heterozygous-homozygous ratio. Here this is 0.13. This is way off. Very different. I mean, as I mentioned a couple slides ago, this heterozygous-homozygous ratio can vary by population, but this is way different from what we expect. So what happened here? You guys have any ideas or any thoughts about what might have happened? So this case actually turned out not to be a patient sample. This was a proficiency sample and the sample was derived from a cell line. So it had a very large number of homozygous variants, which threw this ratio off. But when we first, we're just looking at this quality metrics, we didn't know any of this. And so this, you know, got flagged for review. Let's look at another case. This is a sample from a constitutional whole exome sequencing case. So again, let's start at the top. Let's look at total reads, 79 million. The line reads and on target reads, all looks pretty good. Average coverage, 112, again, in line with what was expected. Region of interest covered, 98%. Again, that's great. Looking at the variant metrics. So the transition-transversion ratio, 2.94. It's pretty close to the three that we generally expect for whole exome data. So I'm fine with that. The only metric that looks slightly off is the heterozygous-homozygous ratio, which again, varies by population. Generally for the population that we sequenced though, it's a little bit higher than what we normally see. It actually was above the 90th percentile for what we usually observe. So this raised a bit of a red flag. So we took a deeper dive into this case. What we did was take a look at the distribution of variant allelic fractions of the variants in the case. So on the left is a control sample, and we see a nice distribution of variant allelic fractions centered around 45 or 46, which is all our heterozygous variants, and then another peak at 100% for the homozygous variants. Our case, which is seen on the right, you see a wide distribution of variant allelic fractions. Not only that, there's a larger number of heterozygous variants than normal, which is what led to our slightly elevated heterozygous to homozygous ratio. So what do you guys think? What could be causing this? The answer actually is contamination. We worked this case up and actually found that contamination happened during DNA extraction due to a loose O-ring on the instrument. So we had to get more sample and re-extract this case before we could report the results. So this brings us to one more important quality metric to keep in mind, which is contamination. So you can have contamination with microbial sequences, for example, and in that case you would probably see unmapped reads, a lot more reads than that map, and if you took those unmapped reads and tried to blast them, those could map to fungal, bacterial, or viral genome databases. In our example, we knew this wasn't the case because we actually had our reads, they were all aligning, and they were all on target. So they're mapping to the human genome reference, and they were on our targeted region. So we didn't think that was the case. But the other possibility, which is what happened in this case, is cross-sample contamination. So you have a different sample, and how you can detect this, there's a couple ways. You can look at allelic ratios, which is what we did in this case. We just visually examined the allelic ratios and saw that the distribution did not look right. Another possibility for calculating or estimating contamination, and this is something that we did in our first examination, is to estimate it by looking at reference reads at homozygous alternate allele sites. And so that can also give you an idea. So far, we've covered primary and secondary analyses, and the data format, and the quality metrics that you can look at to give you a sense if there were any issues during those steps. The final step is tertiary analysis. This is where variants are annotated and interpreted. There's a variety of databases and resources that provide information on population or disease prevalence, on functional annotation of these variants, and or in silico predictions for variants. So this is just an example of some of the many that are out there. Annotation is really key to interpretation. This provides information on the genes and the regions, presence in population or disease databases, predictive in silico scores. All this information is very helpful in filtering the number of variants to review, and also reducing the time needed to review a case and evaluate if this variant should be reported or not, and if reportable, how the variant should be classified. So you can take your VCF file and you can manually, I mean you can set up a pipeline to annotate it with any one of your favorite databases. Listed here are also some very popular tools that actually make that process a lot easier and integrate a lot of the most popular databases to annotate your VCF. There are guidelines for interpretation and reporting of germline and somatic variants. These have been published and you'll hear a little bit more about these guidelines in the case study section, which will be later in the session. Also I just wanted to show this schematic of clinical guidelines that cover the entire next generation sequencing process from library preparation to reporting. In particular the CAP checklist and the bioinformatics guideline by Roy et al. are pretty good references for additional information on the topics covered by this talk. Okay, so the last part of this talk I wanted to very briefly touch on some additional bioinformatics analyses that we routinely do with NGS data, and that's going to be copy number variant detection and RNA sequencing. So copy number variants include deletions, insertions, and duplications of genomic regions, and are important in a number of germline and somatic processes. Approaches for copy number variant detection include read depth, paired end mapping, and split reads. And actually most of the algorithms use one or multiple of these approaches to detect copy number variants. So read depth analysis, basically this uses read depths to infer copy number. So increased copies correlates with increased numbers of reads. The advantage to using this is it can actually provide an exact copy number. This technique does require the use of a paired case or control sample or a pool of samples for normalization. Paired end mapping identifies copy number alterations and can also identify structural variants. And so how this works is by assessing discordantly mapped paired reads. As you see in this example here, our read maps, our paired reads, sorry, maps at a distance significantly different from the average insert size. And this technique can be a little bit limited by low complexity regions, and if the insertions are larger than the average insert size of the library, that can be a little bit difficult to detect. Lastly we'll talk about the split read approach. How this works is it utilizes read pair where we have one read that aligns to the reference genome, so shown here in the purple read aligns perfectly to the reference genome, but the other pair in red is unmapped or partially mapped. What happens then is that the unmapped or partially mapped read is split and then independently mapped, and this allows you to determine the break points of either a copy number alteration or a structural variant. One limitation to this approach is that the detection is dependent on the length of the reads. Okay, the last topic I'm going to touch on in this talk is a very brief overview of RNA sequencing. So RNA sequencing has a lot of great applications out there. It's a great tool that could allow you to examine the abundance and makeup of all the expressed genes in a sample. This data allows you to compare expression among groups of samples. You could use this to identify fusion transcripts. You can detect transcriptional alterations that are caused by variants of uncertain significance. And finally, you can compare DNA and RNA variants, and this is helpful to verify sample identification and to check for sample contamination. There's other applications out there. I'm just listing some of the more commonly used ones. So shown here is general bioinformatics process for RNA sequencing data, and there's a lot of similarities to the process for DNA data that we already talked about today. We start with the raw sequencing reads, which again have their associated quality scores. There could be preprocessing steps for trimming, quality filtering, removing ribosomal RNA for example. These reads are then mapped, and there's actually a few options here for mapping. You can map to the genome, and if you do this, you want to use a gapped aligner, so a aligner that is aware that if your read is crossing exons, it'll be aware and be able to handle the intronic sequence in the genome. The advantage of doing it this way is that you could detect novel transcripts this way. You could, instead of mapping to the genome, you can map to the transcriptome, and you don't detect novel transcripts this way, but you know that it also works pretty well if you're just trying to count reads. Or alternatively, you could do D-level assembly, just where all the reads are assembled. So following that step, then you can move on to variant analysis. So you could employ various variant calling software out there to detect variants from your RNA-seq data. You can also look at fusions or different isoforms, and also you can take your aligned reads and look at run pipelines to annotate your transcripts, or to quantify your transcripts and do differential expression or abundance estimation. So similar ideas, similar sort of workflow as the DNA sequencing. Okay, so that is all I kind of wanted to cover today in this talk. I hope you found this helpful. Thank you for your attention, and feel free to reach out if you have any questions.

bioinformatics

genomic alterations

clinical next-generation sequencing

FASTQ format

SAM file format

variant calling

copy number variant detection

RNA sequencing