Upload New File - R code for single gene PTV-burden association analysis with lipids in the UKBB

UKBB_analysis/PTV_burden.associations.R

+#######################################################
+######### UKBB lipids PTV-burden associations #########
+######### Author: Yunfeng Huang               #########
+#######################################################
+
+# Read phenotype data
+
+library(data.table)
+lipids_pheno <- as.data.frame(fread("~/lipids_heiko/LoF_UKB/2021/lipids_lofs_pheno.txt.gz"))
+
+# single gene burden
+
+library(MASS)
+
+res <- as.data.frame(matrix(NA,120,5))
+colnames(res) <- c("Trait","Gene","Beta","P-value","N_ptv_carriers")
+traits <- c("HDL","LDL","TC","TG")
+genes <- colnames(lipids_pheno)[8:37]
+n <- 0
+for(i in 1:4){
+  for(j in 1:30){
+    n <- n+1
+    sub <- na.omit(dplyr::select(lipids_pheno,paste(tolower(traits[i]),"_res",sep=""),genes[j]))
+    if(ncol(sub)!=2){next}
+    colnames(sub) <- c("resid","gene")
+    n_carrier <- sum(sub$gene>0)
+    g <- rlm(resid~gene,data=sub)
+    res$Trait[n] <- traits[i]
+    res$Gene[n] <- genes[j]
+    res$Beta[n] <- summary(g)$coeff[2,1]
+    res$`P-value`[n] <- pnorm(abs(summary(g)$coeff[2,3]),lower.tail=F)*2
+    res$N_ptv_carriers[n] <- n_carrier
+  }
+  print(i)
+}