Dataset

The data used for validating the pipeline was extracted from the Sequence Read Archive (SRA) toolkit. The samples were extracted from a study from Department of Molecular Biology of Cancer, Medical University of Lodz, Poland (ID: PRJEB27635). The study contains whole exome sequencing data from 8 different melanoma patients according to the following description: Melanoma specimens were obtained during surgical interventions, and grown in vitro in serum-free EGF/bFGF-containing medium as described previously (Sztiller-Sikorska et al., 2012 Melanoma Research). DNA was extracted using DNeasy Blood & Tissue Kit from Qiagen according to manufacturer protocol. DNA samples were subjected to Nextera tagmentation followed by index PCR, first hybridization, capture, washes, elution, second hybridization, capture, washes, elution and enrichment of DNA fragments. The samples were subjected to two rounds of exome enrichment, followed by real-time PCR quantification (KAPA SYBR FAST qPCR Kit) prior to sequencing.

The samples were sequenced using NextSeq 500 sequencing of the cell DNA after in-vitro culture. The different samples extracted can be found in the table below:

Table 1 - Summary of the data contained in study and basic metrics PRJEB27635

The datasets were downloaded in FASTQ format using fastq-dump tool in the SRA toolkit (v. 2.9.1). Documentation and details on how to use this packages can be found here. In general, the data from the study appears to show good general indices, with a good overall GC content, in the lower end of the expected 45-50% range of typical exome sequences (MacKinnon). At the same time, the large base size of around 100 times bigger than a typical human exome, showing good coverage of the reads.

Limitations of the Dataset

Before starting the analysis there are a few notes to point out about the use of this dataset. First, the accuracy of the genomic analysis highly depends on decoding the and subtracting the noise and artifacts generated during the sequencing protocol (Wang et al.). Furthermore, other characteristics, such as coverage, index adapters, number of flow-cells among other parameters are key to avoid false positives. The samples contained in the study PRJEB27635, have limited information about these characteristics, which only allow to make assumptions based on the observations of other parameters. This makes the validation process shown in this part more conversome.

Furthermore, In order to precisely find all the mutations contained in the genome and avoid any bias originated by the sequencing protocol or the genomic analysis, a matched normal sample from the patient is required. The study PRJEB27635 does not contain any normal samples and hence the tumor and normal samples cannot be directly compared to detect the mutations. However, recent developments in the tools allow the use of only the tumor sample (Teer et al.)

Quality Assessment - FastQC Analysis

The data was pre-analysed using FastQC (Version 0.11.7). The program outputs a summary of different metrics of the reads and signals with a flag the scores that are suboptimal for genomic analysis. In general the results show good quality of the reads, with overall good per base sequence quality and Per tile sequence quality for all the patient’s data. Furthermore, the per sequence quality score follows the typical trend for good illumina sequence data as seen in Figure 1. This implies that the data was collected accordingly and it represents the overall field of view of the targeted data. However, FastQC report shows three metrics that are consistently poor for all the samples in the study. These are per base sequence content.

All samples contained an irregular profile of per base sequence content. The per base sequence content represents the relative amount of each base in your genome. Normally, in correctly sequenced data this value should be maintained constant throughout the genome. In the data obtained from the study the difference between A and T, or G and C is greater than 20% in certain position especially concentrated at the end of the library. This results can have several implications in the data to analyse. Often, this errors are caused by the presence of overrepresented sequences such as sequencing adapter dimers use to index and sort the Illumina reads, this can often bias the data and cause unrealistic results (Andrews). This suggest that the adapter primers are contained in the sequences. However, by looking at the Overrepresented sequences metric in report none of the datasets seem to contain any heavily repeated sequences. On the other hand, the low score of the per base sequence content can be attributed to a aggressive trimming of adapters that can introduce bias in the library. This typically occurs at the end of the reads (Andrews). The per base sequence content for the study reads as seen in Figure 2, seem to be extensively trimmed as the sequence content at the end of the reads deviate from the general trend. This is further confirm by the GC content of the data which shows a deviation from the normal distribution of more than 15% of the reads (Andrews).

Sequence Duplication levels were also found to be abnormal, since non-unique sequences make up more than 20% of the total library for all the samples. This can be appreciated by high duplication peaks in the center of Figure 3. FastQC also calculates the overall loss of the sequence for the library after removal of the duplicated sequences, found on the title of Figure 3. For the sample ERR2683866, the non-duplicated sequences represent 65.2% of the total library. This might be due to specific PCR enrichment during the sequencing process (Andrews).

