Clumping

I run the plink clumping algorithm on this GWAS summary statistics to have an insight of the unique association loci. Here are the clumping options I’ve used. Note: GWAS summary file must NOT be zipped.

• p1 is the Significance threshold for an association
• p2 is Secondary significance threshold for clumped SNPs = SNPs that are part of the association, but not pic.
• r2 is the LD threshold for clumping and
• clump-kb is the physical distance threshold for clumping

# good folder
cd /home/y.holtz/VITAMIND_XIA_ET_AL/2_CLUMPING

# Sum stat of each chromosome
for i in $(seq 1 22) ; do
  cat /home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS/GWAS_vitaminD_XiaEtAL.linear | head -1  > file_chr_${i}
  cat /home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS/GWAS_vitaminD_XiaEtAL.linear | grep "^${i} " >> file_chr_${i}
done

# Prepare command
tmp_command="plink \
          --bfile /gpfs/gpfs01/polaris/Q0286/UKBiobank/v2EURu_HM3/ukbEURu_imp_chr{TASK_ID}_v2_HM3_QC\
          --clump file_chr_{TASK_ID} \
          --clump-p1 0.00000001 --clump-p2 0.01 --clump-r2 0 --clump-kb 1000 \
          --out /home/y.holtz/VITAMIND_XIA_ET_AL/2_CLUMPING/GWAS_vitaminD_XiaEtAL_.clump_chr{TASK_ID}"
          
# Send qsub
qsubshcom "$tmp_command" 1 50G clumping_plink 30:00:00 "-array=1-22"

It is more handy to have only one file with all the chromosome. Let’s build it:

# good folder
cd ~/VITAMIND_XIA_ET_AL/2_CLUMPING

# concatenate
cat *clumped | grep "CHR" | head -1 > result_GWAS_VitaminD_XiaEtAl_clumped
cat GWAS*clumped | grep -v "CHR" | sed '/^$/d' | grep -v "NA"  >> result_GWAS_VitaminD_XiaEtAl_clumped

# Transfert it locally
cd /Users/y.holtz/Dropbox/QBI/4_UK_BIOBANK_GWAS_PROJECT/VitaminD-GWAS/0_DATA
scp  y.holtz@delta.imb.uq.edu.au:/home/y.holtz/VITAMIND_XIA_ET_AL/2_CLUMPING/result_GWAS_VitaminD_XiaEtAl_clumped  .

Do I have my 6 loci after clumping?

clump=read.table("0_DATA/result_GWAS_VitaminD_XiaEtAl_clumped", header=T)
clump

FUMA

FUMA is a tool for functional mapping and annotation of genetic associations ref. A FUMA analysis of this GWAS sumstat has been run. The result is accessible here.

• Id: y.holtz@uq.edu.au
• Psw: fumayanholtz

# transfer locally for FUMA?
cd /Users/y.holtz/Desktop
scp y.holtz@delta.imb.uq.edu.au:/home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS/GWAS_vitaminD_XiaEtAL.linear.gz .

A few thoughts on the FUMA analysis:

• It looks like a lot of SNP are missing in their manhattan plot version
• Only one Genomic risk loci is reported. Chromosome 4, on SNP rs113227722 at position 71316300. It’s weird it is not the best SNP of this area. It has a very low MAF: 0.04.
  
• The gene-based test outputs only one significant genes: MUC7 in the chromosome 4 area. This is not what I excpected, not what is described in Xia et al.

–> These results are very weird, I need to investigate further why FUMA does not work properly.

About genes position

• Xia et al. locate the SNP rs3755967 on chromosome 4, position 72,828,262 on their nature paper. I’ve read their publication several times, including supplementary data, but I didn’t see any information concerning the HG version used. After my mail they revised their data and produced a new file, with position 72,609,398

• In my UKB data, this SNP rs3755967 is located at position 72,609,398. As far as I know UKB uses the Hg19 version of Human genome, which is the same as Hg37 if I understood well.

• The dbSNP website tells that rs3755967 position is 71,743,681. It also says it is in the GC gene. It uses the GRCh38 version if I understood well.

• The NCBI platform gives information on the gene GC: - Its position is 71,741,693 to 71,805,520 in HG 3

- Conclusion

Note: cleaning the GWAS sumstat:

# I download the vitamin-D summary statistics from here: https://drive.google.com/drive/folders/0BzYDtCo_doHJRFRKR0ltZHZWZjQ
# Then from local download folder to delta
gzip upload_25HydroxyVitaminD_QC.METAANALYSIS1.txt
scp upload_25HydroxyVitaminD_QC.METAANALYSIS1.txt.gz  y.holtz@delta.imb.uq.edu.au:/home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS

# good to delta folder
cd /home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS

# I miss the frequence of the major allele in this file. I can use the hapmap3 frequencies previously calculated
cat /home/y.holtz/BLOOD_GWAS/1_GWAS/PLINK/all_frequency_1.frq

# --- MA FORMAT --- make the file with this new column, transform it in a format suitable for GSMR.
echo "SNP A1 A2 freq b se p n" > GWAS_vitaminD_XiaEtAL.ma
join <(zcat upload_25HydroxyVitaminD_QC.METAANALYSIS1.txt.gz | sort -k 1,1) <(awk '{print $2, $6}' /home/y.holtz/BLOOD_GWAS/1_GWAS/PLINK/all_frequency*.frq | sort -k 1,1) | awk '{ print $1,$2,$3,$15,$4,$5,$6,$14}' >>  GWAS_vitaminD_XiaEtAL.ma

# --- LINEAR FORMAT --- make a .linear format suitable for PLINK or FUMA
echo "CHR SNP BP A1 TEST NMISS BETA STAT P" > GWAS_vitaminD_XiaEtAL.linear
zcat upload_25HydroxyVitaminD_QC.METAANALYSIS1.txt.gz | grep -v 'Chr' | awk '{ print $12, $1, $13, $2, "ADD", $14, $4, $4*$5, $6}' >> GWAS_vitaminD_XiaEtAL.linear
R  # We need to remove the 10^5 format for 100000 format
data=read.table("GWAS_vitaminD_XiaEtAL.linear", header=T)
data$BP = round(data$BP, 0)
options(scipen=500)
library(dplyr)
data <- data %>% arrange(CHR, BP)
write.table(data, file="GWAS_vitaminD_XiaEtAL.linear", quote=F, row.names=FALSE )

# --- LIGHT FORMAT --- make a light version for github
cat GWAS_vitaminD_XiaEtAL.linear | head -1 >  GWAS_vitaminD_XiaEtAL.linear.light
cat GWAS_vitaminD_XiaEtAL.linear | awk '{if( $9 < 0.1){print $0}}' >>  GWAS_vitaminD_XiaEtAL.linear.light
gzip GWAS_vitaminD_XiaEtAL.linear.light

# Transfer locally?
cd ~/Dropbox/QBI/4_UK_BIOBANK_GWAS_PROJECT/VitaminD-GWAS/0_DATA
scp y.holtz@delta.imb.uq.edu.au:/home/y.holtz/VITAMIND_XIA_ET_AL/1_GWAS/GWAS_vitaminD_XiaEtAL.linear.gz .

# Note: I need to delete 2 SNPs that have a P of 0 in the R output. Because R rounds them to 0
more GWAS_vitaminD_XiaEtAL.linear | awk '{ if($9!=0){print $0}}' | gzip > GWAS_vitaminD_XiaEtAL.linear_IGV.gz