Answers.Body: <p>Effect allele is almost always the reference allele. It theoretically does not have to be, so it's good to not rely on that fact. However, in practice it is typically a safe bet. At first glance, the answer to your question could be this:</p>
  <p>The polygenic score from the PGS catalog is using GRCh37. Your genome build is hg38. I expect that some alleles in the score swapped ref/alt positions between builds, resulting in your 205 cases where the effect allele matches the alt allele.</p>
  <p>To confirm, could you list a few examples of cases that matched as expected and a few examples of cases that didn't match?</p>
  <p>You don't need to deal with this encoding if you want to compute a polygenic score on a cohort. You can match based on the allele itself. (Also, there are existing tools which will calculate polygenic scores for you.) Here is an example (just a snippet, you'd have to iterate through your files, etc. for this to work):</p>
  <pre><code>            if PGS_row scoring file identifier like rsID or position == vcf_row identifier: # When the SNP is in both files
              effect_allele = PGS_row[2] # Wherever the affect allele is in your scoring file
              effect_weight = float(PGS_row[4])
              dosage (of ref allele) = (depends on your VCF formatting)
  
              if effect_allele == vcf_row[3]: # Effect allele matches reference allele
                  running_total += (2 - dosage) * effect_weight # We assume dosage of alt is 2 minus dosage of ref. So if you have 2 alt alleles, 2-2 = 0 and there will be nothing added for this SNP.
                  num_additions += 1
              elif effect_allele == vcf_row[4]: # Effect allele matches alternate allele
                  running_total += dosage * effect_weight
                  num_additions += 1
              else:
                  nomatch = nomatch + 1
  
  # Calculate final output score
  final_output_score = running_total / num_additions
  percent_matched = (num_additions)/((nomatch)+(num_additions))
  </code></pre>
  <p>You shouldn't invert the sign of the &quot;effect_weight&quot; column for these 205 cases where the &quot;other_allele&quot; matches the reference. Not having the effect allele doesn't add the opposite beta to the score, it simply doesn't change it.</p>
  <p>Try this:</p>
  <ol>
  <li>Do you still see those 205 cases if you use a genome build matching the PGS file? (I.e., GRCh37.)</li>
  <li>Do you still see those 205 cases if you use the LiftOver PGS scoring file with hg38 coordinates? (Download it <a href="https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002723/ScoringFiles/Harmonized/PGS002723_hmPOS_GRCh38.txt.gz" rel="nofollow noreferrer">here</a>.)</li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-06-15T07:42:36.333Z
