What should GC coverage bias plot of exome data look like?

Question

I want to see if there is any GC coverage bias in some paired-end Illumina human exome data. Update This is for a variant calling project. I used the CollectMultipleMetrics from GATK4 toolkit with -PROGRAM CollectGcBiasMetrics. An example of the command I used is:

gatk CollectMultipleMetrics -I input.bam -O prefix_dir_name -R input.ref_genome -PROGRAM CollectGcBiasMetrics -PROGRAM CollectAlignmentSummaryMetrics -PROGRAM CollectInsertSizeMetrics -PROGRAM QualityScoreDistribution -PROGRAM MeanQualityByCycle -PROGRAM CollectBaseDistributionByCycle

Update This is my result for the QC coverage bias with one sample highlighted: There seems to be a large bias. I would expect the peak to be around 50% for human. And those weird peaks around 90% are also puzzling to me. On the other hand, I am working with exome data so it might be logical to have a bias towards the higher values of GC% because you have reads on GC-rich exons?

Update I have another plot with the bins produced by the above mentioned tool (.gc_bias.pdf)

Update first part of useful file (gc_bias.detail_metrics) which shows the reads per bin, as suggested by @Maximilian Press. This I will plot.

## htsjdk.samtools.metrics.StringHeader
# CollectMultipleMetrics --INPUT WESP_run/results/aligned/RM8398_MF_S8_MD.sorted.bam --OUTPUT WESP_run/results/aligned/post_alignment_QC/multiple_metrics/RM8398_MF_S8 --PROGRAM CollectAlignmentSummaryMetrics --P
ROGRAM CollectBaseDistributionByCycle --PROGRAM CollectInsertSizeMetrics --PROGRAM MeanQualityByCycle --PROGRAM QualityScoreDistribution --PROGRAM CollectGcBiasMetrics --REFERENCE_SEQUENCE /mnt/scratch_dir/deute
koe/projects/HP/WESPipe/workflow/resources/reference_genomes/hg38/alt_plus/GCA_000001405.15_GRCh38_full_plus_hs38d1_analysis_set.fna --ASSUME_SORTED true --STOP_AFTER 0 --METRIC_ACCUMULATION_LEVEL ALL_READS --IN
CLUDE_UNPAIRED false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_sec
rets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
## htsjdk.samtools.metrics.StringHeader
# Started on: Sat Dec 03 14:25:15 CET 2022

## METRICS CLASS        picard.analysis.GcBiasDetailMetrics
ACCUMULATION_LEVEL      READS_USED      GC      WINDOWS READ_STARTS     MEAN_BASE_QUALITY       NORMALIZED_COVERAGE     ERROR_BAR_WIDTH SAMPLE  LIBRARY READ_GROUP
All Reads       ALL     0       139942  1326    19      0.413544        0.011357                        
All Reads       ALL     1       102412  61      17      0.025996        0.003328                        
All Reads       ALL     2       118957  74      16      0.02715 0.003156                        
All Reads       ALL     3       152181  929     20      0.266429        0.008741                        
All Reads       ALL     4       167732  248     16      0.06453 0.004098                        
All Reads       ALL     5       173319  232     16      0.058421        0.003836                        
All Reads       ALL     6       196365  304     15      0.067567        0.003875                        
All Reads       ALL     7       211217  343     17      0.070875        0.003827                        
All Reads       ALL     8       241125  441     18      0.079822        0.003801                        
All Reads       ALL     9       270272  1068    18      0.172463        0.005277     

I have read a paper that uses the CollectGcBiasMetrics on exome data and they magically produce a plot without bias. I say magically because they don't explain how they use CollectGcBiasMetrics.

Update I added the target coverage plot.

Update Am I missing something? Should I, for instance, use my target regions in instead of the reference genome when using CollectGcBiasMetrics? I will try it, but I am not sure this is the correct thing to do... It is hard to find any workflow for exome data with gatk and particularly CollectGcBiasMetrics, which specifies to use the reference genome.

From the comment @Maximilian Press

Do you have on-instrument diagnostics? Did you run FastQC on the raw reads? What do the other picard metrics look like? Did you manually inspect your alignments (strongly recommended)? Do you actually have coverage of the exome?

I will see if I can get some instrument diagnostics. I did look at the alignments, and as far as I can tell I don't see anything weird? Just that it is very deeply sequenced. I do see some intronic regions, with some reads, but only a few reads compared to the many on the targets, which I think is a consequence of the (too) deep sequencing. The fastqc of the raw reads didn't look troubling, and it did not give me any warnings. The picard metrics I thought did not look out of the ordinary. There is good target coverage. Yes, it is sequenced quite deep. Also, the lab tech did tell me they over-cycled because they wanted to be sure to have all the regions/targets.

Update So this bias might be due to too deep sequencing? And another important question, can I still reliably use this data? If so, I guess I can perform some corrections like suggested here

Any help/suggestions would be very welcome. Thank you in advance.

---

Answer 1

Am I missing something? Should I, for instance, use my target regions instead of the reference genome?

Yes you should do this. Why do you think it's wrong doing that? It should be the first thing to check. And that GC bias might be effected by multiple factors. Most prominent one would be PCR over amplification. Most PCR polymerases give horrible GC bias when library is over cycled. There are many articles published on this.

I also suggest that you look at the number of reads in each bin. Peaks around 90% GC should not puzzle you so much. There is chance that you have only few reads in these bins out of millions of reads.

Answer 2

Just to formalise my comments.

@Supertech's is the right answer.

Background The one thing that is unusual is that the bias peaks at around 75%. The only time I've seen strong GC% bias is for particular groups of bacteria. 75% is about the maximum bias that can be achieved in a protein coding gene and these "extreme bacteria" hit that limit.

I better try and explain that a bit better. The primary selection pressure on a gene is the amino acid composition and the GC% is a second issue, i.e. selection pressure, for some organisms (e.g. certain bacteria, insects ..). Thus amino acids are not mutating, but the triplet codon degeneracy permits a bias in the GC% without alternating a given amino acids. In other words this GC% bias is occurring only for "neutral" mutations (clearly there must be a selection pressure, so they are not really neutral).

Summary and personal thoughts In summary the GC% limit is around 75% for a protein coding gene and beyond that it requires the amino acids to be changed. The peak in @Dandelion's data on that limit is weird because GC% artefacts from Taq errors are not subject to selection pressure: there should not be a peak.

"Contamination" I would avoid as an explanation because in my experience these groups of bacteria are unlikely to invade a lab environment. Sample mix up is a possible explanation. However, I would be cautious about assigning blame because in the chain of command within the wet-lab it could have occurred anywhere.

Official response To reiterate the formal and official answer would be @supertech's response and the 75% GC peak remains a mystery.

---

Please note The second plot has a distribution around GC content around 40% ... There appears a discrepancy between plots. Neither output appears compatible with the GC% of humans, but 40% GC with a reasonable proportion of that data on the 50% mark is better the 75% peak plot. Personally I would investigate the empirical basis for the differences.