How to Trim Sequencing Reads Without Distorting Your Own Data

By the end of this you will have a set of trimmed FASTQ files, a per-sample QC report you can diff against the pre-trim report, and a recorded decision about how aggressive your trimming was and why. You need the raw FASTQs (gzipped, paired-end in most cases), the library prep and instrument model from whoever generated the data, roughly 3x the raw file size in free disk, and a conda or container environment with fastqc, fastp, cutadapt, multiqc, and whichever aligner or quantifier you plan to use. On a laptop, a 30x whole-genome FASTQ pair (around 100 GB gzipped) will take hours to stream through any trimmer, so plan on a machine with at least 8 cores and fast local disk. RNA-seq at 30-50M read pairs per sample is a few minutes per sample.

Before anything else: trimming is a data-cleaning step, not a measurement in itself. Nothing here is a clinical result, and no output of this pipeline tells you anything about your health. Variant calls and expression estimates from personal data need a clinician and, for anything actionable, a CLIA-confirmed test.

1. Inventory the FASTQs and confirm what you have

Read the first few records rather than trusting the filename.

zcat sample_R1.fastq.gz | head -8
zcat sample_R1.fastq.gz | awk 'NR%4==2 {print length($0)}' | head -100000 | sort -n | uniq -c | tail

Three things you want out of this: read length (fixed 150 bp, or already variable because someone trimmed before you), the instrument ID in the header line, and whether quality encoding is modern Phred+33. Anything sequenced in the last decade is Phred+33; if you see characters below !-I range oddities, stop and check.

Instrument matters more than people expect. NovaSeq, NextSeq, and iSeq use two-color chemistry, where “no signal” is called as G. Reads that run past the end of the insert produce long runs of high-quality-looking G, and a generic quality trimmer will not remove them because their Phred scores are fine. On four-color instruments (HiSeq 2500, MiSeq) you get poly-A/poly-N artifacts instead. You need to know which case you are in before picking flags.

Also write down the library chemistry. Whether the RNA was poly-A selected or rRNA-depleted, and which extraction kit was used, changes the composition of the library before a single base is trimmed, including duplication rate and the fraction of intronic and intergenic reads ¹. Trimming will not fix a prep artifact, and mistaking one for the other wastes a day.

2. Run FastQC first and read the adapter panel, not just the traffic lights

fastqc -t 8 -o qc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz

Open the HTML and go straight to four modules:

• Per base sequence quality: where does the median cross Q20, and on which read? Read 2 degrading faster than read 1 is normal.
• Adapter content: does it rise at all, and where does it start climbing?
• Overrepresented sequences: FastQC will often name the adapter for you outright.
• Per base sequence content: a wobbly first 10-12 bases is the random-hexamer priming signature in most RNA-seq libraries, and it is not a quality problem.

The red X on “Per base sequence content” for RNA-seq is expected and you should not hard-trim 12 bases off the 5’ end to make it green. That bias is positional nucleotide composition, not error, and cutting it costs you real signal.

3. Determine the adapter sequence instead of guessing

Most libraries are TruSeq-family, whose adapters share the prefix AGATCGGAAGAGC, or Nextera/tagmentation, which uses CTGTCTCTTATACACATCT. Confirm rather than assume:

zcat sample_R1.fastq.gz | head -4000000 | awk 'NR%4==2' \
  | grep -c "AGATCGGAAGAGC"
zcat sample_R1.fastq.gz | head -4000000 | awk 'NR%4==2' \
  | grep -c "CTGTCTCTTATACACATCT"

Whichever count is materially larger tells you the family. If both are near zero and FastQC shows no adapter content, your insert sizes are comfortably longer than the read length and you have almost no read-through. That is a good position to be in, and it means adapter trimming will be close to a no-op.

For paired-end data you can skip this step entirely by using overlap-based detection, which is what we do in practice (next step). Overlap detection finds the adapter from the data by aligning R1 to the reverse complement of R2 and calling everything past the overlap as adapter. It is more reliable than a hardcoded FASTA because it catches non-standard indexes and custom primers.

