1 Objective

It has become a standard to share and make publicly available the summary-level data when publishing a GWAS study. However, there is not a unified format of the summay-level data. To make a gwas summary data into a proper input file for downstream analysis, such as SBayesRC and COJO, we need to format it and check the quality.

In this practice, we will go through all the aspects in formatting, QC and solving missing information. .

2 Data

The example data used here is the summary statistics from A genome-wide association study with 1,126,563 individuals identifies new risk loci for Alzheimer’s disease from Wightman et al., 2021.

The full data are available to download from https://cncr.nl/research/summary_statistics/. For a quicker computation, I only kept chromosome 22 for the practice.

# read in the data
file.name = "chr22_from_PGCALZ2sumstatsExcluding23andMe.txt"
gwas = data.frame(fread(file.name))
head(gwas)
##   chr PosGRCh37 testedAllele otherAllele           z         p     N
## 1  22  16054713            T           C -0.02557039 0.9796000 63926
## 2  22  16055122            T           G  0.84605261 0.3975234  4305
## 3  22  16199921            T           C  0.22625898 0.8210000  1070
## 4  22  16208755            T           C  2.34046003 0.0192600 63926
## 5  22  16281208            T           G  0.04312737 0.9656000 63926
## 6  22  16289869            C           G  0.46420623 0.6425000 63926



We will start a new data frame for the formatted data.

# make a new table in stardard format
n.snp = nrow(gwas)
formatted.gwas = data.frame(
  matrix(NA, nrow = n.snp, ncol = 8)
)
colnames(formatted.gwas) = c("SNP", "A1", "A2", "freq", "b", "se", "p", "N")

3 SNP ID

The commonly used SNP ID is dbSNP IDs, while sometimes there are only physical position of the SNPs and alleles, such as the example data we are using here.

There are a lot of dbSNP resources, such as in NCBI, Resource bundle from GATK, 1000Genome study or HRC imputation panel.

For convenience, I will borrow the snp.info file from eigen-decomposed LD reference to refill the dbSNP IDs. Some of the SNPs in the summary data don’t exist in this reference data, so their information can’t be refilled in this practice, but it doesn’t matter to my further analysis of profiling PRS using the SBayesRC predictors. As well, only chr22 is used in the practice.

# read in a reference data for missing information
freq.file="chr22_snp.info"
ref.freq = data.frame(fread(freq.file))
head(ref.freq)
##   Chrom          ID GenPos  PhysPos A1 A2   A1Freq     N blk
## 1    22   rs4010437      0 16924779  C  T 0.010000 20000 582
## 2    22 rs139656795      0 16926925  C  G 0.013300 20000 582
## 3    22  rs62042126      0 16929771  G  C 0.020775 20000 582
## 4    22   rs4010372      0 16940346  T  G 0.010000 20000 582
## 5    22 rs143291223      0 16959721  T  C 0.010000 20000 582
## 6    22  rs75959072      0 17011656  C  G 0.012725 20000 582

By matching chromosome, bp position and allelele, we will get the dbSNP ID for the SNPs in our predictor, and put it in the formatted data frame.

colname.SNP = "missing"

colname.chr = "chr"
colname.pos = "PosGRCh37"
colname.A1 = "testedAllele"
colname.A2 = "otherAllele"
if(   (colname.SNP != "missing")  ){
  
  formatted.gwas$SNP = gwas[,colname.SNP]
  
}else{
  
  gwas$A = pmin(toupper(gwas[,colname.A1]), toupper(gwas[,colname.A2]))
  gwas$B = pmax(toupper(gwas[,colname.A1]), toupper(gwas[,colname.A2]))
  gwas$chrbpAB = paste0(gwas[,colname.chr], "_", gwas[,colname.pos], "_", gwas$A, "_", gwas$B)
  
  ref.freq$A = with(ref.freq, pmin(A1, A2))
  ref.freq$B = with(ref.freq, pmax(A1, A2))
  ref.freq$chrbpAB = paste0(ref.freq$Chrom, "_", ref.freq$PhysPos, "_", ref.freq$A, "_", ref.freq$B)
  
  gwas$matched_RSID = ref.freq[match(gwas$chrbpAB , ref.freq$chrbpAB), "ID"] 
  gwas[ !(gwas$chrbpAB %in%ref.freq$chrbpAB), "matched_RSID"]  =  gwas[ !(gwas$chrbpAB %in%ref.freq$chrbpAB), "chrbpAB" ]  
  
  formatted.gwas$SNP = gwas$matched_RSID
  
}