Based on the output of these reports the pipeline was optimized to ensure that any possible adapters that might affect GC sequence content or Per base sequence content are removed, using BBmap or Trimmomatic. Furthermore, given the high content of duplicate reads its necessary to ensure that the library is filtered to remove any bias information that can affect the final outcome. Duplicate reads were filtered using MarkDuplicates package in picard tools.

Sequencing Adapter Cut and Quality Control - Trimmomatic/BBmap

The sequences were analysed using BBmap (v38.36) and Trimmomatic (v0.38), to remove any sequencing adapters related to the reads (Bolger, Lohse and Usadel 2114-2120). Usually sequence data contains adapters used for indexing and sorting the reads in the genome. Normally, is ideal to use the specific adapter-sets according to the technique used to perform the sequencing. TruSeq is commonly used for performing this analysis, however, the original method of extraction of the libraries of study PRJEB27635, are not known which makes the task of assessing the quality of reads a challenge. In order to ensure that the analysis is complete, Ginga implements both Trimmomatic and BBmap trimming tools.

All samples were analysed using Trimmomatic (v0.38). To remove the adapters the Trimmomatic was run using the single end mode (SE) and IlluminaClip was initially performed using the standard parameters (seedMismatches = 2, palindromeClipThreshold = 30 and simpleClipThreshold = 10) and the universal adapters contained in the TruSeq v3 single end adapters which can be found in the table below (Table 1). With the above parameters, no sequences were clipped from the raw initial sequences. The method was repeated, in order to optimize the parameters. seedMismatches was maintained at 2 bp since TruSeq v3 adapters are particularly short and enriching the bias towards mismatched sequences can lead to deletion of exonic information (Bolger, Lohse and Usadel 2114-2120). The other quality control options of Trimmomatic were not used, since the FastQC report of the library shows good overall base per quality and distribution, however they are recommended.

To ensure that the sequences were not containing no other bias, other typical adapters were used, including NextEra [ID:5-6] and TrueSeq [ID:3-4]. AmpliSeq adapter trimming, typically used for creating large cancer panels was also tested [ID:7]. When using adapters 5 and 6, the 265 sequences were successfully clipped, which accounts for below the 0.01% of reads in the library, while the other adapters [ID:3,4,6] where unsuccessful. This is unusually low, and might be due to sequences that are similar to the kmer adapters in the exome.

Table 1 - List of universal and index adapters used against the samples.

Alternatively, to Trimmomatic, Ginga also supports the use of BBmap tools. Using the Decontamination Using Kmers (BBDuk) package it is possible to search the library for the most common kmer adapters and perform trimming automatically. The file containing the list of adapters used can be found in Ginga’s GitHub. BBDuk standard settings were used (ktrim=r k=23 mink=11 hdist=1). Adapters based on pair overlap were also trimmed. The results showed that the algorithm was able to successfully trim 307 reads (0.00%) and 521667 bases (0.01%) for sample ERR2683866, which closely matches the output of trimmomatic.

Reference Alignment, Sort and Indexing - BWA and Samtools

All the samples were aligned with the hg19 human reference genome using Burrows-Wheeler Aligner, BWA (v.0.7.17) mem command . Although, the new updated GRCh38 reference sequence is available, some of the current tools of the pipeline do not support this reference genome. This new genome annotated reference genome will be supported in future versions of Ginga. hg19 is based on the b37 version, with some different updated patches and nomenclature (Auwera). The reference genome was obtained from the Broad Institute FTP Server (McKenna et al. 1297-1303). The alignment was performed using standard conditions, including no clipping and discard options. The output of BWA is a SAM file containing the aligned sequence.

Once the sequences were aligned the output SAM files were sorted using Samtools (v 1.7) sort(Li et al. 2078-2079).The sequences were sorted with no compression level and by default coordinates in the original aligned SAM file. The sorted file was outputted in the compressed BAM file format to reduce the size of the library and simplify the analysis, later in the pipeline. The headers of the processed BAM file can be found in the table below (Table 2). Additionally to this main headers the file contains between 39 clone contigs that were not able to be placed in any specific part of the gene (ChrUN) and 20 additional chr_random tables that contain data related to sequence that is known to be in a particular chromosome, but could not be reliably ordered within the current sequence (Ken et al.). Furthermore, the remaining 9 tables that were not correctly aligned correspond to different haplotypes in the chromosome. Normally, the reference sequence can only display a single haplotype per chromosome, so alternative haplotypes are contained in this file (Ex. chr6_apd_hap) (Ken et al.).