4. Trim with fastp, gently

This is our default for both DNA and RNA:

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o trim/sample_R1.fq.gz -O trim/sample_R2.fq.gz \
  --detect_adapter_for_pe \
  --cut_tail --cut_tail_window_size 4 --cut_tail_mean_quality 15 \
  --trim_poly_g \
  --length_required 36 \
  --n_base_limit 5 \
  --thread 8 \
  --json qc/sample.fastp.json --html qc/sample.fastp.html

What each choice is doing, and why we set it that way:

--detect_adapter_for_pe turns on the overlap-based adapter inference described above. It is off by default for paired-end input in some fastp versions, which surprises people.

--cut_tail with a 4 bp window and mean Q15 is a sliding-window trim from the 3’ end only. We do not use --cut_front on RNA-seq. Quality at the 5’ end is rarely the problem, and trimming there costs alignment anchor.

--trim_poly_g removes the two-color dark-cycle artifact. fastp enables this automatically when it detects a NextSeq/NovaSeq instrument ID, but setting it explicitly documents the intent. On four-color data, use --trim_poly_x instead if you see poly-A tails bleeding in.

--length_required 36 drops reads that trimming shortened below usable length. For splice-aware alignment we would rather discard a 20 bp fragment than let it map ambiguously across the transcriptome. For genome alignment with BWA-MEM you can lower this to 25-30.

Note what is absent: no global quality filter above Q15, no fixed-length head crop, no hard MINLEN 50. This is deliberate. Aggressive quality trimming changes gene expression estimates, and the change is not uniform across genes. Williams et al. compared trimming thresholds across a range and found that trimming at higher Phred cutoffs progressively reduced the number of genes with accurate FPKM estimates, with gentle trimming (around Phred 5) or no quality trimming at all giving the most accurate quantification ². The mechanism is mundane: low-quality bases at read ends get soft-clipped by modern aligners anyway, but shortened reads become harder to place uniquely, and short-transcript and low-expression genes take the hit disproportionately.

The broader evaluation literature agrees on the shape of the tradeoff. Del Fabbro et al. tested trimming across assembly, alignment, and SNP calling and found quality trimming generally helps, but that the benefit saturates and then reverses as thresholds climb, with the optimum depending on the downstream application rather than on the FASTQ alone ³. So the answer to “what threshold” is: mild, and different for assembly than for quantification.

If you are running de novo transcriptome assembly rather than quantification against a reference, trim harder. Assemblers build k-mer graphs, and a single erroneous base spawns a spurious branch. Trimmomatic with SLIDINGWINDOW:4:20 MINLEN:36 ahead of Trinity is the well-trodden path there and the one we would use ⁴.

5. If you must use Trimmomatic, here is the equivalent

Trimmomatic is still the most-cited trimmer and you will inherit pipelines built on it. It runs steps in the order you list them, which is a real footgun: put ILLUMINACLIP first, always, because quality trimming a read before removing its adapter can leave an adapter fragment too short to match.

trimmomatic PE -threads 8 -phred33 \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  trim/sample_R1.fq.gz trim/sample_R1.unpaired.fq.gz \
  trim/sample_R2.fq.gz trim/sample_R2.unpaired.fq.gz \
  ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2:True \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

ILLUMINACLIP:fasta:2:30:10 is seed mismatches:palindrome clip threshold:simple clip threshold. The two extra fields, :2:True, enable palindrome mode keeping both reads of the pair, which matters because the default silently discards the reverse read when read-through is detected. Many old tutorials omit them and then people wonder why R2 shrank.

Trimmomatic is slower than fastp, single-purpose, and produces no QC report. fastp does adapter removal, quality trimming, poly-G handling, length filtering, and before/after QC in one pass over the data, which for a 100 GB WGS FASTQ is the difference between one I/O pass and three. That is the main reason we default to it. Tools in the same generation, such as AfterQC, made the same argument and added per-cycle bias profiling and overlap-based error correction for paired reads ⁵.

6. Verify with a pre/post diff, not a glance

