# INTRODUCTION

High throughput sequencing technologies are standard techniques in biomedical research and personalized medicine. Identification of genetic variations and their association with characterized phenotypes facilitates assessment of the risk of genetic-driven disorders and development of treatments and cures [PMID: 19474294; PMID: 19584936; PMID: 18460330; PMID: 22908194; PMID: 21199653].

The most commonly studied type of genomic sequence variation are single nucleotide polymorphism (SNP)[PMID: 14685227; A global reference for human genetic variation Nature 526, 68–74 (01 October 2015)].
Mutations associated with cancer are often found to be at low frequency (defined in the literature as low allele frequency) [PMID: 29123093; PMID: 29178133; PMID: 23415222; PMID: 21160474]. Previous work has reported mutations, found at low frequencies, that where the results of systematic errors during DNA storage and library preparation, such as "8-oxo-G" and "cytosine deamination", rather than true variants [PMID: 23303777; PMID: 28209900].

Sequencing depth is key to differentiate variants from random errors. In fact, collecting enough reads (to call with enough confidence the majority of variants) is still the bigest chanllenge in variant calling (I made this up from intuition. Is this true?). In order to partilly paliate the lack of enough reads, a local alignment can prevent from discarding up to 50% (I should check this and give exact numbers from our data?) of the reads. However, local alignments have their own problems, as a read could be misplaced, without being exposed by the map-quality assesment provided with the aligner software. On the other side, experimental procedures used to preserve a sample, such as formaldehyde-formalin parafin embedding (FFPE), will introduce errors that tend to accumulate at the ends. In such samples, many reads would not be mapped without soft (or hard) -clipping. Nevertheless, such reads would not show signal from the ends. In order to include such bases, tasmanian evaluates the softcliped region and includes it in the calculations if it has at least X% identity with the reference genome, where X could be adjusted to the analysis. 

A key variable that influences our hability to distinguish variants from random errors is the sequencing depth. As 

Chen and co-authors showed that damage can not be neglected as it leads to wrong variant calling at low allele frequencies, which are often associated with disease.

Due to room-temperature stability and preservation of cellular morphology[PMID: 24845430] cancer tumor samples available for research are often formalin-fixed-parafin-embedded (FFPE) prior to storage. Unfortunately, FFPE treatment is associated with DNA damage and that can be confused with low-frequency  mutations. Ancient DNA [PMID: 20028723] or ultrasonic shearing of the DNA [PMID: 23303777] have been also associated with systematic base changes. At low enough allele frequency, it is difficult to distinguish bases damaged during sample preservation from true somatic variants.

Cytosine deamination is normally a signature of FFPE treated samples[PMID: 31477010]. A previous work found no evidence of immediate deamination of DNA upon FFPE treatment (allele frequencies <2%) compared to the fresh frozen matched tissue [PMID: 29698444]. However, the mentioned study is restricted to the FFPE process itself and did not evaluate the impact of long-term storage of such samples. In practice, we observe  artifacts associated with FFPE treatment and storage.

Experimental protocols to assess the level of damage of DNA have been developed [PMID: 16333309]. PCR-based methods are limited to assessing presence of DNA fragmentation, but cannot measure base-level damage.

Bioinformatic assessment of the damage of a sequenced DNA sample can be used to compute an analytical threshold for acceptance or rejection of low allele fraction variants, leading to higher sensitivity and specificity in associations with phenotype, making a diagnosis or planning a treatment.

In a pioneering study, Costello and co-workers developed a Phred-like score, termed ArtQ, to evaluate the difference on counts of a mutation in its forward or reverse complement forms. Such score, though useful to classify the DNA sample in a particular quality bin, might not be extensible to influence downstream analysis, such as variant calling.

Another similar score, developed by Chen and co-workers, computes the degree of imbalance between variants detected in read1 and read2 in paired-read sequencing [Chen at al.]. This score, termed Global Imbalance Value (GIV) is an alternative to ArtQ and is useful within the same restrictions as the ArtQ.

Tools available for detection of artifacts in ngs data do not incorporate positional bias along the read [Costello et al.], discard reads intersecting repetitive DNA [Chen et al.] or have not been optimized for speed of execution or disk space.

Here, we report a new tool for the analysis of artifacts in high-througpht sequence data and we outline empirically derived criteria to include such analysis into different pipelines.


Moreover, in the science study, repetitive regions were not included in the analysis, leaving aside up to 80% of high quality (Qread>XX and Qmap>XX) data in FFPE samples.

In order to provide a robust and efficient tool to analyze and assess damage from DNA alignemnt data we developed tasmanian... We keep readability as a major component in this tool and any bioinformatician can alter the code without having to spend excesive ammounts of time.



# RESULTS:

