CUSHAW Suite: Parallel and Efficient Algorithms for NGS Read Alignment

After the reference genomes of many organisms are sequenced in this post-genetic era, it has become an extremely important issue that how to do the re-sequencing and assembly for individual genomes from very large amount of reads. In this paper, we will present a re-sequencing tool designed for the Next Generation Sequencing (NGS) data. And these data are composed of a huge amount of short reads which will be aligned onto a reference genome. We modified and implemented the algorithm of Burrows-Wheeler Transform and FM-index to build the genome index of human, and proposed an idea to segment each short read into multiple non-overlapping seeds, which let us align short reads with large Hamming distance. Finally, we used the simulated datasets and real datasets from 1000 Genome Project to demonstrate the performance of our tool on a personal computer, and compared the results with widely used tools, bowtie and SOAPv2.

With the development of new sequencing techniques, whole genomes of many species have become available. This huge amount of data gives rise to new opportunities and challenges. These new sequences provide valuable information on relationships among species, e.g. genome recombination and conservation. One of the principal ways to investigate such information is multiple sequence alignment (MSA). Currently, there is large amount of MSA data on the internet, such as the UCSC genome database, but how to effectively use this information to solve classical and new problems is still an area lacking of exploration. In this thesis, we explored how to use this information in four problems, i.e. sequence orthology search problem, multiple alignment improvement problem, short read mapping problem, and genome rearrangement inference problem. For the first problem, we developed a EM algorithm to iteratively align a query with a multiple alignment database with the information from a phylogeny relating the query species and the species in the multiple alignment. We also infer the query's location in the phylogeny. We showed that by doing alignment and phylogeny inference together, we can improve the accuracies for both problems. For the second problem, we developed an optimization algorithm to iteratively refine the multiple alignment quality. Experiment results showed our algorithm is very stable in term of resulting alignments. The results showed that our method is more accurate than existing methods, i.e. Mafft, Clustal-O, and Mavid, on test data from three sets of species from the UCSC genome database. For the third problem, we developed a model, PhyMap, to align a read to a multiple alignment allowing mismatches and indels. PhyMap computes local alignments of a query sequence against a fixed multiple-genome alignment of closely related species. PhyMap uses a known phylogenetic tree on the species in the multiple alignment to improve the quality of its computed alignments while also estimating the placement of the query on this tree. Both theoretical computation and experiment results show that our model can differentiate between orthologous and paralogous alignments better than other popular short read mapping tools (BWA, BOWTIE and BLAST). For the fourth problem, we gave a simple genome recombination model which can express insertions, deletions, inversions, translocations and inverted translocations on aligned genome segments. We also developed an MCMC algorithm to infer the order of the query segments. We proved that using any Euclidian metrics to measure distance between two sequence orders in the tree optimization goal function will lead to a degenerated solution where the inferred order will be the order of one of the leaf nodes. We also gave a graph-based formulation of the problem which can represent the probability distribution of the order of the query sequences.

A number of next-generation sequencing (NGS) technologies such as Roche/454, Illumina and AB SOLiD have recently become available. These technologies are capable of generating hundreds of thousands or tens of millions of short DNA sequence reads at a relatively low cost. These NGS technologies, now referred as second-generation sequencing (SGS) technologies, are being utilized for de novo sequencing, genome re-sequencing, and whole genome and transcriptome analysis. Now, new generation of sequencers, based on the 'next-next' or third-generation sequencing (TGS) technologies like the Single-Molecule Real-Time (SMRT™) Sequencer, Heliscope™ Single Molecule Sequencer, and the Ion Personal Genome Machine™ are becoming available that are capable of generating longer sequence reads in a shorter time and at even lower costs per instrument run. Ever declining sequencing costs and increased data output and sample throughput for NGS and TGS sequencing technologies enable the plant genomics and breeding community to undertake genotyping-by-sequencing (GBS). Data analysis, storage and management of large-scale second or TGS projects, however, are essential. This article provides an overview of different sequencing technologies with an emphasis on forthcoming TGS technologies and bioinformatics tools required for the latest evolution of DNA sequencing platforms.

Towards standardization of the description and publication of next‐generation sequencing datasets of fungal communities

