## Understanding Genetic Data Structure

(This lab was adapted for statsTeachR by Sara Nunez, Nicholas Reich and Andrea Foulkes from our GWAS tutorial.)

The first step in any statistical analysis is to understand the data your working with. No analysis is useful if the data formats and meanings are misunderstood. Genetic data is complex in that it is not just a simple .csv file that can be read into R. Instead, data are typically organized into either (1) .ped and .map files or (2) .bim, .bed and .fam files. R can be used to read in either set of file types (utilizing external packages), though the second set of files is substantially smaller and thus generally preferable*. In this lab, we describe and explore the elements of the .bim, .bed and .fam files, as well as show how these files can be read into R for analysis.

### PennCath data

In this module series, we utilize the PennCath cohort data, which arises from a Genome-wide association (GWA) study of coronary artery disease (CAD) and cardiovascular risk factors based at the University of Pennsylvania Medical Center. For computational purposes, we will be working with a random subset of the .bim, .bed and fam files, which have been made publically available at https://www.mtholyoke.edu/courses/afoulkes/Data/statsTeachR/.

The complete data set has a total of n = 3850 individuals that were recruited between July 1998 and March 2003. A nested case-control study of European ancestry severe angiographic CAD cases and angiographic normal controls were selected for genome wide genotyping. Samples were genotyped using the Affymetrix 6.0 GeneChip and genotypes were called using the Birdseed calling algorithm. De-identified data used in the tutorial are comprised of n = 1401 individuals with genotype information across 861,473 single-nucleotide polymorphisms (SNPs). Corresponding clinical data, including age, sex, high-density lipoprotein (HDL)-cholesterol (HDL-C), low-density lipoprotein (LDL)-cholesterol, triglycerides (TGs) and CAD status are available as well. HDL-cholesterol, LDL- cholesterol and TGs are all quantitative traits and well- described cardiovascular disease risk factors. Notably, PennCath is one of the core GWA studies nested within the Coronary ARtery DIsease Genome-wide Replication And Meta-analysis (CARDIoGRAM) consortium meta-data, and serves as a representative regional population with no admixture.

source("http://bioconductor.org/biocLite.R")
biocLite("snpStats")
##
##  /var/folders/q5/brz1dtz54w99g1mft67shggm0000gp/T//Rtmp2RngR5/downloaded_packages
library(snpStats)
library(RCurl)

# read in data in two steps:

destfile = "GWAS_data.bed")
destfile = "GWAS_data.fam")

# what class is "data"
class(data)
## [1] "list"
# how many elements are in "data"
length(data)
## [1] 3

You will first notice that when read in this way, the data is of class “list” and has 3 elements in it. Each element corresponds to the three kinds of data that were read into R. We will next explore each component as each serves a different purpose in analysis and will be needed in the upcoming labs.

#### Exercise

Before we look at each piece of data in detail, try looking at the general structure of it using the function str(). You’ll notice that there is information on all three of the objects in the list. What does each element of the list correspond to (e.g. what is it called and what appears to be in it)? Of what class" is the genotypes object?

### The .bim file

The .bim file contains information on each SNP and can be seen by calling the map element in the data object (see tutorial for relationship between .bim and .map file formats).

### The .bed file

The raw version of the .bed file contains a binary version of the genotype data. This is the largest of the three files because it contains every SNP in the study, as well as the genotype at that snp for each individual. When we read this file into R, however, it’s automatically translated into the genotype information from the binary code. R then store this object as a SnpMatrix object that is critical to SNP analysis. This SnpMatrix corresponds to the genotypes element in the data object.

(bed <- data$genotypes) ## A SnpMatrix with 1401 rows and 500000 columns ## Row names: 10002 ... 11596 ## Col names: rs4579145 ... rs946221 #### Exercise What do the columns in this matrix represent? What do the rows represent? How many SNPs are in this matrix? How many individuals are present? What does the element corresponding to the $$ith$$ row and $$jth$$ column represent? ### The .fam file Finally, the .fam file contains the participant identification information, including a row for each individual and six columns, corresponding to “Family ID Number”, “Sample ID Number”,“Paternal ID Number”, “Maternal ID Number”, “Sex”, and “Phenotype”. The .fam file can be called by referring to the .fam element of the data object. fam <- data$fam
head(fam)
##       pedigree member father mother sex affected
## 10002    10002      1      0      0   1        1
## 10004    10004      1      0      0   2        1
## 10005    10005      1      0      0   1        2
## 10007    10007      1      0      0   1        1
## 10008    10008      1      0      0   1        2
## 10009    10009      1      0      0   1        2

Note that not all of these columns contain information depending the nature of the data collection process and study design. In this case, we also have a separate clinical file with several outcomes and covariates that will be read into R in the next section.

#### Exercise

How many observations are there in the fam file? How many males are there(assuming males are coded as 1)? How many females? How many levels of the phenotype variable are there? What does the pedigree" column correspond to? Does it look familar (where else have you seen it)?

### The clinical data

Finally, in some instances you may have acess to a separate phenotype file that contains information on things like additional covariates for each individual and other “outcome variables” you may be interested in (e.g. something other than affected or variable that codes affected differently). The rows of this file represent each subject and the columns correspond to the available covariates and phenotypes. We read this data in separately as it is a simple .txt file.

library(RCurl)
clinicalURL <- getURL("https://www.mtholyoke.edu/courses/afoulkes/Data/statsTeachR/GWAS_clinical.csv")
clinical <- read.csv(text = clinicalURL, colClasses = c("character", rep("factor",
2), rep("numeric", 2)))
rownames(clinical) <- clinical$FamID print(head(clinical)) ## FamID CAD sex age tg hdl ldl ## 10002 10002 1 1 60 NA <NA> <NA> ## 10004 10004 1 2 50 55 23 75 ## 10005 10005 1 1 55 105 37 69 ## 10007 10007 1 1 52 314 54 108 ## 10008 10008 1 1 58 161 40 94 ## 10009 10009 1 1 59 171 46 92 #### Exercise Explore the clinical data. How many complete observations are there? What is the age range of study participants? Which column of this data set links it to all the others? How many people are CAD$$+$$ ($$CAD=1$$)? ### Saving work for following labs… Now that all the data is read into R, we need to save the objects for the rest of the labs. Note that some of the objects must first be renamed for the purpose of future labs. You can rename and save the necessary data by running the following commands: genotype <- data$genotypes
genoBim <- data\$map
colnames(genoBim)<-c("chr", "SNP", "gen.dist", "position", "A1", "A2")

save(genotype,genoBim,clinical,file="lab1_save.RData")