WGS Variant Calling: Variant calling with GATK - Part 1 | Detailed NGS Analysis Workflow

Foreign. Hi guys, welcome back to my channel. I'm so glad to see you all here. So, so far in all my previous videos we have mainly been focusing on RNA SEQ or single cell RNA seq. But in this video we are switching gears and talking about other type of high throughput data I.e. whole genome sequence sequencing. So this video we will be talking about variant calling, which is a very commonly performed analysis on whole genome sequencing or whole exome sequencing data. So today I want to demonstrate how to perform a variant calling that is called variance on whole genome sequencing data using GATK best practice workflow. So the aim of variant calling analysis is to to start with sequencing reads and perform some steps to determine what variants are present in the these reads. So basically the intuition behind this analysis is to generate a file, a variant called format file that is a VCF file that stores variants that are present in the data. So the diagram at the bottom is like a high level overview of the steps that are typically performed in a variant calling analysis. So the DNA sample is first fed into a sequencer which gives out reads which are stored in FAST auto fastql view format and these reads are aligned to a reference genome. So we get aligned reads which are in BAM or SAM file format. And using these aligned reads we employ a variant calling algorithm that essentially calls variants, identifies these variants and store them in a file called variant called format file that is a VCF file. So today we are going to start with sequencing reads and we are going to perform various steps, go into the details of each of these steps and ultimately call variance in our data and generate a VCF file. So since today we are using a GATK best practice workflow, I want to talk about what is gattk? So GATK stands for Genome Analysis Toolkit and is the industry standard toolkit for analyzing and identifying variants. And GATK is basically designed for whole human, whole genome and whole exome sequencing data. So gatki essentially is a package of command line tools that are written in Java. And these tools, they have plethora of tools and these tools can be used individually or can be chained together into complete workflows. And these tools can be used for next generation sequencing, data processing, genotyping and variant discovery, variant filtering and evaluation. GATK also provides end to end workflows which are GATK best practice workflows which are tailored for specific use cases cases. And GATTK is very evolving and adapting to emerging sequencing technologies and has a very active user community and extensive amount of documentation There has never been a time or an occasion where I haven't found a solution to a problem on GATK forum or in their documentation. And so today for this demonstration, I decided to follow and demonstrate a best practice GATK workflow. In the previous video where we were discussing a VCF file format, we did take a look at how different types of variants were represented in a VCF file. So we took a look at how a snip or an insertion or a deletion or even a structural variant was represented in a VCF file. So variant is nothing, but it's an alteration in a DNA sequence and this alteration can be benign or pathogenic or of unknown significance. So one of the criteria to determine the various steps you would essentially perform in a variant calling analysis would be dependent on whether you want to call germline mutations or somatic mutations, whether the variant is a germline variant that you're interested or a somatic. Germline variants or mutations are the ones that are inherited from the parents via the germ cell, that is the sperm and the oocyte, whereas somatic mutations are the mutations that are acquired but are not inherited it so to give you an example of a somatic mutation, let's say a variant is present in a stem cell in an infant. So all the cells that are derived from the stem cell will have that variant, but that variant will not be present on all in all the other cells in that infant. In order to distinguish whether the mutation or a variant is germline or somatic, one has to sequence the tumor sample along with a matched normal. So consider an example of an individual with lung cancer. If we want to identify somatic mutations from this individual along with the tumor biopsy from the lung, one has to also get a matched normal from the blood. So germline mutations are also involved in pathogenesis. But it has been found that somatic variants are typically more involved in a disease and hence they are of special interest. Without having a matched normal, one cannot distinguish whether the variant is germline or somatic because every genome contains tens of thousands of mutations and and so a matched normal from the same donor is necessary to identify whether the mutation is germline or somatic. So for today's demonstration we are going to call germline variants and we are going to use germline Short variant Discovery Gatki Best Practice Workflow I will add the link to this workflow in the description section below. So basically this workflow has three steps. The first Step is data pre processing, which mainly involves getting analysis ready read reads and using these analysis ready reads, we call a variant calling algorithm to call variance, which is the second step. And once we have the variants called, we use these variants and filter them and annotate and predict the effects of these variants on genes. So today we are going to mainly focus on various steps that go in data pre processing and calling variants on analysis ready reads, using Gatkase Haplotype caller and calling raw snips and indels. Once you have called this these RAW variants, you can filter and annotate them. And there are various tools like SNPF or VEP that can be used which annotates and predict the effects of these variants. However, for today's demonstration we are going to just stick with the first two steps. So, talking about data pre processing steps, we start with raw unmapped reads in FASTQ files and we first perform quality control, making sure that there are no adapter sequences. If there are, then essentially we would trim them out and make sure that everything else with the reads checks out. I mean there are no major issues with the reads flagged. And once the quality control checks out, we move on to mapping these reads to a reference genome. And for that we need an algorithm to map the reads to the reference genome. And here we are going to use use BWA mem. So BWA stands for Burrows Wheeler Aligner and it's basically used to map reads against a reference genome. And BWA consists of three algorithms, BWA Backtrack, BWA SW and BWA mem. And each of these algorithms are designed for specific purposes. So BWA BackTrack is designed for Illumina sequence reads up to 100 base pairs. BWA, SW and Mesa are designed for sequences. For longer sequences ranging from 70 base pairs to 1 megabase pair, we are going to go with BWA MEM because it is recommended for high quality queries as it is faster and more accurate. In addition to that, it can also perform gapped alignment for identification of indels and it can effectively map paired entries. It can also perform soft clipping, I.e. it can clip ends if they do not match and it has a model moderate memory requirement. So when aligning raw reads to a reference genome using an aligner like bwame, the reads are. The aligned reads are stored in a SAM or BAM file and sometimes the SAM or the BAM file will miss certain information. The information we are talking about here is a read group. So basically the term read group refers to a set of read that are generated from a single run of a sequencing instrument. So all the reads that are derived from a single library preparation from a single biological sample that are run on a single lane on a flow cell will have the same read group. Now these read group tags are present in the header section of the SAM or BAM files. If you're not sure what the header section of the SAM or BAM file looks like, you can refer to the video where we have spoken about the SAM and BAM file format. I will add the link to that video in the description section below. So this read group tag which is is present in the header section of the SAM or BAM file, will have information on the ID or unique ID that is associated with the reads, the sequencing platform information, DNA preparation library id, sample names, and so on and so forth. So this tag can give us information on the IDs that are associated with a read, what is the sequencing platform used and what sample these reads are derived from. So at this point you would question as to why is it important to have this read group tag in the SAM or BAM file. Because when you use these aligned reads in the downstream processing, a lot of tools gather information from these read groups. So BWA MEM has a parameter called R that allows you to add this read group information in the resulting SAM or BAM files. I will be demonstrating demonstrating how to add this read group information when we move on to the demonstration part of this video. So once we have the aligned reads, next we want to flag the duplicate reads. So during sequencing process the same DNA fragments may be sequenced multiple times. So duplicate reads can arise during sample preparation step, that is Library construction during pcr. So it is important to eliminate these duplicate reads as these are not informative and can be considered as evidence for or against a variant. If we do not remove the duplicate reads, we can risk having over representation in our sequence of areas that are preferentially amplified during pcr. So these duplicate reads are not removed per se from the SAM or BAM files, but these are flagged as duplicate reads. So the second column in a SAM file is known as a bitwise flag and the flag value indicates indicates that the read is flagged as a duplicate in one of the videos. Previously we have spoken about SAM flags and what do they mean? So I will add the link to that video in the description section below. So once these reads are flagged as duplicates, they will then be ignored by all the downstream GATK tools. So for this step we will be using Mark Duplicate Spark tool from gatk, which Performs both the duplicate markings step as well as sorting the alignment files. So the last step in data pre processing is to recalibrate base quality scores. So variant calling algorithms rely heavily on the quality scores that are assigned to individual bases in each sequencing grade. This is because the quality score can tell us how much we trust that particular base at that location. So if we have a base call that has a low quality score, that means that we are not sure sure whether we actually see that base at that location or whether it could be something else. So in other words, we don't trust it as much as we trust the other base base calls that have high quality scores. So basically these quality scores serve as evidence that we have for or against a variant that exists at a particular site. Unfortunately, the scores produced by machines are subjected to various sources of system systemic technical errors which can lead to over or underestimation of these quality scores in the data. These errors are due to physics or chemistry of how sequencing reaction works and some are probably due to manufacturing flaws in the equipment. So base quality score recalibration is a process in which we apply machine learning to model these errors and adjust the quality scores accordingly. So this allows us to get get more accurate base qualities and it improves the accuracy of variant calls. So providing high level overview of how base quality score recalibration works. So we provide the GATK base recalibrator with a set of known variants and GATK base recalibrator analyzes all reads looking for mismatches between the read and the reference and skip those positions that are included in the known variance file. GADG based recalibrator computes statistics on the mismatches based on the reported quality score, the position in the read and the sequencing context and based on the statistics, an empirical quality score is assigned to each mismatch overriding the original reported quality score. Base quality score recalibration step is an optional step, but it's highly recommended in variant calling analysis. In the event that one does not have a node known set of variants available, then these can be generated by by calling variants first without using the base quality score recalibration step and filtering the variance to obtain high confidence set of variants and then using those variants as an input for the base quality score recalibration step. And this process is called as bootstrapping. So the second part of this pipeline is is variant discovery. So we apply a vein calling algorithm on analysis ready reads that we generated in the previous step. So we are going to use HAPLOTYPE caller to call variants from the reads. So the choice of variant calling algorithm has to be dependent on certain criteria like whether one wants to call germline variants or somatic variants, whether the organism we are working with is deployed, how many number of samples you're processing, and whether the samples are pooled or unpulled. So for today's demonstration we are we want to call germline variants and the organism we are working is deployed because it's it is a human sample. We are just working with one sample and hence the choice of variant calling algorithm today is haplotype caller. Haplotype collar can handle multiple samples, but it is not recommended to use haplotype caller when you are trying to analyze molecules more than 100 samples at a time. So for today's demonstration we are using data from Thousand Genomes project and for anyone who hasn't heard about this project, this project is an effort to produce an extensive catalog of human genetic variation that will support future medical research studies. So genomes of over thousand unidentified individuals from around the world were sequenced using next generation sequencing technologies and this DNA data is available for anyone to access and process. So for our demonstration we are using paired reads, whole genome sequencing reads from an individual who was a part of phase three of this project. I will be adding links to download this data in the description section below. So while going through the demonstration there are certain things that I will not go into the details of. So I'm assuming you have basic understandings of Linux system, I'm assuming you know how to write and execute a bash script, I'm assuming you understand Linux file systems and permissions, and I'm assuming you know how to download and install tools using command line and these are the requirements for today's demonstration. So the GAD key runs natively on most if not all flavors of Unix which includes Mac os, Linux and bsd. So it is possible to to get it running on some versions of Windows, but GATK does not provide support or instructions and even does not recommend to run GATK on a Windows system. So for all the Windows user, if you need to run GATK and process this then I will recommend you to get a virtual machine and install a Linux OS in your virtual machine. Because that's what I did when I was in my master's program we were recommended to get a virtual machine and install a Linux OS and perform all the bioinformatics analysis in Linux. So I will recommend you to do that as well. You will need Java 1.8 and for the tools you will need CatK, BWA, FastQC, SAM tools, and MultiQC for today's demonstration. So before we start writing the script, let us set up our directory by creating some folders in it. So this is my project working directory and you can see that it is empty right now. So I'm going to create some folders within it. So folders to hold, align reads, to hold reads, to hold scripts, to hold results, and to hold some data that we generate, the intermediate files that we generate during the pipeline. So now you can see that we have five folders here and these folders are empty right now. But we are soon going to populate these folders with the relevant data within in it. So now moving on to the script. So this is our script and we first start by downloading the data. So within our project working directory you can see that we have subfolders. So this helps to keep the data more organized. So I want to download the reads in the reads folder. So I'm going to provide the path to where I want my files to be downloaded with WGET when I'm downloading it. So I will be providing the path, the prefix for wget. So that's where I want my reads to be downloaded in the reads folder. So now let's save this and go to the terminal and let's run the script. So our script is present within the script folder. So let's go to the script folder. And now let's run. So we might need to give the permission for it to run first. And now let's run the script. So these files can take a couple of minutes as these are larger files. So I will be adding the links to download these files in the description section below. We have finished downloading the reads and the reads are present within the reads folder. So let us go back to our script script after downloading the reads and we are building this pipeline with just one sample. But let's say if you're running and you're trying to create a variant calling pipeline for multiple samples, in that case you will require some files that your pipeline will be reusing multiple times. So these files need to be generated only once. But these are the files that are required by GATK programs. You so some of the files that you will be required to download and create. So first we start by downloading the reference files. This is the reference genome file. And this, the. This is important because we are going to align our reads to the reference genome. So basically we require a Reference genome fasta file. So we will be downloading the reference fasta file. Followed by that we will be creating an index of the reference anna dictionary of the reference and these are some of the intermediate files that GATGPT programs require. And lastly, if you remember, we are going to in the data pre processing step we are going to perform something called as base quality score recalibration. And that step, if you recall, requires a file with known variants. So we are also going to download these files. So we're going to make sure that all these files are downloaded and are ready before we start creating a pipeline so that none of the downstream GATK tools will throw show any error regarding not being able to find these files. So let us write the commands to download and create these files. So we are going to use wget again to download the reference fasta file and we are going to download it at within this folder. So once this gets downloaded, I would also want to unzip or gunzip the file because this file is GZ compressed, so I need to to uncompress this file. So this file will be present at this location and then the name of the file would be this. So I want to gunzip it. So after gunzipping I will have HG38FA FASTA without the GZ extension. Next let us create the index file for the fasta file and for indexing a fasta file we will be using SAM tools faidx command and this is the file that we would want to index. The output of this index file would be dot FAI the same file name with dot FAI extension and it will be created in the same folder. Next we also need to create a dictionary file and for that we will be using get K command create sequence direct dictionary and it requires the reference file. So this is our reference file and the output would be the dictionary file with dot DICT extension. So instead of dot FA it would be dot DICT and again this file needs to be present in the same folder where we have the index file as well as the reference file. So all these files need to be created in one folder. And lastly let us also download the known sites variant sites for base quality score recalibration steps. So let us again use WGET to download this file. This is the file containing the sites for the known variants. And again I want to download this in the same location where I have my reference files. And let us also download the index file for the VCF because that will also be required by the programs and Again, it will go in the same folder. So now you save the script and run the script to create these files in your folders. Again, I already have downloaded and created these files, so I will not be running this section of the code again. But since you are going to build the pipeline for the first time, you will require these files. And these files will only have to be generated once, and they can be used and reused again with the other samples that you are trying to process. Once you have downloaded all your data and you have generated these files, you will not be required to run this block of code again. So I will be putting an if and then block to prevent this from re running again and again every time you run it with separate samples, or let's say if you're running multiple samples, you don't want these files to be recreated or the data to be downloaded every time a different sample is processed. So this saves much of your time and memory in your system. As for the next steps in the effort to make the code appear clean, I'm going to save the parts to to certain files and folders in the variables because this will save us from rewriting the parts again in the script, making the script look very messy. So just to make it more organized and easier to understand, I'm just going to save the parts to some files and folders within the within some variables. So we're going to save the parts to reference and known side files within some variables. So let us start by providing the parts to these files. So this is the reference, this is the path, and the file with the complete path to reference file. So I'm going to save it in a variable called reference. Then we have known sites files as well, and the path for it is this. So this is where the file is saved with the complete path. Now let us also provide the path to the respective folders within our project working directory so that we will not have to write the path again and again. So just going to the terminal and getting the project working directory. And we also do it for reads, We do it for results, And finally for data. This is a very common practice of saving the paths to the variable when writing bash scripts, because this just helps to make the code appear clean cleaner. So now let us begin with the first step that is quality control of our reads, assessing the quality of our reads, making sure that we have good quality reads to begin with and use them for further processing. So we use a tool called FastQC and we provide the path to the reads. And as you know, the Path to the reads is saved in the reeds variable, so we just interpolate that variable and we need to provide the name of the reads. So so since these are paired end reads we will have two files. The first file contains the forward set forward pairs. So we are just going to provide the file name and output it again to the reads folder and we are going to repeat the same process for the second set of reads, that is the reverse reads. Just to make a note here that the tools that I am using the executable of these tools, the path tool to these executables are in my bash RCO bash profile, meaning that these tools are in my path and hence I do not need to provide the complete path to the executable. But if you haven't added the executable to the path then I will recommend you to write the entire path to the executable for this tool. Otherwise it will just throw you an error that it is not able to find this tool if you just just directly try to execute it. So just making a note of that here. So now that we have saved this, let us go back to the terminal and go to our script folder and run the script. Now that this is finished running let us take a look at the folder contents and we can see that in addition to the sequencing reads we have files that are generated ending in.HTML and.zip. so I want to open the files ending in HTML to see the results of the fastqc because these are the FASTQC reports so I'm just going to open them looking at the forward reads first. At first glance this sample these set of reads rather fail for the per tile sequence quality and has been flagged for GC content. But the rest of the metrics are seems to check out. So looking at the basic statistics the the read, the size of the reads are 100 base pairs and the GC content is 40%. So I'm mainly going to check a couple of things starting with whether there are any poor quality bases present in the reads. And here it seems that all the reads check out the per base sequence quality. In addition to that this these set of reads are flat stacked for GC content but did not fail the GC content so I'm not going to be very much worried about that. And lastly I'm also going to check whether there is any adapter content present because if there are adapters present then I will have to trim those out before I align them. Another note that if I do see poor quality bases had it been that this had been flagged or it failed for poor quality bases, I would still not trim out the poor quality bases because as the aligner that I'm using it performs soft clipping, so it basically clips off the bases that do not match. And also I would risk shortening the length of my reeds so that would result in multiple incorrect mappings. So I would not go for trimming the read reads, poor quality basis and risk having shorter reads. So looking at this test report, this assessment report, I'm not going to perform any trimming. And let's go to the reverse reads. Let's go to the second set of reads and it seems that at first glance these set of reads are flagged for per base sequence content and GC content, but the rest of the metrics check out so the per base sequence quality looks great. The GC content is flagged again, but I'm not worried about that very much. Also I do not see any additional peaks here, so I'm not going to be concerned about it. And lastly I do not see any adapter content. So basically overall this report suggests that I do not have to perform trimming mainly to remove any adapter sequences or trimming the poor quality bases. Coming back to the script, I'm just going to add a little comment here saying no trimming required as quality looks okay. Now that the quality of the reads look great, we can move on to the next step where we map these reads to a reference genome and we use BWA aligner to do do that. So BWA also indexes reference because it allows it to be more efficient to search the genome while performing the alignment. So we use BWA index function to generate the index for the reference. And since the reference file with the complete path is stored in the reference variable, we are just going to provide the variable here so that the script can interpolate the path and the file file here and followed by that we can also provide the command to run the alignment. So we are going to use BWAMEM as we previously discussed and we are going to use four threads and since we know that the resulting align reads, the SAM or the BAM files will be missing a read group, we will be providing the read group information here. So we provide the read group information in double quotes where the regroup tag starts with G and this field is tab delimited. So we separate the fields by tab. So the first value would be id and here every read would be having a unique id. But since we just have one sample, I'm just going to provide the name of the sample to be the id Followed by that we provide platform information, which is Illumina and followed by that we provide the sample information. And again, since I'm just running one sample, I'm just going to repeat the sample name here. Followed by that we provide the reference file and provide the reads. So the reads are present in the reads folder and these are the names of the file. And I want to output the align reads in aligned reads folder. So I'm going to provide the sample name and I want the SAM file. So we are using BWAMEM algorithm, we are defining the read group, we are providing the reference, the reads and outputting SAM file. And since we've already done the FASTQC for our reads, I'm not going to rerun it again and hence I am just commenting it out. So now let's save the script and go to the terminal, Go to the script folder and let's run the script. Now that the alignment step is finished running, I first want to show you the output files that are generated at each step in the alignment. So first we index the reference. So let us take a look at what files that are created after the indexing. So if you recall, the first five files were downloaded and generated at the beginning of the script. The next five files are what BWA index function generated. So these were the files essentially that BWA used for alignment. Next let's take a look at the folders. So in the aligned Reads folder we can see that the aligned reads are are present in the form of a dot SAM file. So let us take a look at these reads. So we can use SAM tools view function to view the SAM file. And taking a look at the SAM file we can see see that there is there are various kinds of information present for the alignments. We have previously discussed what a SAM file is and how to make sense of this information. Here, in addition to viewing the alignments, I also want to take a look at the flag stat. The output of the flagster provides us with counts for various categories present in the SAM bitwise flag field. So can you. So here we can get like bunch of information as to how many pairs of reads are present, how many mapped, how many are duplicates, how many are supplementary. And for now you can see that there are no reads mapped as duplicates. And we previously discussed that during generation of the sequencing reads, there could be steps during library preparation where we could have duplicate reads. So it's important to flag these these reads. So the Next step is to mark duplicate reads so that and flag them so that the following programs and the downstream analysis, the programs in GATK would identify these flags and would ignore these reads. So we will utilize gatcase function mark duplicate Spark, which performs both the flagging of the duplicate reads as well as sorting the SAM files. So the input here would be align reads. So basically this file. So I'm going to provide the complete path as well as the file name. And I want to output the deduplicated sorted reads again in the align reads folder. And I'm going to name this as the sample name underscore sorted sorted underscore deduplicated reads dot bam. Again, we have previously run these commands. I'm just going to comment it out because I don't want to run them again. Let's save the script and go to our terminal. And now go to the scripts folder and let's run the script. The mock duplicates step has finished running. So let us take a look at the folder contents. So inside the aligned reads you can see that we have additional files. So we have deduplicated reads, bam and additional files. So let us take a look at the deduplicated reads. So let us go into the aligned reads and take a look at flagstat output for the duplicated reads. The flagstat output indicates that there are these many reads that are marked as duplicates. And the subsequent programs in CATK will ignore these many reads in the following analysis. Now that the duplicate reads have been successfully flagged, the next step in data pre processing is to correct for base quality scores to perform base quality score code recalibration. So if you recall, there are essentially two steps where in the first step we build a machine learning model using the known variance and then use those that model to adjust for the base quality score. So to build the model we will be using GAD K Bayes recalibrator function. The input would be the deduplicated sorted align reads. The reference point would be the reference file. We provide known sites and this file is stored in a variable called known sites and it is going to output a table. So that table will. I want the table to be outputted in the data folder and we will call it recall data table. Once the recall data table is created created, we will use that to adjust the base quality score. So to adjust the base quality score we will be using a function in GATK called apply bqsr. The input again would be the align read sorted deduplicated reads that we generated previously. The reference would be the reference file stored in the ref variable. We will provide the BQSR recall file that we generated in the previous step and the output would again be in the align reads directory. And I will copy the name of the file generated in the previous step, just adding bqsr. So we will be able to distinguish that these are the reads with base quality score recalibrated. So again commenting out mark duplicates step because I do not want to rerun it again. Let's save this and go back to the terminal. Let's go to the script folder and run the script. Once both these steps finish running, that is creating the model and using the model to correct for the base quality scores. Let us check the folder contents to see see what files are created. So in the data table we create the recall data table and this table is essentially used to correct and recalibrate the base quality scores. So the final alignments with base code recalibrator are present in this file and the index for this BAM file is also generated. So basically the reads that are present in this file are analysis ready reads. And these read reads can be used to call variants from. But before going on to that step, there is another step that we can perform, a sort of a QC where we can collect some alignment and insert size metrics. So to collect alignment metrics we can use GAD case function collect alignment summary metrics. The reference used is saved in the reference variable. The input would be the base quality score recalibrated reads and the output would be in the result or rather in the alignment aligned reads folder and we'll call it alignment metrics. To collect insert size metrics, we can use function collect insert size metrics. The input here would be the base the same the base quality score recalibrated align reads. The output would be again, so here I forgot to put the dollar sign. So I want the output to be in the aligned reads folder as well. And I want the output to be in insert size metrics. And I also want histogram and I want it again to be in the align reads folder and we can call it insert size histogram. So basically this histogram would provide us with the distribution of insert sizes across the reads. And before we run this, let us just comment out the previous commands. Save this and let us go to the terminal. Let's go to the script folder and let us run the script. So now that this has finished running, let us take a look at the folder contents. So here you can see that we have alignment metrics, we have insert size metrics, and we have insert size histogram. So taking a look at the histogram here, we can get a distribution of insert sizes for the reads across the reads. So basically this provides us a way to validate the library construction. In addition to that, we can also create a multiQC report to assess the alignment metrics and the insert size metrics. So let us go to the terminal and go to Align reads and use MultiQC to generate a report. Just a quick note before I move on to the MultiQC report. MultiQC is basically a tool that aggregates and summarizes results from the text files or the log files that are generated from other bioinformatics tools. So here we generated alignment metrics and insert metrics from using Cat K functions and the information present in each of these files can be summarized in a report using using a tool like MultiQC. Also, while executing MultiQC, we just provided a dot meaning that to we asked MultiQC to search for the text files or the log files generated by other programs in the current folder, that is in the align reads folder. So in the align read folder it could pick up these two text files I.e. alignment metrics and insert metrics. So basically it's going to create a report using these text2text files. So multiqc generates an HTML report and when we open this, So basically this is how it looks like. So it provides us with bunch of metrics to assess the quality of the alignments as well as the quality of the reads and the number of reads that passed the quality thresholds. And as well as provides us with the way to validate the library construction. So here we get the number of reads that are aligned the insert sizes. We also get the number of reads which are unaligned. We can also get mean read length. We get the insert size. So basically this can be looked at as as a post alignment QC and can help us evaluate the alignment as well well as the library construction. Now that we have finished pre data pre processing steps and we have also collected post alignment qc, we are ready to call variants on our analysis ready reads. So let us use cat k haplotype caller. The age is capital provided with the reference file. The input would be analysis ready reads that is the base quality score, recalibrated reads and the output would be in the result folder, results folder and these are raw variants.vcf I'm going to comment out the previous commands. We save this and we move on to the terminal. So let us go to the Scripts folder and run script. So this steps take usually a longer time to run. So this took around three hours for me and I have access to 16 gigs of RAM. So I'm not going to run this step because I previously run it and I have the results. So I'm just going to show you the output files. So in my results folder I have run rawvarins vcf and I have indexed VCF file which is generated by Haplotype Caller. So now we have the raw variants. So let us take a look at these raw variants. So we can see that there are a lot of variants called and this is a VCF file and we have previously discussed, discussed the VCF file and how data is organized in this file. Now, going back to the script, now that we have called the raw variants, I want the snips to be present in one file and I want the indels to be present in another file. So I'm just going to separate out the variants. So to select specific variants, I'm going to use a function get k function called select variance the reference, I'm going to provide it the reference file. We also need to provide the variance file which is present in the results folder and we need to select variant types. So first I want to extract snips and I want to output it again in the results folder. And I want to repeat the same step for the indels as well. So I'm just going to copy the same command and instead of the snip in the select app, I'm going to type indel and I'm going to change the name of the output file. I'm going to comment out the previous comment and I'm going to save this. And going back to the terminal, let us go to the script folder and let us run the script. So this is finished running and let us take a look at the folder contents again and we can see that we have raw snips, VCF and the indexes and the raw indels, VCF and their index in the result files. So let us take a look at these files. So taking a look at raw SNPs, basically it indicates these are the SNPs that are called from this individual. And taking a look at raw indels, these are the raw indels, these are the indels that are called from this individual. So these are the variants that is a SNP and the indels that are called from this individual. Using the Haplotype Caller variant calling algorithm from GATK and we follow the GATK Best practice workflow. So typically once you have called the variants you would essentially filter out the low quality variants and retain the high confidence variance and you would annotate these variants and predict their effects. There are various tools like Snipe and BE Peak that can help you annotate and predict the effects of these variants. So that brings me to the the end of this video. I hope this video was helpful and informative. I try to keep this lesser on the technical side and more on the purpose behind performing these steps. Again, these best practice workflow are not supposed to be used as is. There are a lot of these pipelines are always has to be tailored according to your data set and your requirements. This was just an idea on what are the various steps that are typically performed in a variant calling analysis and the rationale behind performing these steps. So I will be adding the various links that I have spoken throughout this video in the Description section below and I hope you find this video informative. If you did, please make sure you hit the subscribe button. Like the video, share it and leave your comments under the comment section until next time. See you.