Answers.Text: Thanks a lot. This answer is really helpful. Yes, I checked the scoring file of GRCh37 and there were no mismatches, so it must be as you say that the reference base changed between GRCh37 and GRCh38. I also thought inverting the sign was not OK but now you convinced me. Indeed, it add weight if you have the effect allele. Cheers!
Answers.UserDisplayName: Cristian Riccio
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-06-14T19:42:06.980Z
Answers.LastActivityDate: 2023-06-14T20:36:18.490Z
Answers.LastEditDate: 2023-06-14T20:36:18.490Z
Answers.UserDisplayName: BigMistake
Body: <p>How are the &quot;effect_allele&quot; and &quot;other_allele&quot; bases chosen in a PGS file from the PGS Catalog?
  For example, for PGS002723 (<a href="https://www.pgscatalog.org/score/PGS002723/" rel="nofollow noreferrer">https://www.pgscatalog.org/score/PGS002723/</a>), if I set the &quot;effect_allele&quot; to &quot;ref&quot;, then &quot;ref&quot; is equal to the reference
  base of hg38 in 98.60% of the cases on chromosome 22.
  On chromosome 22, there are 205 PGS variants where the &quot;effect_allele&quot; does not match the
  reference base of hg38. In all of these 205 cases, the &quot;other_allele&quot; matches the hg38
  base.</p>
  <p>Why does PGS Catalog encode the hg38 reference base in 14,477 of cases for chromosome 22
  as &quot;effect_allele&quot; and in 205 of cases as &quot;other_allele&quot;?</p>
  <p>Finally, how can I deal with this encoding if I want to compute the polygenic score on a certain cohort?
  Can I simply invert the sign of the &quot;effect_weight&quot; column for these 205 cases where the &quot;other_allele&quot; matches the reference?</p>
  <p>I paste my reproducible code below:</p>
  <pre><code>suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
  suppressMessages(library(data.table))
  suppressMessages(library(R.utils))
  suppressMessages(library(stringr))
  
  setwd('/cluster/home/criccio/PRS/')
  
  source('scripts/get_hg38_base.R')
  
  if (!dir.exists('output/table'))  dir.create('output/table')
  
  PGS_ID &lt;- 'PGS002723'
  
  # Example to ask questions online
  system(paste0('wget -P output/ --recursive --no-parent ftp://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/'))
  
  scoring_file &lt;- data.table::fread(paste0('output/ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/Harmonized/', PGS_ID, '_hmPOS_GRCh38.txt.gz'))
  data.table::setnames(scoring_file, c('chr_name', 'hm_pos', 'effect_allele', 'other_allele'), c('chr', 'pos', 'ref', 'alt'))
  
  pgs_variants_per_chrom &lt;- scoring_file[, .N, chr]
  data.table::setnames(pgs_variants_per_chrom, 'N', 'nb_pgs_in_scoring_file')
  
  chromosome_number &lt;- 22
  
  scoring_file[chr == chromosome_number, ref_hg38 := sapply(scoring_file[chr == chromosome_number, pos], get_hg38_chr22_base)]
  scoring_file[chr == chromosome_number, ref_hg38 := toupper(ref_hg38)]
  scoring_file[chr == chromosome_number, table(ref == ref_hg38)] # 14477 TRUE and 205 FALSE
  scoring_file[chr == chromosome_number, 100 * mean(ref == ref_hg38, na.rm = TRUE)] # 98.60%
  
  scoring_file_mismatch &lt;- scoring_file[chr == chromosome_number &amp; ref != ref_hg38] 
  scoring_file_mismatch[chr == chromosome_number, mean(alt == ref_hg38)] # 100% match
  data.table::fwrite(scoring_file_mismatch, paste0('output/table/', PGS_ID, '_chr_pos_ref_mismatching_hg38_chr', chromosome_number, '.csv'), sep = ',')
  </code></pre>
  <p>get_hg38_base.R</p>
  <pre><code>library(BSgenome.Hsapiens.UCSC.hg38)
  Hsapiens_chr22 &lt;- strsplit(as.character(Hsapiens$chr22), '')[[1]]
  get_hg38_chr22_base &lt;- function(pos) Hsapiens_chr22[pos]
  </code></pre>
  
Tags: <genomics><gwas>
Title: How are the "effect_allele" and "other_allele" bases chosen in a PGS file from the PGS Catalog?
ThreadText: # How are the "effect_allele" and "other_allele" bases chosen in a PGS file from the PGS Catalog?
  
  (Asked by: username_0 on 2023-06-14)
  
  How are the "effect\_allele" and "other\_allele" bases chosen in a PGS file from the PGS Catalog?
  For example, for PGS002723 (<https://www.pgscatalog.org/score/PGS002723/>), if I set the "effect\_allele" to "ref", then "ref" is equal to the reference
  base of hg38 in 98.60% of the cases on chromosome 22.
  On chromosome 22, there are 205 PGS variants where the "effect\_allele" does not match the
  reference base of hg38. In all of these 205 cases, the "other\_allele" matches the hg38
  base.
  
  Why does PGS Catalog encode the hg38 reference base in 14,477 of cases for chromosome 22
  as "effect\_allele" and in 205 of cases as "other\_allele"?
  
  Finally, how can I deal with this encoding if I want to compute the polygenic score on a certain cohort?
  Can I simply invert the sign of the "effect\_weight" column for these 205 cases where the "other\_allele" matches the reference?
  
  I paste my reproducible code below:
  
  ```
  suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
  suppressMessages(library(data.table))
  suppressMessages(library(R.utils))
  suppressMessages(library(stringr))
  
  setwd('/cluster/home/criccio/PRS/')
  
  source('scripts/get_hg38_base.R')
  
  if (!dir.exists('output/table'))  dir.create('output/table')
  
  PGS_ID <- 'PGS002723'
  
  # Example to ask questions online
  system(paste0('wget -P output/ --recursive --no-parent ftp://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/'))
  
  scoring_file <- data.table::fread(paste0('output/ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/', PGS_ID, '/ScoringFiles/Harmonized/', PGS_ID, '_hmPOS_GRCh38.txt.gz'))
  data.table::setnames(scoring_file, c('chr_name', 'hm_pos', 'effect_allele', 'other_allele'), c('chr', 'pos', 'ref', 'alt'))
  
  pgs_variants_per_chrom <- scoring_file[, .N, chr]
  data.table::setnames(pgs_variants_per_chrom, 'N', 'nb_pgs_in_scoring_file')
  
  chromosome_number <- 22
  
  scoring_file[chr == chromosome_number, ref_hg38 := sapply(scoring_file[chr == chromosome_number, pos], get_hg38_chr22_base)]
  scoring_file[chr == chromosome_number, ref_hg38 := toupper(ref_hg38)]
  scoring_file[chr == chromosome_number, table(ref == ref_hg38)] # 14477 TRUE and 205 FALSE
  scoring_file[chr == chromosome_number, 100 * mean(ref == ref_hg38, na.rm = TRUE)] # 98.60%
  
  scoring_file_mismatch <- scoring_file[chr == chromosome_number & ref != ref_hg38] 
  scoring_file_mismatch[chr == chromosome_number, mean(alt == ref_hg38)] # 100% match
  data.table::fwrite(scoring_file_mismatch, paste0('output/table/', PGS_ID, '_chr_pos_ref_mismatching_hg38_chr', chromosome_number, '.csv'), sep = ',')
  
  ```
  
  get\_hg38\_base.R
  
  ```
  library(BSgenome.Hsapiens.UCSC.hg38)
  Hsapiens_chr22 <- strsplit(as.character(Hsapiens$chr22), '')[[1]]
  get_hg38_chr22_base <- function(pos) Hsapiens_chr22[pos]
  
  ```
  
  **Tags:** genomics, gwas
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-06-14. Score: 2 | Accepted answer
  
  Effect allele is almost always the reference allele. It theoretically does not have to be, so it's good to not rely on that fact. However, in practice it is typically a safe bet. At first glance, the answer to your question could be this:
  
  The polygenic score from the PGS catalog is using GRCh37. Your genome build is hg38. I expect that some alleles in the score swapped ref/alt positions between builds, resulting in your 205 cases where the effect allele matches the alt allele.
  
  To confirm, could you list a few examples of cases that matched as expected and a few examples of cases that didn't match?
  
  You don't need to deal with this encoding if you want to compute a polygenic score on a cohort. You can match based on the allele itself. (Also, there are existing tools which will calculate polygenic scores for you.) Here is an example (just a snippet, you'd have to iterate through your files, etc. for this to work):
  
  ```
              if PGS_row scoring file identifier like rsID or position == vcf_row identifier: # When the SNP is in both files
              effect_allele = PGS_row[2] # Wherever the affect allele is in your scoring file
              effect_weight = float(PGS_row[4])
              dosage (of ref allele) = (depends on your VCF formatting)
  
              if effect_allele == vcf_row[3]: # Effect allele matches reference allele
                  running_total += (2 - dosage) * effect_weight # We assume dosage of alt is 2 minus dosage of ref. So if you have 2 alt alleles, 2-2 = 0 and there will be nothing added for this SNP.
                  num_additions += 1
              elif effect_allele == vcf_row[4]: # Effect allele matches alternate allele
                  running_total += dosage * effect_weight
                  num_additions += 1
              else:
                  nomatch = nomatch + 1
  
  # Calculate final output score
  final_output_score = running_total / num_additions
  percent_matched = (num_additions)/((nomatch)+(num_additions))
  
  ```
  
  You shouldn't invert the sign of the "effect\_weight" column for these 205 cases where the "other\_allele" matches the reference. Not having the effect allele doesn't add the opposite beta to the score, it simply doesn't change it.
  
  Try this:
  
  1. Do you still see those 205 cases if you use a genome build matching the PGS file? (I.e., GRCh37.)
  2. Do you still see those 205 cases if you use the LiftOver PGS scoring file with hg38 coordinates? (Download it [here](https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS002723/ScoringFiles/Harmonized/PGS002723_hmPOS_GRCh38.txt.gz).)
  
  > **Comments:**
  > **username_0**: Thanks a lot. This answer is really helpful. Yes, I checked the scoring file of GRCh37 and there were no mismatches, so it must be as you say that the reference base changed between GRCh37 and GRCh38. I also thought inverting the sign was not OK but now you convinced me. Indeed, it add weight if you have the effect allele. Cheers!  
  
  ---
  
  