fastqc -t 8 -o qc_trim/ trim/sample_R1.fq.gz trim/sample_R2.fq.gz
multiqc -o qc_report/ qc_raw/ qc_trim/ qc/

Read the fastp JSON directly for the numbers that matter:

jq '{before: .summary.before_filtering, after: .summary.after_filtering,
     filtering: .filtering_result, adapter: .adapter_cutting.adapter_trimmed_reads}' \
   qc/sample.fastp.json

Sanity thresholds we apply to RNA-seq:

• Reads passing filter below 90%: investigate before proceeding. Something is wrong upstream.
• Reads passing filter above 99.5% with visible adapter content in raw FastQC: your adapter detection failed. Check step 3.
• Mean read length dropping by more than about 15% of the original: too aggressive. Loosen the quality cutoff.
• too_short_reads above a few percent: your insert size distribution is shorter than your read length, which is a library prep issue, not something trimming solves.

Then check the only metric that is truly decisive: alignment rate. Align a 2M-read subsample both trimmed and untrimmed and compare.

seqtk sample -s42 sample_R1.fastq.gz 2000000 > sub_raw_R1.fq
# ... same for R2 and trimmed ...
STAR --genomeDir idx --readFilesIn sub_raw_R1.fq sub_raw_R2.fq \
     --runThreadN 8 --outSAMtype None --outFileNamePrefix raw_

If trimming moves uniquely-mapped rate by less than a percentage point, trimming is doing nothing for you and you should prefer the gentler setting. If it moves the rate up by five points, you had real adapter contamination and the trimming earned its place.

7. Adjust per downstream application

Variant calling from WGS: adapter removal yes, quality trimming minimal. BWA-MEM soft-clips, and GATK’s base quality score recalibration already models per-cycle and per-context error. Trimming low-quality tails before BQSR removes the very observations BQSR uses to learn the error model. We use adapter and poly-G removal plus --length_required 30 and nothing else.

RNA-seq quantification with Salmon or kallisto: selective alignment handles soft-clipping, so the benefit of quality trimming is small and the risk of shifting estimates is real ². Adapter removal still matters, because adapter sequence is a k-mer that does not exist in your transcriptome index and drags down mapping rate.

Amplicon and 16S data: a different regime entirely. Fixed-length truncation determines whether R1 and R2 still overlap enough to merge, and truncating too far destroys merging while truncating too little drags low-quality bases into the denoiser. FIGARO makes this an optimization problem over expected error and merge rate rather than a guess ⁶. If you are running DADA2, use it.

De novo assembly: harder trimming, as above ⁴.

8. Pin the pipeline and record the decision

Write the exact command and tool versions into the output directory:

fastp --version 2>&1 | tee trim/TRIM_PROVENANCE.txt
echo "$FASTP_CMD" >> trim/TRIM_PROVENANCE.txt
md5sum trim/*.fq.gz >> trim/TRIM_PROVENANCE.txt

For longitudinal data this is the whole ballgame. If you sequence your transcriptome again in eight months and trim the new FASTQs with a different fastp version or a different quality cutoff, any difference you see between timepoints has a trimming component you cannot separate from biology. Freeze the trimming step, then re-run the old samples through the new pipeline whenever you change it. Same rule as any other reproducible build.

Common problems

Adapter content stays flat in post-trim FastQC but mapping rate did not move. You trimmed an adapter that was not there. FastQC’s adapter module only searches a small set of known sequences, so a custom or single-cell prep adapter shows up as an overrepresented sequence instead. Grep for the top overrepresented sequence in your raw reads and pass it explicitly with --adapter_sequence.

R2 output is much smaller than R1 after Trimmomatic. Palindrome mode discarded reverse reads on read-through pairs. Add :2:True to the ILLUMINACLIP argument as shown above.

Reads still contain long G runs after --trim_poly_g. The poly-G trimmer requires a minimum run length (default 10). A read that is half real sequence and half G will be caught; a read with a 6-base G tail will not. Lower with --poly_g_min_len 6 if your post-trim per-base content shows it.

Duplication rate is 60% and you want to trim your way out. You cannot. High duplication in RNA-seq usually comes from low input, high PCR cycles, or a highly skewed transcriptome, all decided before sequencing ¹. Trimming changes nothing here, and deduplicating RNA-seq without UMIs throws away real high-expression signal.

Your DE results changed after you changed trimming parameters. They will. Trimming shifts read counts non-uniformly across genes, and a downstream normalization step then redistributes that shift across the whole sample ². Choose a trimming policy once, apply it to every sample in the comparison, and keep it fixed. If the biological conclusion only survives one trimming setting, it is not a conclusion. The normalization and DE method is a separate lever with its own sensitivities, and you should treat those choices as explicitly as the trimming ones ⁷⁸.

Someone tells you to “trim the data” and means the expression matrix, not the reads. There is a second, unrelated use of the word: trimming extreme values from a normalized expression matrix before fitting a model. In machine learning on omics data, discarding a fraction of the most extreme feature values can improve cross-dataset generalization ⁹. Useful technique, entirely different step, downstream of everything in this guide. Make sure you know which one is being discussed before you change a pipeline.

You have no idea whether any of this improved the biology. Fair. The sequence of checks that resolves it: raw FastQC, trimmed FastQC, subsample alignment rate, then gene body coverage and the fraction of reads in exons from your aligner’s log. If all four are stable across trimming settings, stop tuning and move to quantification ¹⁰.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Marc Sultan, Vyacheslav Amstislavskiy, Thomas Risch, et al. Influence of RNA extraction methods and library selection schemes on RNA-seq data. BMC Genomics, 2014. https://doi.org/10.1186/1471-2164-15-675 ↩ ↩²

2. Claire R. Williams, Alyssa Baccarella, Jay Z. Parrish, et al. Trimming of sequence reads alters RNA-Seq gene expression estimates. BMC Bioinformatics, 2016. https://doi.org/10.1186/s12859-016-0956-2 ↩ ↩² ↩³

3. Cristian Del Fabbro, Simone Scalabrin, Michele Morgante, et al. An Extensive Evaluation of Read Trimming Effects on Illumina NGS Data Analysis. PLoS ONE, 2013. https://doi.org/10.1371/journal.pone.0085024 ↩

4. Steven O. Sewe, Gonçalo Silva, Paulo Sicat, et al. Trimming and Validation of Illumina Short Reads Using Trimmomatic, Trinity Assembly, and Assessment of RNA-Seq Data. Methods in Molecular Biology, 2022. https://doi.org/10.1007/978-1-0716-2067-0_11 ↩ ↩²

5. Shifu Chen, Tanxiao Huang, Yanqing Zhou, et al. AfterQC: automatic filtering, trimming, error removing and quality control for fastq data. BMC Bioinformatics, 2017. https://doi.org/10.1186/s12859-017-1469-3 ↩

6. Michael M. Weinstein, Aishani Prem, Mingda Jin, et al. FIGARO: An efficient and objective tool for optimizing microbiome rRNA gene trimming parameters. 2019. https://doi.org/10.1101/610394 ↩

7. Zihan Cui, Yuhang Liu, Jinfeng Zhang, et al. Super-delta2: an enhanced differential expression analysis procedure for multi-group comparisons of RNA-seq data. Bioinformatics, 2021. https://doi.org/10.1093/bioinformatics/btab155 ↩

8. Jiung-Wen Chen, Lisa Shrestha, George Green, et al. The hitchhikers’ guide to RNA sequencing and functional analysis. Briefings in Bioinformatics, 2023. https://doi.org/10.1093/bib/bbac529 ↩

9. Victor Tkachev, Maxim Sorokin, Constantin Borisov, et al. Flexible Data Trimming Improves Performance of Global Machine Learning Methods in Omics-Based Personalized Oncology. International Journal of Molecular Sciences, 2020. https://doi.org/10.3390/ijms21030713 ↩

10. Canan Külahoglu, Andrea Bräutigam. Quantitative Transcriptome Analysis Using RNA-seq. Methods in Molecular Biology, 2014. https://doi.org/10.1007/978-1-4939-0700-7_5 ↩