BackgroundNext generation sequencing (NGS) is a state of the art technology for microRNA (miRNA) analysis. The quantitative interpretation of the primary output of NGS i.e. the read counts for a miRNA sequence that can vary by several orders of magnitude (1 to 107) remains incompletely understood.FindingsNGS (SOLiD 3 technology) was performed on biopsies from 6 Barrett’s esophagus (BE) and 5 Gastroesophageal Reflux Disease (GERD) patients. Read sequences were aligned to miRBase 18.0. Differential expression analysis was adjusted for false discovery rate of 5%. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed for 36 miRNA in a validation cohort of 47 patients (27 BE and 20 GERD). Correlation coefficients, accuracy, precision and recall of NGS compared to qRT-PCR were calculated. Increase in NGS reads was associated with progressively lower Cq values, p < 0.05. Although absolute quantification between NGS reads and Cq values correlated modestly: -0.38, p = 0.01 for BE and -0.32, p = 0.05 for GERD, relative quantification (fold changes) of miRNA expression between BE &GERD by NGS correlated highly with qRT-PCR 0.86, p = 2.45E-11. Fold change correlations were unaffected when different thresholds of NGS read counts were compared (>1000 vs. <1000, >500 vs. <500 and >100 vs. <100). The accuracy, precision and recall of NGS to label a miRNA as differentially expressed were 0.71, 0.88 and 0.74 respectively.ConclusionAbsolute NGS reads correlated modestly with qRT-PCR but fold changes correlated highly. NGS is robust at relative but not absolute quantification of miRNA levels and accurate for high-throughput identification of differentially expressed miRNA.