Table 2 - Header of the generated BAM file for sample ERR2683866

                  @HD     VN:1.5  SO:coordinate
                  @SQ     SN:chrM LN:16571
                  @SQ     SN:chr1 LN:249250621
                  @SQ     SN:chr2 LN:243199373
                  @SQ     SN:chr3 LN:198022430
                  @SQ     SN:chr4 LN:191154276
                  @SQ     SN:chr5 LN:180915260
                  @SQ     SN:chr6 LN:171115067
                  @SQ     SN:chr7 LN:159138663
                  @SQ     SN:chr8 LN:146364022
                  @SQ     SN:chr9 LN:141213431
                  @SQ     SN:chr10 LN:135534747
                  @SQ     SN:chr11 LN:135006516
                  @SQ     SN:chr12 LN:133851895
                  @SQ     SN:chr13 LN:115169878
                  @SQ     SN:chr14 LN:107349540
                  @SQ     SN:chr15 LN:102531392
                  @SQ     SN:chr16 LN:90354753
                  @SQ     SN:chr17 LN:81195210
                  @SQ     SN:chr18 LN:78077248
                  @SQ     SN:chr19 LN:59128983
                  @SQ     SN:chr20 LN:63025520
                  @SQ     SN:chr21 LN:48129895
                  @SQ     SN:chr22 LN:51304566
                  @SQ     SN:chrX LN:155270560
                  @SQ     SN:chrY LN:59373566


                

An index file was generated for each of the sample using the index function in Samtools. The index file allow aligners to narrow down the potential origin of a query sequence within the genome, saving both time and memory. The index file typically is a .bai file, with the same name as the alined BAM file. This step is key for future package alingers to asses the sequences.

Remove Duplicate Reads - Picard

After assessing the quality of the samples with FastQC, the reports show that more than 20% of the total sample is composed by duplicated and non-unique sequences as seen in Figure 3 (Shlee). The duplicated sequences were search and deleted from the libraries using Picard Tools (v 2.18.14) MarkDuplicates program. The remove sequencing duplicates option was tagged to output the file with the duplicated sequences already removed. Alternatively, it is possible to flag the duplicates instead of removing them from the file by using <-drf> DuplicateRead option (Shlee). The duplication metrics of this analysis were also compiled and can be found in the Ginga Github. A summary of this metrics can be found in Table 3. The average percentage of duplication for the library was 36% which correlates to the previous information reported in the FastQC report. For more information about correct usage of MarkDuplicates please refer to this link.

Table 3 - Duplication Metrics for all the samples in the study.

Base quality assessment - GATK

Variant calling algorithms are heavily sensitive to the quality scores of the sequences. The base quality scores are often subjected to various technical errors due to the sequencing technique. For this reason base recalibration is key to ensure that the quality scores are accurate and that the variant calling analysis is accurate. Base recalibration was perform with BaseRecalibrator inside GATK package (v 4.0.9.0) (McKenna et al. 1297-1303).

Before recalibrating the bases, it is necessary to ensure that the samples must be sorted, indexed and contain the appropriate read groups. Read groups are fields contained in the header of the libraries (@RG) and the represent a set of runs that were obtained in the same run. Due to the unknown specification of the dataset used, it is not possible to determine the appropriate read group, however, for the validation of the pipeline it was decided that each sample was part of a unique read group. To view if the sample contains any read group, Samtools view command can be used by specifying the header group of interest (ex. -H <sample.bam> | grep ‘@RG’). For more information about read-groups and correct usage please refer to (Van deer Auwera).

In addition to the sample input BAM file, BaseRecalibrator has the option to add more additional fields labeled as known sites to calibrate the analysis. The know sites refer to any database of known polymorphic sites for the specific reference sequence. This include a list of know SNPs and Indels sites within the reference genome. During the creation of the model and analysis of reads, this known sites will be avoided and the program assumes that all reference mismatches seen outside this sites are therefore errors and indicative of poor base quality. There are three main files recommended for a correct BaseCalibrator analysis: the dbSNP sites, Mills and 1000 Genome Project gold standard indels and the 1000 Genome Project standard Indels. This datasets can be obtained from the Broad Institute FTP server and are specific for each human reference genome. For this study hg19 reference genome and known sites were used. Base recalibration was performed with standard settings. The output of the analysis is a GATKReport file containing the various recalibration tables.(McKenna et al. 1297-1303) For more information about the format and characteristics of GATKReport format please refer to this link. A summary of the recalibration table by read group can be found in Table 6. The full reports can be found in Ginga Github.

Table 6 - Quality Base Scores for the different samples in the study.

The output BAM file from base recalibration can be used directly to call the somatic mutations in the samples. Additional information might additional be trimmed or flag in order to avoid bias in the final result.