head(formatted.gwas[(grep("rs", formatted.gwas$SNP)),])
##             SNP A1 A2 freq  b se  p  N
## 664   rs4010437 NA NA   NA NA NA NA NA
## 672 rs139656795 NA NA   NA NA NA NA NA
## 683  rs62042126 NA NA   NA NA NA NA NA
## 751   rs4010372 NA NA   NA NA NA NA NA
## 840 rs143291223 NA NA   NA NA NA NA NA
## 891  rs75959072 NA NA   NA NA NA NA NA

Both the target data and reference data used here are from human genome build 37. If the gwas data you are working are is using build 38, make sure you also use a reference data from build 38.

4 Alleles

Naming of alleles can be varied. They could be named as “EffectAllele” vs “NonEffectAllele”, or “A1” vs. “A2”, or “ReferenceAllele” vs. “AlternativeAllele”, etc. We need to make it clear which allele is the effect allele, by reading the README files or online instructions for the data.

In cojo format gwas file, A1 is always the effect allele.

If the alleles are in lower case, make sure you make them uppercase in the formattted data.

formatted.gwas$A1= toupper(gwas[,colname.A1])
formatted.gwas$A2 = toupper(gwas[,colname.A2])

5 p value

It’s not common that p value is not reported in the summary data. The p value can get very small value, and they are often in scientific format. You may need to reformat the P-value column if your software recognises it as a character variable.

colname.p = "p"
formatted.gwas$p = as.numeric(gwas[,colname.p])

head(formatted.gwas)
##               SNP A1 A2 freq  b se         p  N
## 1 22_16054713_C_T  T  C   NA NA NA 0.9796000 NA
## 2 22_16055122_G_T  T  G   NA NA NA 0.3975234 NA
## 3 22_16199921_C_T  T  C   NA NA NA 0.8210000 NA
## 4 22_16208755_C_T  T  C   NA NA NA 0.0192600 NA
## 5 22_16281208_G_T  T  G   NA NA NA 0.9656000 NA
## 6 22_16289869_C_G  C  G   NA NA NA 0.6425000 NA

6 Sample Size

Sample size is sometimes missing in the summary data, when we will have to find it out from the publication.

Sometimes disease data would like to include the sample size as two columns, as cases and controls separately. If that’s the case, we will add them up, or calculate effective sample size. Be aware that, it’s not recommended to use effective sample size in SBayesRC.

If it is NCHROBS (number chromosomes observed) reported in summary data, we will divide the number by 2.

colname.N = "N"

sample.size.names=unlist(strsplit(colname.N, ","))

if(length(sample.size.names)==1){
  if(is.na(as.numeric(colname.N) )== T) {
    if(colname.N == "NCHROBS"){
      formatted.gwas$N = 0.5* (gwas[, colname.N])
    }else{
      formatted.gwas$N = gwas[, colname.N]
    }
  }else{
    formatted.gwas$N = as.numeric(colname.N)
  }
}else{
  formatted.gwas$N= gwas[, sample.size.names[1]] + gwas[,sample.size.names[2]]
}

head(formatted.gwas)
##               SNP A1 A2 freq  b se         p     N
## 1 22_16054713_C_T  T  C   NA NA NA 0.9796000 63926
## 2 22_16055122_G_T  T  G   NA NA NA 0.3975234  4305
## 3 22_16199921_C_T  T  C   NA NA NA 0.8210000  1070
## 4 22_16208755_C_T  T  C   NA NA NA 0.0192600 63926
## 5 22_16281208_G_T  T  G   NA NA NA 0.9656000 63926
## 6 22_16289869_C_G  C  G   NA NA NA 0.6425000 63926



Per SNP sample size can be very different when the data is from a study of meta analysis.

hist(formatted.gwas$N, breaks = 200)

We can do a QC by excluding SNPs with too different sample size. Since the distribution is different for each study, it can be arbitrary to give it a cutline of sample size, or a specified number of SD.

In the example data, SD of sample size is very big, and the distribution is skewed. It’s not a good idea to use SD to defined the threshold.

# example code to use SD as threshold. 

a = mean(formatted.gwas$N)
b = sd(formatted.gwas$N)