Motivation and Objectives The output of a high throughput next generation sequencing (NGS) machine is a collection of short reads, which have to be properly assembled in order to reconstruct the original DNA sequence of the analyzed organism (Metzker, 2010; Earl, 2011). The DNA sequence assembly process is based on aligning and merging these reads for effectively reconstructing the real primary structure of the DNA sample sequence or reference genome. The use of NGS machines results in much larger sets of reads to be assembled, posing new problems for computer scientists and bioinformaticians. In particular, a relevant issue is related with the trade-off between precision of the assembly process and its computational time, stating the need for faster methods that can keep pace with the speed and volume of reads that are generated with NGS. An important step in DNA assembly is the identification of a subset of read pairs that have a high probability of being aligned sequentially in the reconstruction. Such step is often referred to as filtering, and amounts in selecting a significantly smaller subset of the initial set of read pairs (whose dimension is quadratic in the number of initial reads) that can be then processed by an alignment algorithm, usually quite time consuming. The desired effect of filtering is then to quickly filter out from the candidate set of read pairs those that would not provide a good alignment in the following phase. The computation cost of filtering should then be balanced by the speed-up obtained when a smaller set of read pairs is considered for alignment. In this work we propose and test the use of alignment free distances to evaluate the similarity between two short reads as a technique for filtering good read pairs to be assembled. The method operates in constant time in the string length and is tested in its ability to emulate, with a proper level of precision, much more time consuming methods to evaluate the similarity between short DNA sequences, such as the established Needleman-Wunsch edit distance (Needleman, Wunsch, and Christian, 1970), often used in the final step of the assembly procedure. These preliminary experiments show the efficacy of this approach for filtering the promising read pairs - eligible candidates to successfully assemble the entire genome of a given organism. Therefore, the alignment free reads filtering may significantly accelerate the assembly process without a substantial loss in accuracy of the DNA sample sequence reconstruction. Methods ODNA sequence assembly The DNA sequence assembly process is based on the alignment and merging of reads (stretch of sequences) in order to reconstruct the original primary structure of the DNA sample sequences. Given a set of sequences S = {s1, s2, …, sn}, where s ∈ S is a fragment of the primary structure of DNA (read) (e.g. s = {ATTCGA...CTGACT}), assembly is in charge of building the longest sequence from the set S where each pair of consequent reads obey certain similarity conditions. DNA read pairs filtering and Alignment Free Distance This step identifies the promising read pairs in order to reduce the amount of input data given to the real assembly algorithm. We adopted a very quick measure of the similarity between two reads, Alignment Free (AF) based distance (Vinga and Almeida, 2003). AF computes the similarity of two strings based only on the dictionary of their substrings, irrespective of their relative position. As a dictionary we considered the set of 4-mer (sequences composed of 4 different nucleotides) and then built a profile for each read composed by the relative frequencies of each 4-mer in the read. The Euclidean distance between the profiles of two reads was taken as an inverse measure of the similarity of the two reads and thus as an indication that the two reads formed a promising pair to be considered in the assembly phase. AF filtering was then used defining a proper threshold on the AF distance and discarding all the pairs that exhibited a AF distance above the threshold. Computational complexity of AF distance is a constant linearly bounded by the number of k-mers adopted and the length of the strings to be compared. Comparing with other distances: Needleman-Wunsch and “Bowtie” distance Along with AF we considered the well-established Needleman-Wunsch edit distance (NW) and compared them in their ability to identify significant pairs. This comparison was based on the computation of a sort of perfect distance computed after an alignment over an already known sequence has been performed. Such distance, referred to as Bowtie distance (BT), was obtained as follows: a. a large number of reads coming from a known sequence were considered; b. hese reads were aligned over the known sequence using the standard Bowtie algorithm (Langmead et al, 2009); c. any two reads received a maximum BT distance if their alignment did not intersect over the reference sequence, else they received a distance inversely proportional to their intersection over the sequence (e.g., they would have BT distance equal to 0 if they were aligned one on top (or inside) of the other by the Bowtie algorithm). By construction we assumed BT distance to be the reference distance, e.g., the distance that expressed the best possible alignments - being based on the knowledge of the reference sequence - and tested the correlation of AF and NW with BT; moreover, we verified the ability of AF and NW to predict that a given read pair had BT distance above or below a given threshold. Results and Discussion For our test we considered the E.Coli genome and a set of reads from this genome reads obtained by Roche 454 sequencing machine. Reads have average length of ~235 nucleotides and standard deviation of approx. 10 (the large majority of them having length in the interval 225-245). Reads were aligned with the reference sequence with Bowtie and then 100,000 were sampled at random according to their alignment along the sequence. Reads were considered both forward and reversed, giving rise to a total of 200,0002 read pairs. All 620,798 read pairs with BT distance &lt; 1 were considered for the experiments; then, out of the remaining pairs, 233,099 were sampled at random. A total of 853,897 read pairs composed the working data set. For all these reads, NW distance and AF distance over the 4-mer were computed. AF, NW and BT distances were all normalized between 0 (maximal similarity) and 1 (maximal dissimilarity). The first interesting results was that the correlation between distances showed that AF approximates BT somehow better than NW: we obtained a correlation coefficient of 0.761 for AF and BT, compared with a smaller 0.706 when NW and BT were considered (coherently, correlation between AF and NW is 0,721). The second interesting results was obtained when we compared the ability of AF and NW to predict whether BT was above or below a given threshold. We defined a threshold predictor for a given function F2 based on function F1 and on a given pair a1, a2 as follows: if (F1 &lt; a1) then predict (F2 &lt; a2), else predict (F2 ≥ a2). To a given pair (a1, a2), we associated the measure of True Positive rate (TP) (percentage of cases where (F1 &lt; a1) and (F2 &lt; a2) and of True Negative rate (TN) (percentage of cases where (F1 ≥ a1) and (F2 ≥ a2); analogously we defined False Positive rate (FP) and False Negative rate (FN). For each (a1, a2) with both values ranging from 0 to 1, we then computed, with step 0.05, the positive and negative error rates taking AF as a predictor of BT and NW as a predictor of BT. Part of the results are summarized in the charts of Figure 1, that show for 3 different levels of a2 (0.1, 0.2, and 0.3) the precision of the predictors (y-axis) when the value of a1 is changed (x-axis), both when AF is used as a predictor of BT (charts in the first row) and when NW is used as a predictor of BT (charts in second row). Similar results are obtained also for other levels of a2, here omitted for brevity. The curves bring to light very clearly how AF is a very good threshold predictor for BT for the considered data; despite its light computational complexity, it appears to perform significantly better than the more complex NW edit distance when its ability to support a threshold predictor is considered. Acknowledgements The authors are partially supported by the FLAGSHIP ”InterOmics” project (PB.P05) funded by the Italian MIUR and CNR institutions, and by the cooperative programme 2010–2012 between the National Research Council of Italy (CNR) and the Polish Academy of Sciences (PAN). References Earl D et al (2011); Assemblathon 1: A competitive assessment of de novo short read assembly methods; Genome Research, 21 Langmead B, Trapnell C, et al (2009); Ultrafast and memory-efficient alignment of short DNA sequences to the human genome; Genome Biology, 10:R25 Metzker ML (2010); Sequencing technologies — the next generation; Nat Rev Genet., 11(1) Needleman SB, Wunsch CD. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology 48(3): 443–53 Vinga S, Almeida J (2003); Alignment-free sequence comparison—a review, Bioinformatics 19 (4), 513-523 Note: Figures and tables are available in PDF version only.

Background: The structural annotation of genes serves as an important resource in the genomics era. Next-generation sequencing (NGS)-based RNA-seq has been used to reconstruct transcriptomes by assembling short reads; however, coupled alternative RNA processing events occurring over long distances cannot be captured. Although the more recently developed third-generation sequencing (TGS) technologies can generate longer reads, challenges exist in full-length transcriptome (FLT) reconstruction due to their higher sequencing error rate and relatively lower throughput. Results: Here, to address the above-stated technical challenges and provide practical FLT annotation in rat, we developed a hybrid-sequencing workflow for reconstructing FLT and applied it to rat hippocampus. Our workflow combines high sequencing accuracy from NGS and long read length from TGS, and has two main features: (1) cDNA normalization and size fractionation to facilitate the profiling of diverse transcripts at an affordable throughput; and (2) a computational approach to correct TGS sequencing errors using NGS reads. The error correction substantially improved the read accuracy with respect to gene expression levels. To determine high-confidence FLT, we used independent 5´ and 3´-end profiling data to validate and filter for a final FLT set in rat hippocampus, which contains 28,268 transcripts derived from 6,380 RefSeq gene loci and 849 unannotated loci. Comparing our FLT to RefSeq annotation, we found that our FLT displayed a higher diversity and more than 90% of our FLT isoforms were not annotated in RefSeq. Moreover, we discovered significant co-occurrence of alternative RNA processing events between transcript 5´-ends and 3´-ends, alternative 5´ and 3´ splice sites, and multiple retained intron. Integrating with polysome profiling and ribosome footprinting data, we predicted active or inactive translational status of transcript isoforms, and reconstructed an open reading frame (ORF)-eome in rat hippocampus, which was largely supported by mass spectrometry-based proteomics datasets. Compared to RefSeq-annotated ORFs, only less than one-third of our ORFome were annotated in RefSeq, whereas the majority were ORF variants, and about one-tenth were novel ORFs. Finally, we identify isoforms with specific subcellular localization patterns in neurons. Conclusions: Altogether, our data expand the current rat gene annotation and provide a rich resource for future functional studies. The developed hybrid-sequencing workflow and downstream analysis can be easily adapted to uncover the complex transcriptome/ORFeome landscapes expressed in other tissues from other organisms. Full article: PMID: 31676752 PMCID: PMC6825209 DOI: 10.1038/s41467-019-13037-0

Collection and storage of HLA NGS genotyping data for the 17th International HLA and Immunogenetics Workshop

BackgroundWith the advent of Next-Generation Sequencing technologies (NGS), a large amount of short read data has been generated. If a reference genome is not available, the assembly of a template sequence is usually challenging because of repeats and the short length of reads. When NGS reads cannot be mapped onto a reference genome alignment-based methods are not applicable. However it is still possible to study the evolutionary relationship of unassembled genomes based on NGS data.ResultsWe present a parameter-free alignment-free method, called , based on variable-length patterns, for the direct comparison of sets of NGS reads. We define a similarity measure using variable-length patterns, as well as reverses and reverse-complements, along with their statistical and syntactical properties. We evaluate several alignment-free statistics on the comparison of NGS reads coming from simulated and real genomes. In almost all simulations our method outperforms all other statistics. The performance gain becomes more evident when real genomes are used.ConclusionThe new alignment-free statistic is highly successful in discriminating related genomes based on NGS reads data. In almost all experiments, it outperforms traditional alignment-free statistics that are based on fixed length patterns.

In medical and genomic research, the Next Generation Sequencing (NGS) has a major role. Presently, NGS data are produced at the rate of 10 TB a day and challenge the storage and data processing capacities. These huge datasets are being used by a wide sort of applications such as customized cancer healing and precision medicine. NGS technologies offer prospects for understanding unidentified species and complex syndrome. To utilize genomic data for such applications, the genomic data in the form of short reads produced by NGS initially has to be assembled into whole genome sequence. And then, the sequences have to be compared for similarity and variation discovery which will be useful for analyzing and arriving at health-related solutions. In this paper, the fundamental methods for short read alignment such as assembly-based and alignment-based methods are discussed. Followed by which, the different ways to compare the sequences to check the alignment for similarity/dissimilarity discovery are discussed. This comparative analysis report can be utilized for health-related medical decisions.

### Widespread RNA and DNA Sequence Differences in the Human Transcriptome Li et al Science . 2011;333:53–58. New genomic studies show that the transfer of genetic information from DNA to RNA within humans is not very faithful. If true, frequent, or purposeful, then the implications for human health and disease are significant, but the burden of proof required is significant. Commentary: Molecular biology and genetics has a “central dogma.” First postulated by Francis Crick1,2 in 1958, shortly after the discovery of the double helix, it states that the molecular transfer or flow of genetic information occurs only as DNA→RNA→ protein. This argument is intrinsic to our understanding of genotypes, phenotypes, and human disease. The central “dogma,” a word whose precise meaning Crick acknowledged he did not precisely know3 and that could have been given the less glamorous moniker “hypothesis,” is frequently misunderstood as also implying high accuracy or even inflexibility.2 We now know that genetic information flow can be sometimes corrupted and can also occur in other directions in particular biological systems. A recent study by Li and colleagues,4 from the laboratory of Vivian Cheung at the University of Pennsylvania, now demonstrates apparent widespread DNA→RNA sequence differences in humans. The experiments they conducted involved sequencing reverse-transcribed RNA from cultured B cells isolated from 27 different humans and comparing them with their host genome sequences. Surprisingly, after appropriate quality controls, they discovered numerous base positions at which the RNA sequence was different from the DNA sequence in multiple molecules, a phenomenon they refer to as RDD (RNA-DNA difference). Specifically, they identified 28766 RDDs in the 27 individuals at 10210 exonic sites within 4741 genes. Each individual had high-quality RNA-based sequence data at 9.8 million sites: so, RDD was detected at a rate of 1.1×10−4 per exonic site (base) across human genes. Interestingly, they observed all 12 types …

The advances of next-generation sequencing technology have facilitated metagenomics research that attempts to determine directly the whole collection of genetic material within an environmental sample (i.e. the metagenome). Identification of genes directly from short reads has become an important yet challenging problem in annotating metagenomes, since the assembly of metagenomes is often not available. Gene predictors developed for whole genomes (e.g. Glimmer) and recently developed for metagenomic sequences (e.g. MetaGene) show a significant decrease in performance as the sequencing error rates increase, or as reads get shorter. We have developed a novel gene prediction method FragGeneScan, which combines sequencing error models and codon usages in a hidden Markov model to improve the prediction of protein-coding region in short reads. The performance of FragGeneScan was comparable to Glimmer and MetaGene for complete genomes. But for short reads, FragGeneScan consistently outperformed MetaGene (accuracy improved ∼62% for reads of 400 bases with 1% sequencing errors, and ∼18% for short reads of 100 bases that are error free). When applied to metagenomes, FragGeneScan recovered substantially more genes than MetaGene predicted (>90% of the genes identified by homology search), and many novel genes with no homologs in current protein sequence database.

Contigs assembled from the second generation sequencing short reads may contain misassemblies, and thus complicate downstream analysis or even lead to incorrect analysis results. Fortunately, with more and more sequenced species available, it becomes possible to use the reference genome of a closely related species to detect misassemblies. In addition, long reads of the third generation sequencing technology have been more and more widely used, and can also help detect misassemblies. Here, we introduce ReMILO, a reference assisted misassembly detection algorithm that uses both short reads and PacBio SMRT long reads. ReMILO aligns the initial short reads to both the contigs and reference genome, and then constructs a novel data structure called red-black multipositional de Bruijn graph to detect misassemblies. In addition, ReMILO also aligns the contigs to long reads and find their differences from the long reads to detect more misassemblies. In our performance test on short read assemblies of human chromosome 14 data, ReMILO can detect 41.8-77.9% extensive misassemblies and 33.6-54.5% local misassemblies. On hybrid short and long read assemblies of S.pastorianus data, ReMILO can also detect 60.6-70.9% extensive misassemblies and 28.6-54.0% local misassemblies. The ReMILO software can be downloaded for free under Artistic License 2.0 from this site: https://github.com/songc001/remilo. baoe@bjtu.edu.cn. Supplementary data are available at Bioinformatics online.

Sensitive detection of internal tandem duplication (ITD) mutation of FLT3 is very important in acute myeloid leukemia. To increase detection sensitivity of FLT3-ITD, we developed new detection algorithm using next generation sequencing (NGS) data. We validated results using nested polymerase chain reaction (PCR) methods. We compared results of NGS data, nested PCR and conventional PCR methods. First, using whole exome sequencing data of 83 AML patients, we applied calling algorithm for FLT3-ITD. Briefly, to detect ITDs with NGS data, the reads are aligned to a reference sequence (UCSC hg19), with BWA which is a read aligner allowing soft-clipping. Some reads can be an indication of the occurrence of ITD and BWA aligns those reads as soft-clipped. Second, we deigned two types of primer for Nested PCR. The first primer was targeted wildly for between exon14 and exon15 of FLT3 gene. Nested PCR primer was deigned to target previously reported regions which are frequently occurred ITD mutation. PCR reactions of two steps were performed using the PCR primers sequentially. In these 83 patients, FLT3-ITD was positive only in 7 patients when tested by conventional PCR methods. When NGS detection method was applied, this resulted in positive FLT3-ITD in 11 patients (11/83, 13%). When validation was performed using nested PCR, FLT3-ITD was confirmed in all of 11 patients. Nested PCR detected additional 4 patient positive for FLT3-ITD in this population. For 68 patients, FLT3-ITD was negative by both NGS and nested PCR method. Overall, NGS method improved sensitivity of FLT3-ITD detection by 57% in this population. And the concordance rate of NGS method and nested PCR was 95.2% (79/83). Then we investigated clinical significance of sensitive FLT3-ITD detection. For this, we performed nested PCR and conventional PCR at the same time in 238 AML patients to detect FLT3-ITD. Positive rate for FLT3-ITD was 20% (48/238) and 10% (24/238) by nested PCR and conventional PCR respectively. When survival analysis was performed, among patients with negative FLT3-ITD result by conventional PCR, patients who showed positive for FLT3-ITD by nested PCR had shorter overall survival compared to those who showed negative for FLT3-ITD by nested PCR. (p = 0.03). This implies that sensitive FLT3-ITD detection using nested PCR is clinically meaningful. Diagnosis of FLT3-ITD is very important genetic factor, leading a therapeutic direction for AML patient. Here we report that we have developed alternative more sensitive detection methods for FLT3-ITD based on nested PCR and NGS. Sensitive detection of FLT3-ITD was clinically meaningful, suggesting that these methods should be incorporated in a future clinical practice. Also, we want to note that, NGS method is capable of quantifying FLT3-ITD size and amount in AML patients. Citation Format: Daeyoon Kim, Yoojin Hong, Youngil Koh, Sung-Soo Yoon, Choong-Hyun Sun, Kwang-Sung Ahn, Seungmook Lee, Hongseok Yun, Suyeon Lee. Improved sensitive detection method of FLT3 (FMS-like tyrosine kinase) internal tandem duplication (ITD) mutation using next-generation sequencing technology and nested PCR. [abstract]. In: Proceedings of the 107th Annual Meeting of the American Association for Cancer Research; 2016 Apr 16-20; New Orleans, LA. Philadelphia (PA): AACR; Cancer Res 2016;76(14 Suppl):Abstract nr 3186.

A Recurrent T-Cell Receptor γ Rearrangement of the Vγ9-JγP Type Represents the Major Clone of γδ T-Cells in Normal Individuals and Is Highly Suppressed in Lymphoproliferative Disorders