Mutation Calling - Mutect2

Mutation calling is one of the core functions of the Ginga pipeline. The somatic mutations found in the samples are responsible for the modifications in the final translated proteins and effectively the neoantigen sequences that Ginga aims to detect. Unfortunately, due to the lack of ethical approval to access other patient’s dataset and the concerns regarding the data protection of whole exome sequences, the results of this validation were perform without a reference matched normal sample. The lack of a matched normal reduces the mutation calling accuracy of this validation results and requires extensive optimization of the parameters. However, this analysis is now supported in the latest version of Mutect2 (v 4.0.9.0) and can be performed by running the single tumor sample mode in addition to using databases of known polymorphic sites (Kalatskaya et al.).

The mutation calling was done with Mutect2 in GATK (v 4.0.9.0), which currently is under development (Beta Version). The samples were initially analysed in single tumor sample mode against the hg19 human reference genome with no germline resource or PON. This first results identify an unrealistic number of mutations, for example a total of 198,760 mutations were identified for sample ERR2683866. These results (Table 7) are expected as no matched normal sample was used to perform the analysis and the known sites were not included.

Table 7 - First 10 mutations found in the output VCF file of mutect after analysis sample ERR2683866 using Mutect2.

                  CHROM	POS	ID	REF	ALT	QUAL	FILTER
                  chr1	830181	.	A	G	.	.
                  chr1	871334	.	G	T	.	.
                  chr1	876499	.	A	G	.	.
                  chr1	877831	.	T	C	.	.
                  chr1	878314	.	G	C	.	.
                  chr1	880238	.	A	G	.	.
                  chr1	881627	.	G	A	.	.
                  chr1	883625	.	A	G	.	.
                  chr1	884101	.	A	C	.	.
                  chr1	887560	.	A	C	.	.

                

The analysis was then optimized by using the germline resources from the gnomad, in order to consider the known sites of exome SNPs and Indels in the analysis.

Filter Mutations - GATK

Due to the lack of matched normal to perform the mutation calling analysis, there is a higher risk to select false positives in the pipeline. Some of these artifacts are generated by the presence of cross-sample contamination in the original FASTQ file and are due to sequencing technique. For this reason, it is important to reduce the large pool of mutations generated and filter them accordingly. This can be done by using GetPileupSummaries, CalculateContamination and FilterMutectCalls from GATK (v 4.0.9.0) package (McKenna et al. 1297-1303).

The data was first processed GetPileupSummaries, to summarize the reads that support reference genome and list of known allele sites. This can estimate the overall cross-sample contamination. The BAM files used prior to mutation calling were assessed against the hg19 human reference genome and two sets of known sites from the Gnomad Library only containing biallelic variants. The whole exome known sites (v 2.0.2) was filtered using GATK (v 4.0.9.0), SelectVariants to restrict the alleleles to only biallelic. The maximum population of allele frequency was set to 0.2 and the minimum mapping quality to 50. Furthermore, due to the interest in search on the whole exome sequence for the neoantigen variants, the whole Gnomad biallelic variants were selected as the interval of the analys. The output consists in a table containing 6 different columns: contig, position, ref_count, alt_count, other_alt_count, allele_frequency. The alt_count field contains the reads that are supported by the germline resource.

The pileup summaries table was then used to calculate the fraction of cross-sample contamination in the tumor reads and output the contamination table that can be used to filter the mutations. The resulting table contains the estimates for contamination and error in the sample. Furthermore, CalculateContamination also supports an alternate mode which uses the matched normal sample to improve on the accuracy of the analysis.The results show that the contamination levels for the samples are relatively low.

Table 8 - Contamination Table for sample ERR2683866.

The mutations outputted by Mutect2 were filtered using FilterMutectCalls. This applies a preset of thresholds and filters to remove false positives from the data. It uses the contamination table as an additional input. FilterMutectCalls outputs a VCF and index file, containing the true positives with the PASS label and the false positives with the FILTER field.

Gene-Based Annotation - Annovar

The gene-based annotation was performed using Annovar (v. 2018Apr16) annotate_variation.pl script. Gene-based annotation was used to annotate the mutations of the discover in the samples and the coding_change.pl was change to translate the annotated sequences into the corresponding coding proteins. Additionally, filter-based mutation was used to filter mutations from the Cosmic Database and the International Cancer Genome Consortium (ICGC) database.

