Detrimental effects of duplicate reads and low complexity regions on RNA- and ChIP-seq data
BMC Bioinformatics volume 16, Article number: S10 (2015)
Adapter trimming and removal of duplicate reads are common practices in next-generation sequencing pipelines. Sequencing reads ambiguously mapped to repetitive and low complexity regions can also be problematic for accurate assessment of the biological signal, yet their impact on sequencing data has not received much attention. We investigate how trimming the adapters, removing duplicates, and filtering out reads overlapping low complexity regions influence the significance of biological signal in RNA- and ChIP-seq experiments.
We assessed the effect of data processing steps on the alignment statistics and the functional enrichment analysis results of RNA- and ChIP-seq data. We compared differentially processed RNA-seq data with matching microarray data on the same patient samples to determine whether changes in pre-processing improved correlation between the two. We have developed a simple tool to remove low complexity regions, RepeatSoaker, available at https://github.com/mdozmorov/RepeatSoaker, and tested its effect on the alignment statistics and the results of the enrichment analyses.
Both adapter trimming and duplicate removal moderately improved the strength of biological signals in RNA-seq and ChIP-seq data. Aggressive filtering of reads overlapping with low complexity regions, as defined by RepeatMasker, further improved the strength of biological signals, and the correlation between RNA-seq and microarray gene expression data.
Adapter trimming and duplicates removal, coupled with filtering out reads overlapping low complexity regions, is shown to increase the quality and reliability of detecting biological signals in RNA-seq and ChIP-seq data.
Next generation sequencing (NGS) technology is primarily based on massively parallel sequencing of millions of short reads from DNA/RNA samples, although the read length has been increasing [1, 2]. The costs of NGS have rapidly dropped  and, consequently, there has been a relatively rapid shift from the use of microarrays to RNA-seq data to study transcription. This increased reliance on NGS necessitates examination of analysis steps that may affect the quality of the data and the interpretations drawn from it.
Although different types of NGS experiments and library preparation protocols dictate downstream processing steps, removing sequence adapters used to construct the short read library  as well as removing low quality bases from short reads  are the typical quality control steps. This is followed by removing duplicate reads, which can arise during library amplification by polymerase chain reaction . The rationale behind this step is that such duplicate reads may lead to erroneous conclusions regarding the true level of biological signal, e.g., variant detection in DNA-seq data , gene expression in RNA-seq data , quantification of gene rearrangements , and in ChIP-seq data . Two schools of thought have emerged in the field regarding how to best address this issue. The first is a widely used practice to remove all duplicate or low complexity reads from the dataset, presuming they are a source of potential bias [7, 8, 10–12]. The second believes that these duplicates may be true unique observations and their removal introduces bias on its own . Although the effect of duplicate reads has been investigated in DNA sequencing , the question of how duplicates affect biological signals from gene expression in RNA-seq experiments and motif detection in ChIP-seq experiments remains open-ended.
The presence of low complexity [14, 15] and repetitive elements  in the reference genome received less attention in how they may affect the conclusions of biological studies. Such regions complicate alignment because short reads originating from them can be mapped to multiple locations, making their interpretation challenging [17, 18] or even impossible. This problem is not small, as eukaryotic genomes can be very rich in repeats; for example, some have estimated that the human genome contains ~47% repetitive regions . Although recent findings from the ENCODE project suggest that the genome is pervasively transcribed [20, 21], RNA molecules originating from low complexity and repetitive regions are a potential problem for analysis of RNA-seq data because they may promiscuously align throughout the genome . This also causes a problem for motif detection within protein-DNA interaction regions, identified by ChIP-seq [18, 22]. The use of paired-end sequencing, such as implemented by Illumina technologies, helps to control for proper mapping of reads by identifying the discordant read pairs whose mapped loci deviate from the expected orientation and insert size. The problem of multi-mapped reads has also been investigated in detection of structural variants, with the general consensus being to ignore them [23, 24].
In this work, we systematically investigated the effect of adapter trimming, duplicate reads removal, and filtering out reads overlapping low complexity regions upon biological signal in RNA-seq and ChIP-seq experiments. At each processing step, we used an indirect measurement of the strength of the biological signals by performing pathway- and gene ontology enrichment analyses of genes detected from RNA-seq data, and transcription factor binding sites enrichment analyses in peaks detected from ChIP-seq data. Our rationale here is that, if a processing step leads to a more significant enrichment p-value, that processing step is likely to positively influence biological signal.
Removal of reads overlapping low complexity regions (referred hereafter as "low complexity reads") received less attention. For example, this step has been included into a set of quality control steps in PRINSEQ tool . To allow isolated testing of the effect of low complexity reads removal, we introduce a simple post-alignment filtering tool, RepeatSoaker, that filters out reads overlapping with a user provided template file which contains genomic coordinates of low complexity regions. Designed to be aligner-independent, RepeatSoaker processes the aligned data in BAM format, removes low complexity reads, and outputs a cleaned BAM file and filtering statistics. RepeatSoaker is a straightforward method to remove alignment artifacts from NGS data, designed to eliminate potential false positive reads in quantifying transcript expression. Extendable to any other sequencing technology where low complexity reads may introduce bias, such as ChIP-seq, we envision RepeatSoaker could become a standard step helping to better structure reproducible NGS pipelines.
We applied a combination of adapter trimming, duplicate removal, and filtering out low complexity reads with RepeatSoaker to RNA-seq and ChIP-seq experiments, and investigated how each step affects the results of downstream enrichment analysis. Our results show that adapter trimming increases the significance of gene ontology and pathway enrichment analyses in RNA-seq data, and strengthens motif detection in ChIP-seq data. The duplicate removal step, despite decreasing the number of reads, further helps to increase the significance of biological signals, especially when coupled with adapter trimming. Filtering out low complexity reads with RepeatSoaker has minor effect on the total number of reads, yet this step had a positive effect on the detection of biological signals. Our study suggests that adapter trimming and duplicates removal are important steps in detecting stronger biological signals within RNA-seq and ChIP-seq data, and optional filtering out of reads overlapping low complexity regions will further increase the quality of conclusions.
The effect of adapter trimming, duplicate removal and filtering out reads overlapping low complexity region was investigated in 2 types of data: RNA-seq and ChIP-seq.
Human RNA-seq data were obtained from 60 Sjogren's syndrome patients and 30 healthy controls (unpublished data). In short, peripheral blood was collected in PaxGene tubes (BD Diagnostics; Franklin Lakes, NJ) and RNA was isolated using standard protocols (Qiagen, Inc., Valencia, CA). Globin transcripts were removed using GlobinClear (Life Technologies; Grand Island, NY) and samples were prepared for sequencing using the NuGen ENCORE complete kit (San Carlos, CA). The paired-end 100 bp-long sequencing was performed using the Illumina HiSeq 2000 employing standard procedures. For the 54 samples used in RNA-seq experiment, matching gene expression profiling was performed using Illumina HumanWG-6 v3.0 BeadChips. Pearson's correlation coefficient was used to test the agreement between human RNA expression estimates measured by RNA-Seq (log2 FPKM counts) and microarrays (log2 intensities) in the cohort of patients diagnosed with Sjogren's syndrome.
Mouse RNA-seq data were obtained from 3 wild type mice and 3 mice expressing the D842V mutant form of PDGFRα (platelet-derived growth factor receptor alpha), which was previously described . Cell suspensions from neonatal dermis were prepared from 3 day old mouse skin after separating dermis/epidermis and digesting the dermis with 0.35% collagenase type 1 for 60 minutes. Subsequently, Nestin-GFP+ singlets were sorted with a MoFlo XDP cell sorter (Beckman Coulter). RNA was isolated from 2-3 million cells using RNAeasy kit (Qiagen). cDNA libraries were prepared with NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs) according to the manufacturer's protocol. In short, mRNA was isolated from 1µg purified total RNA with oligo dT beads, and fragmented. Then, first and second strand cDNA were synthesized, followed by purification using Agencourt AMPure XP beads (Beckman Coulter). The second strand cDNAs were end-repaired, A-tailed, and adaptor-ligated. Size-selected DNA with Agencourt AMPure XP beads was enriched by 13-cycle PCR each with index primers, and again purified using the beads. The each indexed library was analyzed by Agilent 2200 TapeStation system (Agilent), and RNA integrity number equivalent (RINe) were ranged from 9.2 to 9.6. Then libraries were quantified with Qubit 2.0 Fluorometer (Thermo Fisher Scientific), and pooled for sequencing.
ChIP-seq data were obtained from 10 systemic lupus erythematosus patients and 10 healthy controls of European descent (unpublished data). Briefly, all nucleated cells were isolated from human blood using PolyPrep (Sigma, Deisenhofen, Germany) density gradient medium. Proteins were cross-linked to the DNA using formaldehyde and protein-DNA complexes were immunoprecipitated using a polyclonal rabbit anti-PU.1 (Spi-1) antibody (sc-352, Santa Cruz Biotechnology, Santa Cruz, California, USA). Individual sequencing libraries were prepared for each individual using the Illumina ChIP-Seq DNA Sample Prep Kit (Illumina, San Diego, California, USA). Sequencing was done using the Illumina HiSeq 2000 platform with 5 samples per lane. The case-control samples were sequenced on the same lane, e.g., 3 cases+2 controls in one lane.
Quality of raw sequence data was assessed using FASTQC. Adapter trimming was performed using Trimmomatic v0.30 program . The reads were cropped to the length of 70 bp, the adapters were trimmed, and bases having quality below 20 on Phred 33 scale were also cut. The reads with length less than 10 bases were discarded. Only the paired reads were used for subsequent analyses. Duplicates removal was performed using PICARD MarkDuplicates tool. The summary statistics were obtained using SAMTOOLS FLAGSTAT command .
To remove reads overlapping any of the regions identified by RepeatMasker [28, 29], we implemented RepeatSoaker (https://github.com/mdozmorov/RepeatSoaker). RepeatSoaker utilizes a user-provided list of genomic coordinates of low complexity regions in BED format . Note that the coordinates should correspond to the organism and genome assembly of the original BAM files. We provide a Makefile that automates generation of BED files containing genomic coordinates of all regions identified by RepeatMasker for GRCh37/hg19 and NCBI37/mm9 genomes.
RepeatSoaker provides flexibility to set a threshold for deciding whether a read should be kept or filtered due to its overlap with a low complexity region. A user can set the percent overlap threshold, e.g., 75%, 50%, 25%. A read overlapping with a low complexity region more than the threshold, e.g., more than 75%, is filtered. We tested the effect of removing reads overlapping with low complexity regions more than 75%, 50%, 25%, and 0%. A 0% threshold indicates that a read is filtered if it is immediately proximal to a repeat region.
For the human and mouse RNA-seq data, raw FASTQ files were aligned to the human and mouse genomes (hg19/mm9, respectively) using TOPHAT . The read counts per gene or transcript were generated using HTSEQ-COUNT. Differentially expressed (DE) transcripts were determined using DESeq R package with a false discovery rate (FDR) q-value of <0.05 and a fold change of >2 or <0.5. All data manipulations were performed in the R/Bioconductor environment .
For the human ChIP-seq data, raw FASTQ files were aligned to the human hg19 genome using bowtie2 . The PU.1 binding peaks were called using MACS2 , and the consensus motifs were detected using MEME-ChIP suite .
Systematic testing of data preprocessing steps
To elucidate the effects of adapter removal, elimination of duplicates, and filtering of low complexity regions, we performed systematic testing of sequencing data with and without applying these three preprocessing steps (Figure 1). At each step, we compared the alignment statistics, the number of differentially expressed genes (RNA-seq) or identified transcription factor binding peaks (ChIP-seq), and the results of Gene Ontology, KEGG and Reactome pathway enrichment analyses (RNA-seq) and motif enrichment analyses (ChIP-seq). We also compared combinations of data preprocessing steps, e.g., how duplicate removal affected trimmed and untrimmed data. At each comparison, we evaluated how the data preprocessing steps affected biological signals as judged by the functional enrichment analysis.
Adapter removal increases the data quality and the strength of biological signals
Before investigating the effects of low complexity region filtering, we assessed how adapter- and duplicate removal affected quality of alignment of the sequencing data. Adapter trimming increased the number of total reads in RNA-seq data, due to the fact that Trimmomatic had cut low quality bases from the middle portion of some of the reads, thus splitting some of them into multiple shorter reads, which still survived the minimum length (10 bp) threshold (Table 1 Additional Files 1 and 2). Such cuts resulted in more reads with mates mapped to a different chromosome (referred hereafter as "mismapped reads"). Adapter trimming of ChIP-seq data slightly decreased the total number of reads, as compared with unprocessed data (Additional File 3). However, the percentage of properly paired reads increased while the percent of singletons and mismapped reads decreased (Table 2). In summary, the adapter-trimming step altered the properties of sequencing data depending on the quality of the unprocessed data.
To identify how the adapter-trimming step affects biological signals in RNA-seq data, we evaluated the total number of differentially expressed genes before and after trimming, and their functional enrichment results. The total number of differentially expressed genes remained virtually unchanged, as well as the order of enrichment p-values for KEGG, GO, and Reactome pathway enrichment results, and Pearson's correlation coefficient with the microarray gene expression data (Table 1 Additional Files 4 and 5). In the case of ChIP-seq experiments, we evaluated the number and the significance of detected motifs within ChIP-seq binding peaks. As we have evaluated genome-wide binding of the PU.1 protein, the most significant motif enriched in the detected peaks was SPI1 (aka PU.1, Additional File 6). The total number and the significance of the detection of this motif increased following the adapter trimming step. These results support the notion that adapter trimming alone increases the strength of biological signals within it .
Removing duplicates negatively affects alignment statistics but strongly improves the detection of biologically relevant signal
Removing duplicates had the greatest effect on the alignment statistics, decreasing the total number of reads by ~40% in RNA-seq and ChIP-seq data, as compared with the unprocessed data (Additional Files 1, 2, 3). This step also decreased the percentage of singletons (Tables 1 and 2). Yet, this step had a slight negative effect on the percentage of properly paired reads at the expense of increased percentage of mismapped reads. Overall, the effect of duplicates removal on the alignment statistics was detrimental.
Despite the lower number of short reads remaining after duplicates removal, the number of differentially expressed genes increased (Table 1). Although the lists of differentially expressed genes before and after duplicate removal did not show complete overlap (Additional Files 4 and 5), the order and the significance of the enriched KEGG, GOs and Reactome pathways remained similar indicating that the biological signal was retained and became more significant, as the enrichment p-values were also increased (Table 1). We did not observe any increase in correlation of RNA-seq and microarray gene expression data.
In contrast to a larger number of differentially expressed genes when removing duplicates from the data, the number of detected motifs in ChIP-seq data decreased. However, the significance of the detected PU.1 motif increased (Table 2 Additional File 6). In summary, these results suggest that removal of duplicates, although decreasing the total number of reads, increases the significance of biological signals in both RNA-seq and ChIP-seq data.
Trimming adapters coupled with duplicates removal synergistically improves the quality of sequencing data
Adapter trimming coupled with duplicates removal led to an increase in the percentage of mismapped reads in RNA-seq data (Table 1). Consequently, the percentage of properly paired reads has decreased. However, the number of differentially expressed genes was larger than in the unprocessed data and their enrichment analyses results, except for the adapter-trimming step, were comparably significant.
Applied to ChIP-seq data, trimming the adapters coupled with removing duplicates improved overlap statistics, making them better than the unprocessed data (Table 2). Notably, despite lower total number of detected motifs, the detection of the PU.1 motif was the most significant (Additional File 6). Overall, our results suggest that both adapter trimming and duplicates removal steps help to emphasize biological signal in RNA-seq and ChIP-seq data.
Removing low complexity regions improves detection of true biological signal
We used different stringency thresholds for filtering reads overlapping low complexity regions. A threshold was defined by the percent of overlap of a read with low complexity regions. For example, 75% threshold indicates that a read will be filtered only if it overlaps with a low complexity region at least 75% of the read's total length. A special case of 0% threshold indicates that a read can be located side-by-side with a low complexity region to be considered for removal. Thus, lowering percent overlap threshold corresponds to a more stringent threshold for filtering out of reads overlapping low complexity regions, with 0% threshold indicating the most aggressive filtering.
Applying RepeatSoaker to RNA-seq and ChIP-seq data decreased the number of reads by ~3%, independently of the threshold used. However, the alignment statistics of RNA-seq and ChIP-seq data improved in all conditions, when compared with the corresponding data preprocessing steps. For example, comparing the alignment statistics for unprocessed and processed (filtered) data shows not only increase in the percentage of properly paired reads, but also in decrease in mismapped reads (Additional Files 1, 2, 3). Increasing the stringency of filtering of the reads overlapping low complexity regions further increased the alignment statistics. This suggests that simply removing reads overlapping low complexity regions is an essential step in improving the quality of RNA- and ChIP-seq data.
Reflecting the improvement in the alignment statistics, the significance of the GO/KEGG/Reactome enrichments (RNA-seq) and the PU.1 motif enrichments (ChIP-seq) have increased upon RepeatSoaking the data (Table 2), while the order of the enrichments remained virtually unchanged. Moreover, more stringent removal of reads overlapping low complexity regions further increased the significance of the functional enrichments. Correlation of gene expression from RNA-seq and microarray data was also increased. This observation suggests that aggressive removal of reads overlapping low complexity regions (RepeatSoaker threshold 0%) aids in emphasizing true biological signal within the data.
Filtering out low complexity reads affects low expressed metabolism-related genes
We investigated the effect of RepeatSoaker on the number of differentially expressed (DE) genes in RNA-seq data with and without duplicates (Figure 2). The total number of DE genes decreased with more stringent removal of reads with RepeatSoaker. Each condition, from no RepeatSoaker to the most stringent RepeatSoaker threshold 0%, has a unique set of genes not detected as differentially expressed in other conditions ("leaves" of the Venn diagram, Figure 2). These genes showed comparable fold change distributions as the "core" gene set detected as DE in all conditions (Figure 3A). However, the expression level of those genes was lower (Figure 3B), making them susceptible to lose their DE status upon removal of reads with RepeatSoaker. We further investigated the biological significance of those genes. These were predominantly predicted (as opposed to known) genes and metabolism-related genes, as can be seen from their enrichment analysis results (Additional File 7). This observation suggests that RepeatSoaker removes biological "noise" from the data while increasing the signal that reflects the underlying biology of the experiment.
Removing duplicates decreases overall expression level but retains fold changes
Lastly, we investigated the effect of data preprocessing steps upon the expression and fold change levels in RNA-seq data. Removing duplicates had decreased the overall level of expression (Figure 4A), as can be expected from losing ~40% of reads. The effect of other processing steps on gene expression level was negligible, which is reflected in the virtually unchanged Pearson's correlation coefficient of RNA-seq and microarray gene expression data.
To compare the effect of preprocessing steps on fold changes, we compared the distribution of fold changes at each step (Figure 4B), and investigated gene-by-gene fold change differences. As pre-processing steps were applied to both healthy and diseased groups of samples, thus uniformly changing read counts of the differentially expressed genes in both groups, the fold changes remained stable (Figure 4B). Expectedly, removing duplicates, with or without adapter trimming, decreased the maximum, but not overall, fold change, as would be expected from removing ~40% of reads. Similarly, filtering out reads overlapping and touching low complexity regions (RepeatSoaker threshold 0%) also slightly decreased the maximum fold change. Overall, our results suggest that data preprocessing steps, although affecting overall gene expression level, retain biological signal within the data, as reflected by relatively unchanged fold changes.
Our findings suggest some general guidelines for RNA-seq and ChIP-seq data processing. We show that each of the three steps, adapter trimming, duplicate removal, and filtering out reads overlapping low complexity regions, increases the strength of biological signals within the data, as shown by the more significant functional enrichment p-values. Our results emphasize the need of removing reads overlapping low complexity regions in order to improve biological signals in RNA-seq and ChIP-seq data. Overall, our analysis suggests that all three steps should be an integral part of NGS data processing pipelines in order to obtain better insights into the biology.
Although it is a general consensus that trimming the adapters improves the quality of sequencing data , our experience with applying adapter trimming to RNA-seq and ChIP-seq data has been mixed. Trimming the adapters did not significantly improve the significance of biological signals of RNA-seq data (Table 1), in contrast with what we observed for ChIP-seq data. This may be attributed to the use of different aligners, TopHat and Bowtie2, used for RNA-seq and ChIP-seq data, respectively. The former, TopHat, has been designed to deal with unmapped portions of the reads, such as adapters , and therefore the data processed with it may not be notably improved by the adapter trimming step. Our results warrant the testing of other adapter trimming tools, each reported to have different effect on data quality .
Our motivation to investigate the effect of reads overlapping low complexity regions came from our and others  empirical observations that such reads have multiple alignments within the reference genome and tend to pile up within low complexity regions, such as centromeres and telomeres. We hypothesized that such pileups may negatively affect the detection of true gene expression level when summarizing read counts into FPKM measurements. Furthermore, such pileups were picked up as the strongest peaks in ChIP-seq experiments, ultimately affecting motif enrichment analysis. To this end we have developed the RepeatSoaker tool that filters out reads overlapping low complexity regions. In addition to using low complexity regions defined by the RepeatMasker program , a user has an option to provide his/her own list of genomic coordinates of any other regions, such as the Duke excluded regions or the DAC blacklisted regions, defined by the ENCODE project and obtainable from the UCSC genome database , or the high-depth coverage regions defined by Pickrell et.al. by scanning the 1000 Genomes data . This ability of RepeatSoaker to use any list of genomic regions as a "mask" further empowers a user to ignore reads not only in the low complexity regions, but in other uninteresting regions, such as ribosomal genes in RNA-seq data, or even completely mask out non-exome portions of the genome.
One limitation of our study is that we cannot judge the biological significance of our findings other than by indirect assessment of the number of differentially expressed genes (RNA-seq), the total number of motifs (ChIP-seq), and the results of the enrichment analyses. Our hypothesis here was that, if a processing step improves the quality of the data, it should improve the significance of enriched gene ontologies, KEGG and Reactome pathways (for RNA-seq data) and detected motifs (for ChIP-seq data). Although we observed improvements in the enrichment analyses results after each data preprocessing steps, future studies warrant investigation of expression level of genes removed by data preprocessing steps, and by RepeatSoaker, by other techniques, such as polymerase chain reaction, or by direct comparison with gene expression changes measured by microarray technology.
In conclusion, we recommend adapter trimming, duplicates removal, and filtering out reads overlapping low complexity regions as data preprocessing steps of RNA-seq and ChIP-seq data. Our comprehensive comparison suggests that these data preprocessing steps will help to emphasize true biological signals in sequencing studies.
Huddleston J, Ranade S, Malig M, Antonacci F, Chaisson M, Hon L, et al: Reconstructing complex regions of genomes using long-read sequencing technology. Genome Res. 2014, 24 (4): 688-696. 10.1101/gr.168450.113.
Howorka S, Cheley S, Bayley H: Sequence-specific detection of individual DNA strands using engineered nanopores. Nat Biotechnol. 2001, 19 (7): 636-639. 10.1038/90236.
DNA sequencing costs. [http://www.genome.gov/sequencingcosts/]
Walsh P, Lu X, Carroll J: An Analysis of Next Generation Sequence Clipping Tools. Collaborative European Research Conference CERC 2013. 2013
Del Fabbro C, Scalabrin S, Morgante M, Giorgi FM: An extensive evaluation of read trimming effects on Illumina NGS data analysis. PLoS One. 2013, 8 (12): e85024-10.1371/journal.pone.0085024.
Turner S, Armstrong LL, Bradford Y, Carlson CS, Crawford DC, Crenshaw AT, et al: Quality control procedures for genome-wide association studies. Current protocols in human genetics / editorial board. Edited by: Jonathan L Haines. 2011, Chapter 1:Unit11.19:
DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, et al: A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011, 43 (5): 491-498. 10.1038/ng.806.
How PCR Duplicates Arise in Next-Generation Sequencing. [http://www.cureffi.org/2012/12/11/how-pcr-duplicates-arise-in-next-generation-sequencing/]
Abel HJ, Al-Kateb H, Cottrell CE, Bredemeyer AJ, Pritchard CC, Grossmann AH, et al: Detection of Gene Rearrangements in Targeted Clinical Next-Generation Sequencing. J Mol Diagn. 2014, 16 (4): 405-417. 10.1016/j.jmoldx.2014.03.006.
Chen Y, Negre N, Li Q, Mieczkowska JO, Slattery M, Liu T, et al: Systematic evaluation of factors influencing ChIP-seq fidelity. Nat Methods. 2012, 9 (6): 609-614. 10.1038/nmeth.1985.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al: The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009, 25 (16): 2078-2079. 10.1093/bioinformatics/btp352.
Furey TS: ChIP-seq and beyond: new and improved methodologies to detect and characterize protein-DNA interactions. Nat Rev Genet. 2012, 13 (12): 840-852. 10.1038/nrg3306.
Zhou W, Chen T, Zhao H, Eterovic AK, Meric-Bernstam F, Mills GB, Chen K: Bias from removing read duplication in ultra-deep sequencing experiments. Bioinformatics. 2014
Majewski J, Ott J: Distribution and characterization of regulatory elements in the human genome. Genome Res. 2002, 12 (12): 1827-1836. 10.1101/gr.606402.
DePristo MA, Zilversmit MM, Hartl DL: On the abundance, amino acid composition, and evolutionary dynamics of low-complexity regions in proteins. Gene. 2006, 378: 19-30.
Kapitonov VV, Jurka J: A universal classification of eukaryotic transposable elements implemented in Repbase. Nat Rev Genet. 2008, 9 (5): 411-412. 10.1038/nrg2165-c1.
Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN: RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics. 2010, 26 (4): 493-500. 10.1093/bioinformatics/btp692.
Chung D, Kuan PF, Li B, Sanalkumar R, Liang K, Bresnick EH, et al: Discovering transcription factor binding sites in highly repetitive regions of genomes with multi-read analysis of ChIP-Seq data. PLoS Comput Biol. 2011, 7 (7): e1002111-10.1371/journal.pcbi.1002111.
Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, et al: Initial sequencing and analysis of the human genome. Nature. 2001, 409 (6822): 860-921. 10.1038/35057062.
1000 Genomes Project Consortium, Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, et al: An integrated map of genetic variation from 1,092 human genomes. Nature. 2012, 491 (7422): 56-65. 10.1038/nature11632.
Hangauer MJ, Vaughn IW, McManus MT: Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs. PLoS Genet. 2013, 9 (6): e1003569-10.1371/journal.pgen.1003569.
Pickrell JK, Gaffney DJ, Gilad Y, Pritchard JK: False positive peaks in ChIP-seq and other sequencing-based functional assays caused by unannotated high copy number regions. Bioinformatics. 2011, 27 (15): 2144-2146. 10.1093/bioinformatics/btr354.
Medvedev P, Stanciu M, Brudno M: Computational methods for discovering structural variation with next-generation sequencing. Nat Methods. 2009, 6 (11): S13-S20.
Lee H, Popodi E, Foster PL, Tang H: Detection of structural variants involving repetitive regions in the reference genome. J Comput Biol. 2014, 21 (3): 219-233. 10.1089/cmb.2013.0129.
Schmieder R, Edwards R: Quality control and preprocessing of metagenomic datasets. Bioinformatics. 2011, 27 (6): 863-864. 10.1093/bioinformatics/btr026.
Olson LE, Soriano P: Increased PDGFRalpha activation disrupts connective tissue development and drives systemic fibrosis. Dev Cell. 2009, 16 (2): 303-313. 10.1016/j.devcel.2008.12.003.
Bolger AM, Lohse M, Usadel B: Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics. 2014, 30 (15): 2114-2120. 10.1093/bioinformatics/btu170.
Chen N: Using RepeatMasker to identify repetitive elements in genomic sequences. Curr Protoc Bioinformatics. Edited by: Andreas D Baxevanis. 2004, Chapter 4:Unit 4.10:
RepeatMasker Open-3.0. [http://www.repeatmasker.org]
FAQ: BED format. BED, [http://genome.ucsc.edu/goldenPath/help/customTrack.html]
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012, 7 (3): 562-578. 10.1038/nprot.2012.016.
Team RDC: R: A Language and Environment for Statistical Computing. 2013
Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10 (3): R25-10.1186/gb-2009-10-3-r25.
Feng J, Liu T, Zhang Y: Using MACS to identify peaks from ChIP-Seq data. Curr Protoc Bioinformatics. 2011, Chapter 2:Unit 2.14:
Machanick P, Bailey TL: MEME-ChIP: motif analysis of large DNA datasets. Bioinformatics. 2011, 27 (12): 1696-1697. 10.1093/bioinformatics/btr189.
Rosenbloom KR, Sloan CA, Malladi VS, Dreszer TR, Learned K, Kirkup VM, et al: ENCODE data in the UCSC Genome Browser: year 5 update. Nucleic Acids Res. 2013, 41 (Database issue): D56-D63.
We thank the two anonymous reviewers for their comments and critique. The content of this article is solely the responsibility of the authors and does not necessarily represent the official views of any agency involved in the funding of the work. The authors would like to acknowledge the National Institute of Arthritis and Musculoskeletal and Skin Diseases (a subaward from grant # P30 AR053483 to Dozmorov MG), an Institutional Development Award (IDeA) from the National Institute of General Medical Sciences (a subaward from grant # P30GM103510 to Dozmorov MG), the National Center for Advancing Translational Sciences (CTSA award No. UL1TR000058), NIH/NIGMS P30GM110766, NIH/NIGMS P20GM103456 and NIH/NIAMS R01AR065953 to Adrianto I, NIH/NIAMS P50AR060804, NIH/NIAMS R01AR065953 and NIH/NIGMS P20GM103456 to Lessard CJ, NIH/NIGMS P20GM103636 to Wren JD and Olson LE, and the National Science Foundation division of Advanced CyberInfrastructure # ACI-1345426 to Wren JD.
Funding for publication
Publication fees were paid primarily from NIH #P20GM103636.
This article has been published as part of BMC Bioinformatics Volume 16 Supplement 13, 2015: Proceedings of the 12th Annual MCBIOS Conference. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/16/S13.
The authors declare that they have no competing interests.
Conceived and designed the experiments: JDW, MGD. Provided the data: CM, KLS, LEO, TI, WMF, CJL. Analyzed the data: MGD, IA, CBG, SBG. Wrote the first draft of the manuscript: MGD. Made critical revisions and approved final version: JDW, IA, CL, CBG, EG. All authors reviewed and approved the final manuscript.
Electronic supplementary material
Additional File 1: Alignment statistics for the non-trimmed RNA-seq data in each data preprocessing step. Each worksheet contains step-specific alignment statistics. "no_remDup"/"remDup" in the worksheet name indicate whether the data has duplicates kept/removed, respectively. "All"/"0.75" etc., indicate whether the reads overlapping low complexity regions were kept/removed at the specified threshold, respectively (XLSX 148 KB)
Additional File 2: Alignment statistics for the trimmed RNA-seq data in each data preprocessing step. Each worksheet contains step-specific alignment statistics. "no_remDup"/"remDup" in the worksheet name indicate whether the data has duplicates kept/removed, respectively. "All"/"0.75" etc., indicate whether the reads overlapping low complexity regions were kept/removed at the specified threshold, respectively. (XLSX 149 KB)
Additional File 4: Differentially expressed genes in the non-trimmed RNA-seq data, and the results of KEGG/GO/Reactome pathway enrichment analyses identified in each data preprocessing step. Each worksheet contains step-specific alignment statistics, abbreviations are the same as in the Additional File 1. (XLSX 1 MB)
Additional File 5: Differentially expressed genes in the trimmed RNA-seq data, and the results of KEGG/GO/Reactome pathway enrichment analyses identified in each data preprocessing step. Each worksheet contains step-specific alignment statistics, abbreviations are the same as in the Additional File 2. (XLSX 1 MB)
Additional File 6: Motifs and their detection E-values identified in the ChIP-seq data in each data preprocessing step. (XLSX 70 KB)
Additional File 7: Investigation of the biology of genes lost after filtering out low complexity reads with RepeatSoaker. (PDF 3 MB)
About this article
Cite this article
Dozmorov, M.G., Adrianto, I., Giles, C.B. et al. Detrimental effects of duplicate reads and low complexity regions on RNA- and ChIP-seq data. BMC Bioinformatics 16 (Suppl 13), S10 (2015). https://doi.org/10.1186/1471-2105-16-S13-S10