Multiqc: Theoretical distribution for GC content

[mikelove]

This is an open-ended feature request. I'm suggesting we specify a format for describing a theoretical distribution which could be plotted over FASTQC's Per Sequence GC Content module.

I'm willing to contribute code separately as an R package to generate this distribution and write it to the appropriate location in the appropriate format.

[ewels]

This should be very doable. My suggestion would be to approach it in a few ways:

1. Have a new file format detected - eg. fastqc_theoretical_gc.txt
   
   
   
   • MultiQC would find and parse this if it's there, then plot as a single background line.
   
   
   • If multiple files are found it would overwrite, using the last one found
   
   
   • Would probably make sense to use the same data format as FastQC - example here
   
   
2. Bundle a few common curves with MultiQC
   
   
   
   • eg. Human / Mouse, genome / transcriptome
   
   
   • Create a new MultiQC config group - fastqc_config > theoretical_gc where the desired line would be specified.
   
   
   • This could be global or in the working directory multiqc_config.yaml
   
   
   • A fastqc_theoretical_gc.txt file would take priority over this if found.
   
   

Does this sound sensible? This way people don't have to run any scripts if they're using a common organism, but it'd be easy to run your new script / package thing to have custom lines.

@mikelove - do you think you could generate those common curves for me? 😁

Phil

[ewels]

ps. For bonus points, if we can put this together quickly then it'll ship in the v0.9 release tomorrow 😉

[mikelove]

Sure not a problem. I can generate human and mouse 0-100 % GC for the most recent Ensembl UCSC genome and transcriptome. These curves don't change much across versions. I'll have to think about what subset of the transcriptome would be most appropriate (e.g. for mRNA-seq, which is what I see most often). And it will be quite easy for someone who wants something else to generate it via a little R package + CL script I write.

[mikelove]

ha! ok I'll see how much i can do today :)

[mikelove]

btw, i've been promoting MultiQC at every rna-seq talk I give. "do you use MultiQC? no? ok, you need to use this from now on."

[ewels]

Sounds great! It should be pretty easy for me to add to MultiQC, so I'll get on it as soon as I get some data files from you.

Thanks for the promotion! I'm running out of excited emoji to add to your comments..

[mikelove]

Wondering how to scale the counts (not a problem for percent, but if user switches tab...). A reasonable default would be 30e6 of 100bp reads for RNA-seq. What about for whole genome seq?

A little more complex would be for MultiQC to take the preset counts and scale them to the median of the # of reads across the samples.

[ewels]

The lazy option would be to just show it for the percentage tab. But yeah, should be simple enough to scale for counts to the sample with the highest read count. Better than trying to guess people's read count numbers - that's guaranteed to mess up someone's report somewhere.

[ewels]

So yeah - a single file with the percentage of the genome within each GC percentage bin would be ideal I think.

[mikelove]

Ok will have in the next hour.

[ewels]

Great - I'll start working now using this as an example for the parsing code: fastqc_theoretical_gc.txt

(Will switch out with your data when available).

[ewels]

Ok, some starting code pushed. Only change in the file format is that I added a header line with the reference name, as follows:

# FastQC theoretical GC content curve: Human Genome
0   0.005311768
1   0.004108502
2   0.004060371
[..]

(see test file).

No built-in reference file stuff yet, only parsing of any files called fastqc_theoretical_gc.txt. Can add in the built-in references tomorrow (it's getting late here and I need to head home).

Shout if you find any bugs or would like it to behave differently!

Note that it doesn't work on the secondary _Counts_ plot yet - I need to add to the core Additional data series code to make this work.

Phil

[mikelove]

Ok gotcha. I'll try it out.

[mikelove]

I've written an R package which can be installed with

devtools::install_github("mikelove/fastqcTheoreticalGC")

and I generated some example files here:

https://github.com/mikelove/fastqcTheoreticalGC/tree/master/inst/script

And the code used to generate those files is here:

https://github.com/mikelove/fastqcTheoreticalGC/blob/master/inst/script/human_mouse.R

[ewels]

Ok, I've added these files and set up the MultiQC config file loading system as described above. This is all documented here.

Thanks for the idea and super fast work @mikelove! Cool new feature :) Let me know if you run into any problems with it.

Phil

[ewels]

Apologies, I'd forgotten about the counts tab. New code written to handle dataseries across multiple dataset plots, counts calculation added. Should now show up in both plots.