SDS script

[ABThijssen]

Hi! Thanks for this great paper. I have a couple of other traits that were not included in the paper that I want to check for evolutionary signatures.
I ran into an issue trying to get the tSDS scores with your code. Specifically, this part of the code seems to assume that there are sds columns in the gwas right? Did you merge that somewhere outside of this script?

sds=fread("/lustre/home/acct-bmelgn/bmelgn-3/ldsc/SDS.txt",data.table=F)
sds_rho=c()
setwd("/lustre/home/acct-bmelgn/bmelgn-3/ldsc/all_gwas/")
write.csv(c("name","rho_true","p"),"sds_rho.csv")
tmp=list.files(pattern="*.gz")
for (name in tmp){
gwas=fread(name,data.table=FALSE)
gwas=gwas[which(gwas$sds!=0),]
x=match(gwas$SNP,sds$ID)
gwas$freq=sds[x,"DAF"]

[WeiCSong]

Thanks for your message. Yes, in my analysis i first add the SDS column to the gwas summary data. I downloaded SDS from the original paper, mapped to summary data by rsid, and check if the derived allele (recorded in SDS data) is trait-increasing (if not, SDS sign was flipped). Note that GWAS A1 is not always the derived allele.

[ABThijssen]

Thank you so much for your reply. I added the SDS column and flipped the SDS sign to match trait-increasing alleles. As a test I work with the schizophrenia gwas, to try to reproduce the results you have in the supp table 3. However, I get a pretty different value for rho_true (0.05 instead of -0.35). Below is the complete code, any idea what is going wrong?

library(data.table)
library(tidyverse)

sds_rho=c()

#####Read in SDS and GWAS
sds <- fread("SDS/SDS_UK10K_n3195_release_Sep_19_2016.tab",data.table=F)
gwas <- fread("data/sumstats/PGC3_SCZ_wave3.european.autosome.public.v3.vcf.tsv.gz",data.table=F)

######Add sds column to the GWAS, check if derived allele is the trait increasing allele, otherwise flip SDS sign
gwas <- left_join(gwas, sds, by = "ID")
gwas <- gwas %>% mutate(sds = ifelse(DA == A1, SDS, -SDS)) %>% select(-c("SDS"))

gwas=gwas[which(gwas$sds!=0),]
x=match(gwas$ID,sds$ID)
gwas$freq=sds[x,"DAF"]

######re-normalized per DAF bin
gwas$SDS=0
for(i in seq(0.05,0.94,0.01)){
x=gwas[which(gwas$freq>i&gwas$freq<i+0.01),"sds"]
x=(x-mean(x))/sd(x)
gwas[which(gwas$freq>i&gwas$freq<i+0.01),"SDS"]=x
}

########spearman rho
gwas=gwas[order(gwas$PVAL),]
res=c()
for(i in seq(1,nrow(gwas),1000)){
start=i
end=min(i+1000,nrow(gwas))
s=median(gwas[start:end,"SDS"])
res=rbind(res,c(i,s))
}
rho_true=cor(res[,1],res[,2],method="spearman")

[WeiCSong]

check if the derived allele (recorded in SDS data) is trait-increasing (if not, SDS sign was flipped)

despite DA==A1, you also need to check the direction of beta. tSDS is the SDS corresponded to trait-increasing allele, and tSDS>0 means that the allele that literally "increase" the trait was favored over the past 2000 years. If DA==A1 and beta<0, you still need to flip the direction.

[ABThijssen]

Thanks so much, I understand now what you meant. Unfortunately, the rho that I get is still way off (it is now 0.26 instead of -0.35). This is the current code:

sds <- fread("SDS/SDS_UK10K_n3195_release_Sep_19_2016.tab",data.table=F)
gwas <- fread("data/sumstats/PGC3_SCZ_wave3.european.autosome.public.v3.vcf.tsv.gz",data.table=F)

####Add sds column to the GWAS, and flip sign if necessary so it will reflect the trait-increasing allele
gwas <- left_join(gwas, sds, by = "ID")
gwas <- gwas %>%
mutate(sds = case_when(
DA == A1 & BETA > 0 ~ SDS,
DA == A1 & BETA < 0 ~ -SDS,
DA == A2 & BETA > 0 ~ -SDS,
DA == A2 & BETA < 0 ~ SDS)) %>%
select(-c("SDS"))

gwas=gwas[which(gwas$sds!=0),]
x=match(gwas$ID,sds$ID)
gwas$freq=sds[x,"DAF"]

######re-normalized per DAF bin
gwas$SDS=0
for(i in seq(0.05,0.94,0.01)){
x=gwas[which(gwas$freq>i&gwas$freq<i+0.01),"sds"]
x=(x-mean(x))/sd(x)
gwas[which(gwas$freq>i&gwas$freq<i+0.01),"SDS"]=x
}

########spearman rho
gwas=gwas[order(gwas$PVAL),]
res=c()
for(i in seq(1,nrow(gwas),1000)){
start=i
end=min(i+1000,nrow(gwas))
s=median(gwas[start:end,"SDS"])
res=rbind(res,c(i,s))
}
rho_true=cor(res[,1],res[,2],method="spearman")

[WeiCSong]

could you share the merged gwas summary statistics so that i could compare it with the one i generated?

[WeiCSong]

Ah, my bad, sorry for my carelessness. This version of code sorted the p value in the ascending order, so larger bin id actually means less significance. Then, positive rho -> larger the SDS larger the p-value -> significant SNPs has small SDS-> negative adaptation. On the contrary, if we want negative adaptation as negative rho, p-value should be in descending order, like figure 4A of SDS paper and figure 3c of my paper.

[ABThijssen]

Thanks so much, that seems to work!