Annovar uses the predetermined file format VCF4 for the inputs. The annotated mutations from Mutect2 or other GATK packages commonly found in VCF format can be converted using the perl code convert2annovar.pl. This tool adds three additional rows to the VCF format, zygosity status, genotype quality and read depth. The genotype calling files from the samples extracted after filtering the contamination with FilterMutectCalls were converted to VCF4 format.

Before annotating the samples with annovar, the FASTA sequences for all annotated transcripts in RefSeq Gene hg19, were downloaded using the built-in --downdb argument in the annotated_variation.pl package. This argument allows you to access files directly from the UCSC Genome Browser Annotation Database. Furthermore, the reference gRNA files were created using the retrieve_seq_from_fasta.pl package in Annovar. The VCF4 were annotated using the core function of Annovar, the annotate_variation.pl program. Standard parameters were selected and the hg19 refGene was used as a reference. The annotated variations are outputted in two different VCF files. The variant_function VCF file, contains the information regarding the nature of the variation and the position in the exome. This file summarises if the variant is exonic, splicing, or if it hit a non-coding RNA genes, as well as the gene containing the mutation. The second file, exonic_variant_function VCF file, contains the amino acid changes as a result of the exonic variant and furthermore identifies the type of functional consequences caused by the variant (frameshift insertion, deletion,frameshift block substitution...). The exonic_variant_function VCF annotations, were translated into protein sequences through the coding_change.pl tool in Annovar. This outputs a FASTA format containing the wild type and mutated protein sequences under a header containing the original mutation, gene affected and amino acid sequence change. All the annotation results can be found in Ginga’s Github.

Sort Annotated Proteins - neoExtract

The proteins annotated in Annovar contain both the mutated sequence as well as the wild type protein. Moreover, since neoantigen sequences are typically 8 to 11 amino acids long, the major part of protein sequence which doesn’t contain any functional variation associated with the cancer mutation of the sample, has no implications on the end library of neoantigens. For determining the final neoantigens, only a short amino acid sequence around the mutation is required, in the case of SNVs. neoExtract was used to sort the annotated proteins, removing the wildtype from the Annovar FA output file, and only retains a short amino acid sequence of specified length around the mutation. The coding change exonic variant function, was processed by neoExtract, using standard parameters. The output FASTA file contains all the possible combinations of neoantigen candidates derived from the mutated proteins. neoExtract can also process other functional consequences of the variant such as insertions, deletions and frameshift. The first 10 results from the

MHC-I-Peptide binding Affinity - NetMHC

In order to determine if the neontaingen candidates can be presented in the surface of of antigen presenting cells such as dendritic cells, which are the precursors of the immune response, the neoantigens need to have a high binding affinity to the MHC-I complexes on the cell surface. This MHC-I complexes are specific for cytotoxic cytokine releasing CD8+ T cells, that can target and lysate the tumor cells. In order to shortlist the best neoantigen candidates, the neoantigens etracted can be rank according to their binding affinity to MHC-I complexes. This was performed using netMHC (v 4.0)(Andreatta and Nielsen 511-517). The FASTA file containing the candidate neoantigens was input and 8 to 14mer peptide len gth was selected to obtain the widest library of neoantigens. The species/loci chosen was the standear HLA supertype representative. To test the system only alleles of class HLA-A were selected, including HLA-A0101, HLA-A0201, HLA-A0301, HLA-A2402, HLA-A2601. The options to sort by predicted affinity and save output in XLS format, were also selected in order to be able to easily postprocess the neoantigens. The results of the analysis can be found in the Ginga Github.

Additionally, to netMHC Ginga also includes the support to mixMHCpred, a software that adds mass-spectrometry analysis of the MHC-I alleles to improve the prediction of MHC-I-peptide binding (Bassani-Sternberg M et al.). However, this is currently under development.

Rank Binding Affinity - neoSearch

Although, netMHC can sort the peptide sequences by their ranking affinity this data is difficult to organis. Furthermore, sorting the peptides according to their binding affinity and not the parent sequence that they originate from is currently not possible. neoSearch sorts the XLS file from netMHC and returns the ranked neoantigens in a CSV file that can be easily read and processed to create high throughput therapeutics. The final expected results of the pipeline can be found in the Ginga Github and a brief summary of the results can be found in Table 9.

neoSearch has also the option to obtain the original exome mutated sequence and index position, of the mutation that coded the specific neoantigen. neoSearch uses the protein index and the VCF file outputted by Mutect2 to find the original coding sequence for the neoantigen candidate. The length of the sequence extracted can be adjusted accordingly, to detect the . The purpose of this trace back is to

Table 9 - First 19 neoantigens discovered for sample ERR2683866 for the HLA-A0101 allele. The SB label indicates the binding affinity (SB).