N.outlier.snps = formatted.gwas[which(formatted.gwas$N > (a+1*b) | formatted.gwas$N < (a-1*b) ),"SNP"]



Another way of QC is to exclude the SNPs with very low sample sizes. At the same time, we don’t want to lose too many SNPs.

a = mean(formatted.gwas$N)
b = sd(formatted.gwas$N)

N.outlier.snps = formatted.gwas[which(formatted.gwas$N < quantile(formatted.gwas$N, 0.1) ),"SNP"]

By excluding the bottom 10% quantile, we will exclude 15080, and the lowest sample size will be 7.4004^{4}.

7 Allele Frequency

ref.freq$gwasA1 = formatted.gwas[match(ref.freq$ID, formatted.gwas$SNP),"A1"]
ref.freq$sign = sign(( as.numeric(ref.freq$gwasA1  == ref.freq$A1))-0.5)
ref.freq$gwasA1freq= abs(as.numeric(ref.freq$gwasA1  != ref.freq$A1) - ref.freq$A1Freq)



Allele frequency is sometimes separately included for cases and controls. If that’s the case, we will do a calculation with the number of cases and controls.

If allele frequency is available, we can compare it to the reference data in a plot, and remove the outliers.

If allele frequency is missing, we can borrow that information from a reference data. This is what happened in our toy data.

colname.freq = "missing"

allele.frequency.names=unlist(strsplit(colname.freq, ","))

if(length(allele.frequency.names) == 2) {
  
  ## Allele frequency is sometimes separately included for cases and controls. 
  freq.case = allele.frequency.names[1]
  freq.control = allele.frequency.names[2]
  N.case = sample.size.names[1]
  N.control = sample.size.names[2]
  
  formatted.gwas$freq = ( (gwas[,freq.case] * gwas[,N.case]) + (gwas[,freq.control] * gwas[,N.control]) ) /(gwas[,N.case] + gwas[,N.control])
  ref.freq$freq.in.gwas = formatted.gwas[match(ref.freq$ID, formatted.gwas$SNP) , "freq"]
  freq.plot = ggplot(data = ref.freq, aes(x = gwasA1freq, y =freq.in.gwas )) + geom_point(size = 0.2) + xlab("AF in reference data") + ylab("AF in GWAS summary statistic")
  ggsave(paste0(file.name, "_AF_plot.png"), freq.plot, height = 8, width = 8)
  
}else{
  
  ## allele frequency is sometimes missing
  if(colname.freq == "missing"){
    formatted.gwas$freq = ref.freq[match( formatted.gwas$SNP, ref.freq$ID),"gwasA1freq"]
    
  }else{
    
    # if freq is available, we will compare it with reference data and make a plot
    formatted.gwas$freq = as.numeric( gwas[,colname.freq] )
    ref.freq$freq.in.gwas = formatted.gwas[match(ref.freq$ID, formatted.gwas$SNP) , "freq"]
    freq.plot = ggplot(data = ref.freq, aes(x = gwasA1freq, y =freq.in.gwas )) + 
      geom_point(size = 0.2) + 
      xlab("AF in reference data") + 
      ylab("AF in GWAS summary statistic")
    ggsave(paste0(file.name, "_AF_plot.png"), freq.plot, height = 8, width = 8)
    
  }
}

head(formatted.gwas[(grep("rs", formatted.gwas$SNP)),])
##             SNP A1 A2     freq  b se         p      N
## 664   rs4010437  C  T 0.010000 NA NA 0.0557100   7409
## 672 rs139656795  C  G 0.013300 NA NA 0.4451000  74004
## 683  rs62042126  G  C 0.020775 NA NA 0.8018050 145339
## 751   rs4010372  T  G 0.010000 NA NA 0.4560000  74004
## 840 rs143291223  T  C 0.010000 NA NA 0.4588000  74004
## 891  rs75959072  C  G 0.012725 NA NA 0.9132982 145339



Here is an example plot from another dataset. You can use a difference threshold 0.2 to exclude outlying SNPs.

This is not the case of our practical data, so I’m just giving an example code here:

freq.outlier = ref.freq[which(abs(ref.freq$gwasA1freq - ref.freq$freq.in.gwas) > 0.2), "ID"]

8 Effect Size

Effect size is usually named as “effect” or “b”. If it’s OR, odds ratio, we will do a log transformation.

If effect size is missing, and only z score is supplied, we will calculate effect size with z score, allele frequency and sample size.