How Tasmanian compares to other tools in the analysis of NGS missmatches? Picard collectSequencingArtifactMetrics [pmid:...] does not consider the positional information of the missmatches. It provides summarized tables containing error rates and phred-scale quality scores for each type of missmatches. Damage estimator [pmid...] wraps samtools mpileup [pmid:...] and provides detailed information about the position of the missmatch within each of the reads in paired end sequencing data. Nevertheless, second strand is not included in the analysis, leaving aside half of the available data. MPileup is a very powerfull tool for producing binnary VCF files containing all genomic positions. This information is produced at the expense of very intensive computation but is necessary to call variants. However, such information includes irrelevant details for assessment of randomly distributed, systematic errors. Hence, the time consumed by mpileup alone is not justified in this case.

Tasmanian tool compares the aligned data to the reference genome and counts each of the 12 substitution and 4 matches at every position along each of the reads in a paired end setting. As the data should be randomly distributed, the frequency of substitutions/matches should be proportional to the counts. Nevertheless, different experiments could lead to different systematic biases and the user can also compute the frequencies, normalized by the reference base at the indicated position.
Often times, the analysis of mismatches is restricted to genomic positions such as non-repetitive regions. Other times different regions could contain different types of missmatches, depending on sample preparation (e.g. target capturing different regions in different samples). While tools like bedtools [pmid:...] can split reads provided a file containing the regions of interest, such information is not sufficient to identify the position of a missmatch in the read. To solve that issue, tasmanian provides am optimized routine to include such information in the analysis, provided a file with the genomic positions of interests.
When a read overlaps with a defined genomic location, tasmanian can split the read and mask the intersecting and the complementary bases into different subreads, keeping the information about the position of such bases in the read. This procedure is done "on the fly" with an insignificant addition of time and computation power.
The user can also execute the analysis in a debugging mode, which will allow to serialize two complementary masked bam files.
After comparing our tool to available and well established alternatives, we found unexpected behaviours in simulated data. to assess such "errors" or ... we provide a suite with minimal sets of data to manually and minimalistically simulate artifacts to assess the results provided by the software and to ensure the absence of errors or even report otherwise.

We applied tasmanian to the analysis of systematic positional artifacts in samples differing in the solution used during acustic fragmentation of the DNA (water vs TE buffer?). Our results are consistent with previously published results [pmid...science Chen], with more G to T counts in read 1 for water sheared samples and more C to A counts in read 2. This data do not show a positional bias.

We also applied tasmanian to FFPE data and compared samples originated from protocols including or not an enzymatic repair step. Sorprisingly, we found that while repair has a profound impact in regions with unique sequences, it did not show any effect within repetitive regions. While the mechanism of such difference is beyond the scope of this study, such observation would not have been possible without dissecting the regions of the reads, keeping their position unaffected.

Tasmanian tool can read sam files. However, the recommended use is to pipe the text from the standard input. This way, files are stored as binaries and read with samtools or sambamba and filters and/or subsampling can be executed on the fly and piped into tasmanian.

The output consists in tables with counts for all detected observed_reference "observations" and each row corresponds to the base position within the read, for both read1 and read2.



# DISCUSSION:

Analysis of missmatches usually preceedes calling variants. It can be used as a first screen to decide if the sample is suitable to further variant calling analysis or it can be used as a cryteria to assign a qualitative measure to potential variants associated with desease.
While the extent of systematic missmatches due to artifacts associated to DNA manipulation have very low frequency, it has a big impact in variant calling where frequencies are as low or even lower than potential systematic evaluated artifacts. Hence, half of the data discarded by damage estimator will have a profound impact in further assesment of variant calling results.  <!-- ANALYZE DATA WITH  DIFFERENT AMOUNTS OF READS AND ASSESS THE IMPACT IN THE RESULTS -->



# Datasets:
For the TOPMed Phase 3 and CCDG projects, library core technicians used 600ng of DNA starting material with a minimum starting amount of 450ng. The library construction technicians fragmented the DNA with the Covaris LE220 Focused Ultrasonicator. A mean fragment size of ~375 bp was achieved. The library core technicians then prepared KAPA Hyper PCR-free libraries (Roche) using Perkin Elmer SciClone NGS (96-well configuration). The dual-same index system from Illumina was utilized to ensure DNA molecule integrity and data quality, which allows identification and elimination of chimeric molecules formed by cross-sample "index hopping". The library construction technicians assessed the libraries for quality and quantity using the HT DNA Hi Sens Dual Protocol Assay with the DNA Extended Range LabChip (GX, GX Touch HT) on the Perkin Elmer LabChip GX instrument using the manufacturer’s instructions. The library construction team passed undiluted libraries to the loading team for qPCR, pooling, and sequencing.
