Regression output tables using lactate SNPs on lactate levels in ProtecT as an example

Use MR Base to find the SNPs linked to lactate

• rs762523 on ch22
• rs1260326 on ch02

---

Move the files to a folder in my home directory using WinSCP and use gunzip to unzip ch22 and ch02

Pathway to folder: /panfs/panasas01/sscm/ca16591/ProtecT_genetic_Vanessa

---

Use gcta to convert the (MaCH) genotype data to Plink format

gcta64 –dosage-mach-gz data_chr22.step2.mldose.gz data_chr22.step2.mlinfo.gz –make-bed –out test

---

Write a bash script to call Plink

I called the script lactate_snp_extract.sh

#!/bin/bash
#PBS -l walltime=00:05:00,nodes=1:ppn=1
#PBS -o snps.txt
#PBS -j oe
# Set the name of the job
#PBS -N lactate_snp22_extract
echo Running on host `hostname`
echo Time is `date`
echo Directory is `pwd`
echo PBS job ID is $PBS_JOBID
echo This jobs runs on the following machines:
echo `cat $PBS_NODEFILE | uniq`

# First change the MaCH files into Plink format with gcta. First get all in the zip format to use the --dosage-mach-gz option
# gcta64 --dosage-mach-gz data_chr22.step2.mldose.gz data_chr22.step2.mlinfo.gz --make-bed --out test

module add apps/plink-1.90
plink --bfile /panfs/panasas01/sscm/ca16591/ProtecT_genetic_Vanessa/test --snp rs762523 --recode A --out /panfs/panasas01/sscm/ca16591/ProtecT_genetic_Vanessa/Protect_lactate_22snps 

---

Import the results for ch22 to R and select rs762523_A

---

Merge the genetic data for rs762523_A with the ProtecT metabolite data

---

Recode the rs762523_A such that the G is the risk allele

---

Tabulation of cases and controls by rs762523_A

	0	1	2
Case	132	406	348
Control	77	283	263

---

Run a simple linear regression of genotype (rs762523) on log-transformed lactate

rs762523_A_on_Lac = lm(Lac ~ rs762523_A, data = rs762523_A)
summary(rs762523_A_on_Lac)

## 
## Call:
## lm(formula = Lac ~ rs762523_A, data = rs762523_A)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95446 -0.22900 -0.05605  0.17879  1.92147 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.54901    0.02403  22.848   <2e-16 ***
## rs762523_A1  0.03767    0.02743   1.373    0.170    
## rs762523_A2  0.02273    0.02784   0.816    0.414    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3474 on 1506 degrees of freedom
## Multiple R-squared:  0.001327,   Adjusted R-squared:  4.543e-07 
## F-statistic:     1 on 2 and 1506 DF,  p-value: 0.368

---

A basic table for regression output

library(xtable)
print(xtable(summary(rs762523_A_on_Lac)), type='html')

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	0.5490	0.0240	22.85	0.0000
rs762523_A1	0.0377	0.0274	1.37	0.1699
rs762523_A2	0.0227	0.0278	0.82	0.4143

---

Another way to make a table of regression output

library(sjPlot)

## Visit http://strengejacke.de/sjPlot for package-vignettes.

library(htmltools)
library(sjmisc)
library(magrittr)

table =lm(Lac ~ rs762523_A, data = rs762523_A) %>%
sjt.lm(no.output=TRUE, show.se = TRUE) %>% 
  return() %>% .[["knitr"]]  %>%  asis_output

		Lac
		B	CI	std. Error	p
(Intercept)		0.55	0.50 – 0.60	0.02	<.001
rs762523_A
1		0.04	-0.02 – 0.09	0.03	.170
2		0.02	-0.03 – 0.08	0.03	.414
Observations		1509
R2 / adj. R2		.001 / .000

---

---