where:
- z is the z-score.
- f is the allele frequency.
- N is the sample size.

colname.b = "missing" 
colname.z = "z"


if(colname.b == "missing" & colname.z != "missing"){
  formatted.gwas$z = gwas[,colname.z]
  formatted.gwas$b = ( formatted.gwas$z )/( ((2 *(formatted.gwas$freq) * (1 - formatted.gwas$freq) * ( formatted.gwas$N + (formatted.gwas$z)^2)))^0.5 )
}  else if(colname.b == "OR" | colname.b == "odds_ratio"){
  formatted.gwas$b = log(gwas[,colname.b])
}else{
  formatted.gwas$b = gwas[,colname.b]
}

head(formatted.gwas[(grep("rs", formatted.gwas$SNP)),])
##             SNP A1 A2     freq            b se         p      N          z
## 664   rs4010437  C  T 0.010000  0.157929336 NA 0.0557100   7409  1.9132973
## 672 rs139656795  C  G 0.013300 -0.017326411 NA 0.4451000  74004 -0.7636095
## 683  rs62042126  G  C 0.020775  0.003264196 NA 0.8018050 145339  0.2510118
## 751   rs4010372  T  G 0.010000 -0.019474057 NA 0.4560000  74004 -0.7454495
## 840 rs143291223  T  C 0.010000 -0.019353227 NA 0.4588000  74004 -0.7408243
## 891  rs75959072  C  G 0.012725  0.001801737 NA 0.9132982 145339  0.1088793

9 Standard error of effect size

If se is missing, we calculate it from effect size, allele frequency and sample size using the formula:

where:
- z is the z-score.
- f is the allele frequency.
- N is the sample size.

colname.se = "missing"

if(colname.se == "missing"){
  
  formatted.gwas$z = sign(formatted.gwas$b) * abs( qnorm(formatted.gwas$p/2) )
  formatted.gwas$se =1/sqrt(2 * formatted.gwas$freq *(1 - formatted.gwas$freq) *(formatted.gwas$N +  (formatted.gwas$z ^ 2)  )  )
  
}else{
  
  formatted.gwas$se = gwas[,colname.se]
  
}

head(formatted.gwas[(grep("rs", formatted.gwas$SNP)),])
##             SNP A1 A2     freq            b         se         p      N
## 664   rs4010437  C  T 0.010000  0.157929336 0.08254302 0.0557100   7409
## 672 rs139656795  C  G 0.013300 -0.017326411 0.02269015 0.4451000  74004
## 683  rs62042126  G  C 0.020775  0.003264196 0.01300415 0.8018050 145339
## 751   rs4010372  T  G 0.010000 -0.019474057 0.02612391 0.4560000  74004
## 840 rs143291223  T  C 0.010000 -0.019353227 0.02612391 0.4588000  74004
## 891  rs75959072  C  G 0.012725  0.001801737 0.01654802 0.9132982 145339
##              z
## 664  1.9132973
## 672 -0.7636095
## 683  0.2510118
## 751 -0.7454495
## 840 -0.7408243
## 891  0.1088793

10 output into cojo format

We can exclude the SNPs with missing information, outlying sample size or distinct allele frequency compared to the LD reference dataset, before write it into a cojo file.

output = paste0(file.name, ".ma")

formatted.gwas = formatted.gwas[which(is.na(formatted.gwas$freq) == F) ,]

formatted.gwas = formatted.gwas[which(!formatted.gwas$SNP %in% N.outlier.snps) ,]

formatted.gwas = formatted.gwas[,1:8]

write.table(formatted.gwas, file=output, quote = F, sep ="\t", row.names = F)

This is the first few lines of the formatted and QCed data:

head(formatted.gwas)
##             SNP A1 A2     freq            b          se         p      N
## 672 rs139656795  C  G 0.013300 -0.017326411 0.022690147 0.4451000  74004
## 683  rs62042126  G  C 0.020775  0.003264196 0.013004153 0.8018050 145339
## 751   rs4010372  T  G 0.010000 -0.019474057 0.026123910 0.4560000  74004
## 840 rs143291223  T  C 0.010000 -0.019353227 0.026123911 0.4588000  74004
## 891  rs75959072  C  G 0.012725  0.001801737 0.016548016 0.9132982 145339
## 952 rs145889691  C  A 0.082775  0.005568544 0.006731399 0.4080959 145339