Answers.Body: <blockquote>
  <p>Am I missing something? Should I, for instance, use my target regions
  instead of the reference genome?</p>
  </blockquote>
  <p>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.</p>
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-30T03:40:19.663Z
Answers.LastActivityDate: 2022-12-30T03:40:19.663Z
Answers.UserDisplayName: Supertech
Answers.Body: <p>Just to formalise my comments.</p>
  <p>@Supertech's is the right answer.</p>
  <p><strong>Background</strong> 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 <em>protein coding gene</em> and these &quot;extreme bacteria&quot; hit that limit.</p>
  <p>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 &quot;neutral&quot; mutations (clearly there must be a selection pressure, so they are not really neutral).</p>
  <p><strong>Summary and personal thoughts</strong> 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 <em>Taq</em> errors are not subject to selection pressure: there should not be a peak.</p>
  <p>&quot;Contamination&quot; 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.</p>
  <p><strong>Official response</strong> To reiterate the formal and official answer would be @supertech's response and the 75% GC peak remains a mystery.</p>
  <hr />
  <p><strong>Please note</strong> 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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-30T18:35:40.143Z
Answers.LastActivityDate: 2022-12-30T22:50:09.297Z
Answers.LastEditDate: 2022-12-30T22:50:09.297Z
Answers.UserDisplayName: M__
Body: <p>I want to see if there is any GC coverage bias in some paired-end Illumina human exome data. <strong>Update</strong> This is for a variant calling project. I used the <code>CollectMultipleMetrics</code> from <a href="https://gatk.broadinstitute.org/hc/en-us/articles/360037594031-CollectMultipleMetrics-Picard-" rel="nofollow noreferrer">GATK4 toolkit</a> with <code>-PROGRAM</code> <a href="https://gatk.broadinstitute.org/hc/en-us/articles/360036801531-CollectGcBiasMetrics-Picard-" rel="nofollow noreferrer">CollectGcBiasMetrics</a>. An example of the command I used is:</p>
  <pre><code>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
  </code></pre>
  <p><strong>Update</strong> This is my result for the QC coverage bias with one sample highlighted:
  <a href="https://i.sstatic.net/WgnZn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/WgnZn.png" alt="GC coverage bias plot" /></a>
  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?</p>
  <p><strong>Update</strong> I have another plot with the bins produced by the above mentioned tool (.gc_bias.pdf)
  <a href="https://i.sstatic.net/tVCcb.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/tVCcb.png" alt="enter image description here" /></a></p>
  <p><strong>Update</strong> first part of useful file (gc_bias.detail_metrics) which shows the reads per bin, as suggested by @Maximilian Press. This I will plot.</p>
  <pre><code>## 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     
  </code></pre>
  <p>I have read a <a href="https://hereditasjournal.biomedcentral.com/articles/10.1186/s41065-021-00171-3#Sec14" rel="nofollow noreferrer">paper</a> 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.</p>
  <p><strong>Update</strong> I added the target coverage plot.
  <a href="https://i.sstatic.net/UNqTb.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/UNqTb.png" alt="Target region coverage" /></a></p>
  <p><strong>Update</strong> Am I missing something? Should I, for instance, use my target regions in instead of the reference genome when using <code>CollectGcBiasMetrics</code>? 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 <code>CollectGcBiasMetrics</code>, which specifies to use the reference genome.</p>
  <p>From the comment @Maximilian Press</p>
  <blockquote>
  <p>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?</p>
  </blockquote>
  <p>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.</p>
  <p><strong>Update</strong> 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 <a href="https://support.illumina.com/content/dam/illumina-support/help/Illumina_DRAGEN_Bio_IT_Platform_v3_7_1000000141465/Content/SW/Informatics/Dragen/BiasCorrections_fDG_dtSW.htm" rel="nofollow noreferrer">here</a></p>
  <p>Any help/suggestions would be very welcome. Thank you in advance.</p>
  <hr />
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-30T14:42:51.737Z
Comments.Text: A bit of upvoting would be cool, its very good QC indeed and clearly presented. @supertech's answer makes sense to me. The GC% bias in living organisms could only represent bacteria IMO, I've seen this in soil bacteria so this has to be a lab artefact.
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-30T17:14:25.567Z
Comments.Text: Thank you @M__. About the bacterial IMO (which I think means a contamination?). There is a separate run (another day a protocol slightly different) that shows the same bias. So either the lab tech is very sloppy, or the protocol was not correct (too many PCR cycles or something).
Comments.UserDisplayName: Dandelion
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-12-31T00:51:46.017Z
Comments.Text: Thank you for incorporating the comments (and your answers to the comments) into your question. This is how comments should be used, and makes my job as a moderator much easier.
Comments.UserDisplayName: gringer
Tags: <gatk><illumina><human-genome><coverage><exome>
Title: What should GC coverage bias plot of exome data look like?
ThreadText: # What should GC coverage bias plot of exome data look like?
  
  (Asked by: username_0 on 2022-12-29)
  
  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](https://gatk.broadinstitute.org/hc/en-us/articles/360037594031-CollectMultipleMetrics-Picard-) with `-PROGRAM` [CollectGcBiasMetrics](https://gatk.broadinstitute.org/hc/en-us/articles/360036801531-CollectGcBiasMetrics-Picard-). 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:
  [![GC coverage bias plot](https://i.sstatic.net/WgnZn.png)](https://i.sstatic.net/WgnZn.png)
  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)
  [![enter image description here](https://i.sstatic.net/tVCcb.png)](https://i.sstatic.net/tVCcb.png)
  
  **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](https://hereditasjournal.biomedcentral.com/articles/10.1186/s41065-021-00171-3#Sec14) 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.
  [![Target region coverage](https://i.sstatic.net/UNqTb.png)](https://i.sstatic.net/UNqTb.png)
  
  **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](https://support.illumina.com/content/dam/illumina-support/help/Illumina_DRAGEN_Bio_IT_Platform_v3_7_1000000141465/Content/SW/Informatics/Dragen/BiasCorrections_fDG_dtSW.htm)
  
  Any help/suggestions would be very welcome. Thank you in advance.
  
  ---
  
  **Tags:** gatk, illumina, human-genome, coverage, exome
  
  > **Comments:**
  > **username_1**: A bit of upvoting would be cool, its very good QC indeed and clearly presented. @username_3's answer makes sense to me. The GC% bias in living organisms could only represent bacteria IMO, I've seen this in soil bacteria so this has to be a lab artefact.  
  > **username_0**: Thank you @username_1. About the bacterial IMO (which I think means a contamination?). There is a separate run (another day a protocol slightly different) that shows the same bias. So either the lab tech is very sloppy, or the protocol was not correct (too many PCR cycles or something).  
  > **username_2**: Thank you for incorporating the comments (and your answers to the comments) into your question. This is how comments should be used, and makes my job as a moderator much easier.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2022-12-30. Score: 2 | Accepted answer
  
  > 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 by username_1 on 2022-12-30. Score: 1
  
  Just to formalise my comments.
  
  @username_3'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 @username_0'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 @username_3'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.
  
  
  ---
  
  