Multiqc: AdapterRemoval output confusing vs FLASH
Description of bug:
We tested using AdapterRemoval and FLASH to merge reads for a CRISPR analysis pipeline and generated MultiQC outputs for both. It was confusing to see them side-by-side as AdapterRemoval looked like it had 2x the reads, but really I think some of the numbers weren't getting subtracted from one another.
AdapterRemoval output
It seems that AdapterRemoval is showing total reads when it should be showing:
- Full-length collapsed pairs (as is)
- uncollapsed pairs = Retained read pairs - Full length collapsed pairs
FLASH output (on same samples)
File that triggers the error:
Please drag and drop (and upload to the GitHub issue) an input file that I can use to replicate the error.
Here are some log files:
adapterremoval_logs.zip
flash_logs.zip
MultiQC run details (please complete the following):
- Command used to run MultiQC:
multiqc -f $rtitle $rfilename --config $multiqc_config . -m fastqc -m adapterRemoval -m flash - MultiQC Version: MultiQC v1.7
- Operating System: Ubuntu
- Python Version:
(base) root@671e64af0899:/# /opt/conda/envs/nf-core-crisprvar-1.0dev/bin/python --version
Python 2.7.15
- Method of MultiQC installation: conda, this is the environment.yml
olgabot
All 15 comments
Ok, so if I understand you correctly, we should just not plot Retained read pairs? As this is already a sum of the other categories?
ewels
on 18 Aug 2019
It could be, however the plot is indicating number of "Full-length collapsed pairs" and "uncollapsed pairs," and "singleton R{1,2}" are all separate. Here's some log files run on the same sample for reference.
FLASH output
Here's the critical output:
[FLASH] Read combination statistics:
[FLASH] Total pairs: 67277
[FLASH] Combined pairs: 66658
[FLASH] Innie pairs: 65235 (97.87% of combined)
[FLASH] Outie pairs: 1423 (2.13% of combined)
[FLASH] Uncombined pairs: 619
[FLASH] Percent combined: 99.08%
Total Pairs > Combined pairs
67 277 > 66 658
Full output:
Sun 18 Aug - 12:07 î‚° ~ î‚°
î‚° grep FLASH /Users/olgabot/Downloads/mNGplate3/crispresso2_vanilla/cripresso/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/CRISPResso_RUNNING_LOG.txt
Running FLASH command: flash mNGplate3_sorted_A1_POLR1A-C_S97_R1_001.fastq.gz mNGplate3_sorted_A1_POLR1A-C_S97_R2_001.fastq.gz --min-overlap 10 --max-overlap 100 --allow-outies -z -d /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001 -o out >>/tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/CRISPResso_RUNNING_LOG.txt 2>&1
[FLASH] Starting FLASH v1.2.11
[FLASH] Fast Length Adjustment of SHort reads
[FLASH]
[FLASH] Input files:
[FLASH] mNGplate3_sorted_A1_POLR1A-C_S97_R1_001.fastq.gz
[FLASH] mNGplate3_sorted_A1_POLR1A-C_S97_R2_001.fastq.gz
[FLASH]
[FLASH] Output files:
[FLASH] /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/out.extendedFrags.fastq.gz
[FLASH] /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/out.notCombined_1.fastq.gz
[FLASH] /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/out.notCombined_2.fastq.gz
[FLASH] /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/out.hist
[FLASH] /tmp/nxf.BJPjRHsyZR/mNGplate3_sorted_A1_POLR1A-C/CRISPResso_on_mNGplate3_sorted_A1_POLR1A-C_S97_R1_001_mNGplate3_sorted_A1_POLR1A-C_S97_R2_001/out.histogram
[FLASH]
[FLASH] Parameters:
[FLASH] Min overlap: 10
[FLASH] Max overlap: 100
[FLASH] Max mismatch density: 0.250000
[FLASH] Allow "outie" pairs: true
[FLASH] Cap mismatch quals: false
[FLASH] Combiner threads: 64
[FLASH] Input format: FASTQ, phred_offset=33
[FLASH] Output format: FASTQ, phred_offset=33, gzip
[FLASH]
[FLASH] Starting reader and writer threads
[FLASH] Starting 64 combiner threads
[FLASH] Processed 25000 read pairs
[FLASH] Processed 50000 read pairs
[FLASH] Processed 67277 read pairs
[FLASH]
[FLASH] Read combination statistics:
[FLASH] Total pairs: 67277
[FLASH] Combined pairs: 66658
[FLASH] Innie pairs: 65235 (97.87% of combined)
[FLASH] Outie pairs: 1423 (2.13% of combined)
[FLASH] Uncombined pairs: 619
[FLASH] Percent combined: 99.08%
[FLASH]
[FLASH] Writing histogram files.
[FLASH] WARNING: An unexpectedly high proportion of combined pairs (99.70%)
increasing this parameter. (As-is, FLASH is penalizing overlaps longer than
[FLASH]
[FLASH] FLASH v1.2.11 complete!
[FLASH] 1.611 seconds elapsed
[FLASH] Finished with 1 warning (see above)
AdapterRemoval output
[Trimming statistics]
Total number of read pairs: 67277
Number of unaligned read pairs: 385
Number of well aligned read pairs: 66892
Number of discarded mate 1 reads: 0
Number of singleton mate 1 reads: 11
Number of discarded mate 2 reads: 11
Number of singleton mate 2 reads: 0
Number of reads with adapters[1]: 2878
Number of full-length collapsed pairs: 66849
Number of truncated collapsed pairs: 0
Number of retained reads: 67694
Number of retained nucleotides: 22884082
Average length of retained reads: 338.052
Number of retained reads > Total number of read pairs > Number of well aligned read pairs > Number of full-length collapsed pairs
67 694 > 67 277 > 66 892 > 66 849
Sun 18 Aug - 12:07 î‚° ~ î‚°
î‚° head -n 45 /Users/olgabot/Downloads/mNGplate3/adapterremoval/adapterremoval/mNGplate3_sorted_A1_POLR1A-C.settings
AdapterRemoval ver. 2.3.0
Trimming of paired-end reads
[Adapter sequences]
Adapter1[1]: AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG
Adapter2[1]: AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
[Adapter trimming]
RNG seed: 4242927527
Alignment shift value: 2
Global mismatch threshold: 0.333333
Quality format (input): Phred+33
Quality score max (input): 41
Quality format (output): Phred+33
Quality score max (output): 41
Mate-number separator (input): '/'
Trimming 5p: 0
Trimming 3p: 0
Trimming Ns: Yes
Trimming Phred scores <= 2: Yes
Trimming using sliding windows: No
Minimum genomic length: 15
Maximum genomic length: 4294967295
Collapse overlapping reads: Yes
Minimum overlap (in case of collapse): 11
[Trimming statistics]
Total number of read pairs: 67277
Number of unaligned read pairs: 385
Number of well aligned read pairs: 66892
Number of discarded mate 1 reads: 0
Number of singleton mate 1 reads: 11
Number of discarded mate 2 reads: 11
Number of singleton mate 2 reads: 0
Number of reads with adapters[1]: 2878
Number of full-length collapsed pairs: 66849
Number of truncated collapsed pairs: 0
Number of retained reads: 67694
Number of retained nucleotides: 22884082
Average length of retained reads: 338.052
Seems like they match in their "Total number of read pairs" (67 277) and 2 significant figures for "full-length collapsed pairs" (66 658 for FLASH, 66 849 for AdapterRemoval), so I'd be inclined to think those are comparable. It doesn't seem like for this report, AdapterRemoval is treating the separate R1 and R2 of paired-end reads separately.
Anyway, I can make the PR to change this if this is making any sense.
olgabot
on 18 Aug 2019
if this is making any sense.
Not really, but I have a crying baby in the background which is making it difficult to think ðŸ˜
If you think you know how to fix it then, yeah - please put in a PR! That would be fab 😀
ewels
on 18 Aug 2019
@alexanderscholz - do you have any feedback on this issue? You wrote the AdapterRemoval MultiQC module. I'll look into resolving this now but would be good to have your thoughts.
Phil
ewels
on 27 May 2020
@jfy133 - we're using AdapterRemoval v2 in https://github.com/nf-core/eager too - what are your thoughts on this one?
apeltzer
on 27 May 2020
@alexhbnr and I have already noticed that but Alex hasn't had time to fix and I still haven't learnt python ;).
We had decided the following: Retained reads should definitely be a separate plot. This is the only metric that matches the actual number of reads in the output file.
All other metrics can go together as they represent the actual processing of the reads themselves.
jfy133
on 27 May 2020
I submitted the pull request (#1204) with my initial changes and am happy if someone wants to think about the implemented changes and improve them.
It's a poor man's fix since it doesn't pull out the retained reads in its own plot but instead treats them as reads that haven't been changed, i.e. no singletons or collapsed reads. This get rid of the problem that the length of the bar plots are twice as long as they should be but more or less sacrifices the retained reads metric.
alexhbnr
on 27 May 2020
Nice, thanks @alexhbnr! The plot definitely looks a lot better, I'm sure that it is a step in the right direction. However the comment from @olgabot makes me think that we can't combine many of these categories together in a single stacked bar plot.
This is similar to work I was doing a couple of days in the Cutadapt module for e7ed1052b207c2396e69438b1134a39ed7358afa - I parsed a tonne of fields and then realised that I couldn't plot most of them as the categories didn't stack properly.
Here I'm not even sure that the % reads trimmed in the General Statistics is really trustworthy, given that @olgabot's example has _more_ retained reads (67694) than were inputted (67277)..
I also don't think that the % Trimmed Reads in this sense is actually very interesting. What I would think that people actually want to know is what percentage of their data was lost. For Cutadapt I did this as the % of bases removed.
Phil
ewels
on 28 May 2020
@MikkelSchubert - hello! As author of AdapterRemoval, I was wondering if you could help us?
We currently have a MultiQC module for AdapterRemoval, but we think that the bar plot it produces is wrong / misleading in the way that it stacks categories - see an example report: multiqc_report.html
Are you able to shed some light for us on how the different categories reported by AdapterRemoval are related, and which we would be able to use to create a stacked bar plot representing the total read count?
If you have any views on what would be best reported in MultiQC, that would also be very helpful.
Many thanks in advance,
Phil
ewels
on 28 May 2020
Hey Phil - just some clarifications (but I'm sure Mikkel will be better here), the inputted reads you've written above is _read pairs_, whereas the retained reads represent each unique entry in the FASTQ file i.e. will include merged reads, but also unmerged reads as a form of 'singletons', which is why you can have more (because 1 un-merged read pair turns in 1).
The way I would personally plot this would be as follows (based off of Olga's stats table), if I understand the statistics correctly:
[Trimming statistics]
Input
Maybe just for general stats for validation e.g. compared to FASTQC
Total number of read pairs: 67277
Merging
stacked bar; could also have percentage option
Number of unaligned read pairs: 385
Number of well aligned read pairs: 66892
[This basically reports the number of input reads pairs that overlap, prior 'collapsing'; if the input is paired end and 'collapsing' is selected by the user - therefore unaligned read pairs are those that did not overlap so the insert size is bigger than the sequencing cycles]
Quality filtering
stacked bar; percentage option wouldn't make sense here - see comment
Number of full-length collapsed pairs: 66849 (merged reads, no adapter so untrimmed)
Number of truncated collapsed pairs: 0 (merged reads, adapter trimmed)
Number of singleton mate 1 reads: 11 (unmerged forward reads)
Number of singleton mate 2 reads: 0 (unmerged reverse reads)
Number of discarded mate 1 reads: 0 (low-quality removed reads, e.g. length filter)
Number of discarded mate 2 reads: 11 (low-quality removed reads, e.g. length filter)
[This would represent the categories that could be included as single entries in the 'Retained Reads category; however, this totals is 66860, meaning we are still missing some reads, something I've never been able to work out - only thing I can think of is that these are single reads without one of a pair in the input FASTQ file itself; but I'm not sure if that should really exist coming off demultiplexing ]
Output reads
barplot
Number of retained reads: 67694
jfy133
on 28 May 2020
I agree with @jfy133, separate plots would make more sense. However, for the suggested Quality filtering, the six categories included into it (full-length collapsed, truncated collapsed, singleton mate 1 and 2, and discarded mate 1 and 2) do not sum up to the number of retained reads.
For the example above, retained reads (67694) > collapsed reads (66849) + singletons (22).
So it would make sense to get feedback from @MikkelSchubert how he actually calculates these.
alexhbnr
on 28 May 2020
Would could always calculate the difference and make a new plot category called "Mystery reads" 😆
ewels
on 28 May 2020
I apologize for the confusing output from AdapterRemoval; it is something I've been wanting to rework, but never got around to. It is not clear to me exactly what you want to show, so I'll just try to explain the output you are seeing:
As has been noted, the numbers below show either pairs or individual reads. You need to keep this in mind, since a pair accounts for two individual reads.
Total number of read pairs: 67277
The total number of input read pairs.
Number of unaligned read pairs: 385
Number of well aligned read pairs: 66892
Statistics about the alignment pair-wise alignments; the sum of these numbers correspond to the number of (input) reads.
Number of discarded mate 1 reads: 0
Number of singleton mate 1 reads: 11
Number of discarded mate 2 reads: 11
Number of singleton mate 2 reads: 0
Statistics for discarded reads and "singleton" reads; this gets a bit complicated, especially when you have merged reads. Firstly, a singleton read is a read in which the mate (but not the read itself) was discarded, meaning that a pair that was discarded increments both "discarded mate 1" and "discarded mate 2", while a pair with one discarded mate increments "discarded mate A" and "singleton mate B", where A != B. In addition, merged reads that are discarded are counted as if the pair was discarded (increments both "discarded mate" values).
Number of reads with adapters[1]: 2878
This number is the number of reads with adapter sequence that got trimmed; for normal paired reads (where the read lengths are identical) this will always be an even number.
Number of full-length collapsed pairs: 66849
Number of truncated collapsed pairs: 0
The number of merged reads, with "truncated" pairs being those merged reads that were trimmed following merging (due to low quality bases or Ns). Note that these numbers do not include discarded merged reads.
Number of retained reads: 67694
Number of retained nucleotides: 22884082
The total number of output reads (and how many nucleotides that amounts to); a simple way to calculate this is 2 * input - sum(discarded) - sum(merged), since each input pair accounts for 2 reads, the discarded reads are excluded, and merging transforms two input reads into one output read; 2 * 67277 - (11 + 0) - (66849 + 0) = 67694.
As has been noted, this matches the sum of reads in the output files (excluding discarded reads), where a merged pair has been transformed into a single read.
MikkelSchubert
on 28 May 2020
Many thanks @MikkelSchubert - this is brilliant! Any chance that this text could be put in the AdapterRemoval docs somewhere (https://adapterremoval.readthedocs.io/en/latest/)? I'm happy to open a PR pasting your issue text if you'd like.
Need to do some head scratching now and figure out is the most important information / what people actually want to know in the MultiQC report. @olgabot - what were you interested in? I guess the number of reads that were collapsed?
Phil
ewels
on 28 May 2020
In addition - I have a request from @aidaanva that the percentage of reads could also be calculated and reported specifically in the general stats table (as gives you an idea of how well your PE data has merged - this is a good indicator for ancient DNA)
I guess this would be Number of well aligned read pairs / Total number of read pairs, if I understand Mikkel correctly. @apeltzer is that how you calculated in ReportTable.
jfy133
on 2 Jun 2020
Related issues
berguner
·
5Comments
anoronh4
·
5Comments
choudharis2
·
7Comments
adirot
·
5Comments
cgirardot
·
4Comments
Most helpful comment
I apologize for the confusing output from AdapterRemoval; it is something I've been wanting to rework, but never got around to. It is not clear to me exactly what you want to show, so I'll just try to explain the output you are seeing:
As has been noted, the numbers below show either pairs or individual reads. You need to keep this in mind, since a pair accounts for two individual reads.
The total number of input read pairs.
Statistics about the alignment pair-wise alignments; the sum of these numbers correspond to the number of (input) reads.
Statistics for discarded reads and "singleton" reads; this gets a bit complicated, especially when you have merged reads. Firstly, a singleton read is a read in which the mate (but not the read itself) was discarded, meaning that a pair that was discarded increments both "discarded mate 1" and "discarded mate 2", while a pair with one discarded mate increments "discarded mate A" and "singleton mate B", where A != B. In addition, merged reads that are discarded are counted as if the pair was discarded (increments both "discarded mate" values).
This number is the number of reads with adapter sequence that got trimmed; for normal paired reads (where the read lengths are identical) this will always be an even number.
The number of merged reads, with "truncated" pairs being those merged reads that were trimmed following merging (due to low quality bases or Ns). Note that these numbers do not include discarded merged reads.
The total number of output reads (and how many nucleotides that amounts to); a simple way to calculate this is 2 * input - sum(discarded) - sum(merged), since each input pair accounts for 2 reads, the discarded reads are excluded, and merging transforms two input reads into one output read; 2 * 67277 - (11 + 0) - (66849 + 0) = 67694.
As has been noted, this matches the sum of reads in the output files (excluding discarded reads), where a merged pair has been transformed into a single read.