This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision Next revision Both sides next revision | ||
keller_and_evans_lab:gscan_db_ga_p [2017/02/13 21:35] scott /* Phenotypes */ |
keller_and_evans_lab:gscan_db_ga_p [2017/02/15 13:51] scott /* ARIC */ |
||
---|---|---|---|
Line 576: | Line 576: | ||
====== ARIC ====== | ====== ARIC ====== | ||
- | |||
- | (Hannah/ | ||
- | |||
- | ===== ID Mapping ===== | ||
Line 883: | Line 879: | ||
====== eMERGE ====== | ====== eMERGE ====== | ||
- | (Hannah/ | ||
===== Phenotypes ===== | ===== Phenotypes ===== | ||
- | Description of phenotypes | + | |
+ | |||
+ | options(stringsAsFactors=F) | ||
+ | |||
+ | |||
+ | ### eMERGE is broken into different consent classes. We can conduct analyses on hmb, hmb-gso-nic, | ||
+ | |||
+ | |||
+ | emerge.hmb <- read.table("/ | ||
+ | emerge.hmb.genos <- read.table("/ | ||
+ | |||
+ | |||
+ | emerge.hmb.gso.nic <- read.table("/ | ||
+ | emerge.hmb.gso.nic.genos <- read.table("/ | ||
+ | |||
+ | |||
+ | emerge.hmb.gso <- read.table("/ | ||
+ | emerge.hmb.gso.genos <- read.table("/ | ||
+ | |||
+ | |||
+ | ### Merge all files above according to SUBJID, which is used in the | ||
+ | ### genotype files. | ||
+ | |||
+ | emerge <- merge(emerge.hmb, | ||
+ | emerge <- merge(d, emerge.hmb.gso.nic, | ||
+ | |||
+ | ### SMOKING INITIATION | ||
+ | ### | ||
+ | ### The eMERGE variable name is SMOKING_STATUS | ||
+ | ### C65108 = never smoker | ||
+ | ### C67147 = current smoker | ||
+ | ### C67148 = past smoker | ||
+ | ### C67151 = Unknown if ever smoked | ||
+ | ### | ||
+ | ### Descriptives: | ||
+ | ### | ||
+ | ### table(emerge$SMOKING_STATUS) | ||
+ | ### | ||
+ | ### C65108 C67147 C67148 C67151 | ||
+ | ### | ||
+ | |||
+ | si <- emerge$SMOKING_STATUS | ||
+ | si[si == " | ||
+ | si[si == " | ||
+ | si[si != 1 & si != 2] <- NA | ||
+ | |||
+ | ### SMOKING Cessation | ||
+ | ### | ||
+ | ### Current == 2 & Former == 1 in GSCAN. This is already the case for these data. | ||
+ | |||
+ | sc <- emerge$SMOKING_STATUS | ||
+ | sc[sc == " | ||
+ | sc[sc == " | ||
+ | sc[sc != 1 & sc != 2] <- NA | ||
+ | |||
+ | |||
+ | ### eMERGE age variable is tricky because there is no obvious age at | ||
+ | ### assessment. We will use their " | ||
+ | ### approximation. | ||
+ | ### 1=1900-1919; | ||
+ | ### | ||
+ | ### Descriptives: | ||
+ | ### | ||
+ | ### table(emerge$DECADE_BIRTH) | ||
+ | ### | ||
+ | ### | ||
+ | ### | ||
+ | birthyear <- emerge$DECADE_BIRTH | ||
+ | birthyear[birthyear == " | ||
+ | birthyear[birthyear == " | ||
+ | |||
+ | |||
+ | ### SEX | ||
+ | sex <- emerge$SEX | ||
+ | sex[sex == " | ||
+ | sex[sex == " | ||
+ | |||
+ | |||
+ | ### Scott decided not to correct for additional case-control variables | ||
+ | ### given what appears to be a highly complex sample and uncertainty | ||
+ | ### about the best course of action to account for disease status in | ||
+ | ### conducting smoking analyses. | ||
+ | |||
+ | phenotypes <- data.frame(fid = emerge$SUBJID, | ||
+ | iid = emerge$SUBJID, | ||
+ | patid = " | ||
+ | matid = " | ||
+ | sex = sex, | ||
+ | si = si, | ||
+ | sc = sc) | ||
+ | |||
+ | phenotypes[is.na(phenotypes)] <- " | ||
+ | |||
+ | write.table(phenotypes, | ||
+ | "/ | ||
+ | row.names=F, | ||
+ | quote = F, | ||
+ | sep=" | ||
+ | |||
+ | |||
+ | covariates | ||
+ | iid = emerge$SUBJID, | ||
+ | patid = " | ||
+ | matid = " | ||
+ | sex = sex, | ||
+ | birthyear = birthyear) | ||
+ | |||
+ | covariates[is.na(covariates)] <- " | ||
+ | |||
+ | write.table(covariates, | ||
+ | "/ | ||
+ | row.names=F, | ||
+ | quote = F, | ||
+ | sep=" | ||
Line 1040: | Line 1149: | ||
- | beagess_data <- read.table("/ | + | options(stringsAsFactors=F) |
- | geno_data <- read.table("/ | + | |
- | geno_data <- geno_data[-c(6)] | + | |
- | geno_data$V3[geno_data$V3 ==" | + | |
- | geno_data$V4[geno_data$V4 ==" | + | |
- | ### IMPORTANT: All columns in original phenotype file are shifted one column to the the right, so the label of the column does not match the data in that column!!! | + | beagess |
- | si <- beagess_data$BMI | + | |
- | sc <- beagess_data$BMI | + | |
- | age <- beagess_data$sex | + | |
- | beagess_data | + | genos <- read.table("/ |
- | colnames(beagess_data)[13] <- " | + | |
- | ### SMOKING INITIATION | + | ### The file reads into R incorrectly because of a weird trailing tab |
+ | ### character in the data file, so use the below code to shift column | ||
+ | ### names to the correct column. | ||
+ | names(beagess)[1] <- " | ||
+ | beagess$dbGaP_Subject_ID <- row.names(beagess) | ||
+ | beagess$assocEABEvsCO <- NULL | ||
+ | names(beagess) <- c(names(beagess)[2: | ||
+ | |||
+ | |||
+ | ### SMOKING INITIATION | ||
### | ### | ||
### BEAGESS variable name is " | ### BEAGESS variable name is " | ||
Line 1062: | Line 1172: | ||
### -99 = not consented | ### -99 = not consented | ||
### -9 = Missing | ### -9 = Missing | ||
- | ### 0 = Never | + | ### 0 = Never |
### 1 = Former | ### 1 = Former | ||
### 2 = Current | ### 2 = Current | ||
- | ### | + | ### |
- | ### Descriptives: | + | ### Descriptives: |
### | ### | ||
### > table(beagess_data$BMI) | ### > table(beagess_data$BMI) | ||
- | ### | + | ### |
- | ### -99 | + | ### -99 |
- | ### 494 2051 1515 2288 575 | + | ### 494 2051 1515 2288 575 |
### | ### | ||
### > summary(beagess_data$BMI) | ### > summary(beagess_data$BMI) | ||
- | ### Min. 1st Qu. Median | + | ### Min. 1st Qu. Median |
- | ### -99.000 | + | ### -99.000 |
- | beagess_data$BMI[beagess_data$BMI | + | si <- beagess$cig_smk_status |
- | beagess_data$BMI[beagess_data$BMI | + | si[si == 1 | si == 2] <- 2 |
- | beagess_data$BMI[beagess_data$BMI == 1] = 2 | + | si[si == 0] <- 1 |
- | beagess_data$BMI[beagess_data$BMI == 0] = 1 | + | si[si != 1 & si != 2] <- NA |
- | ### SMOKING Cessation | + | ### SMOKING Cessation |
### | ### | ||
### BEAGESS variable name is " | ### BEAGESS variable name is " | ||
Line 1089: | Line 1199: | ||
### -99 = not consented | ### -99 = not consented | ||
### -9 = Missing | ### -9 = Missing | ||
- | ### 0 = Never | + | ### 0 = Never |
### 1 = Former | ### 1 = Former | ||
### 2 = Current | ### 2 = Current | ||
- | ### | ||
- | ### Descriptives: | ||
### | ### | ||
- | ### > table(beagess_data$BMI) | + | ### Descriptives: |
- | ### | + | ### table(beagess$cig_smk_status) |
- | ### -99 | + | ### -99 |
- | ### 494 2051 1515 2288 575 | + | ### 494 2051 1515 2288 575 |
### | ### | ||
- | ### > summary(beagess_data$BMI) | + | ### Current == 2 & Former == 1 in GSCAN. This is already the case for these data. |
- | ### Min. 1st Qu. Median | + | |
- | ### -99.000 | + | |
- | beagess_data$sc[beagess_data$sc == -99] = " | + | sc <- beagess$cig_smk_status |
- | beagess_data$sc[beagess_data$sc == -9] = " | + | sc[sc == 0 | sc == -9 | sc == -99] <- NA |
- | beagess_data$sc[beagess_data$sc == 0] = " | + | |
- | ### AGE | + | ### BEAGESS variable name is " |
- | ### | + | |
- | ### BEAGESS variable name is " | + | |
### " | ### " | ||
- | ### | + | ### |
### -9 = Missing | ### -9 = Missing | ||
### 1 = 15-29 years of age | ### 1 = 15-29 years of age | ||
Line 1129: | Line 1232: | ||
### 14 = 90-100 years of age | ### 14 = 90-100 years of age | ||
### | ### | ||
- | ### Descriptives: | + | ### Descriptives: |
### | ### | ||
- | ### > table(beagess_data$sex) | + | ### > table(beagess_data$agegpcat) |
- | ### | + | |
- | ### -9 1 2 3 4 5 6 7 8 9 | + | |
- | ### 27 | + | |
### | ### | ||
- | ### > summary(beagess_data$BMI) | + | ### -9 1 2 3 4 5 6 7 8 9 |
- | ### Min. 1st Qu. Median | + | ### 27 |
- | ### -9.000 | + | ### |
+ | ### This is fine, but change the -9 to NA | ||
+ | beagess$agegpcat[beagess$agegpcat == -9] <- NA | ||
- | beagess_data$sex[beagess_data$sex== 1] = 22 | ||
- | beagess_data$sex[beagess_data$sex== 2] = 32 | ||
- | beagess_data$sex[beagess_data$sex== 3] = 37 | ||
- | beagess_data$sex[beagess_data$sex== 4] = 42 | ||
- | beagess_data$sex[beagess_data$sex== 5] = 47 | ||
- | beagess_data$sex[beagess_data$sex== 6] = 52 | ||
- | beagess_data$sex[beagess_data$sex== 7] = 57 | ||
- | beagess_data$sex[beagess_data$sex== 8] = 62 | ||
- | beagess_data$sex[beagess_data$sex== 9] = 67 | ||
- | beagess_data$sex[beagess_data$sex== 10] = 72 | ||
- | beagess_data$sex[beagess_data$sex== 11] = 77 | ||
- | beagess_data$sex[beagess_data$sex== 12] = 82 | ||
- | beagess_data$sex[beagess_data$sex== 13] = 87 | ||
- | beagess_data$sex[beagess_data$sex== 14] = 95 | ||
- | beagess_data$sex[beagess_data$sex== -9] = " | ||
- | beagess_data$sex <- as.numeric(beagess_data$sex) | ||
- | ### rename SUBJID column to dbGaP_Subject_ID in genotype file so it matches phenotype file where SUBJID is misslabelled as dbGaP_Subject_ID | ||
- | colnames(geno_data) <- c(" | ||
- | ### merge geno and pheno files | + | ### looks like BEAGESS uses their own internal " |
- | phen <- merge(geno_data,beagess_data, by="dbGaP_Subject_ID", | + | ### genotype files, so we'll use that in our pedigree |
- | phen <- phen[,c(1:5,7,10,17)] | + | phenotypes |
+ | iid = beagess$SUBJECT_ID, | ||
+ | | ||
+ | | ||
+ | sex = beagess$sex, | ||
+ | si = si, | ||
+ | sc = sc) | ||
- | phen <- phen[,c(2, | + | phenotypes[is.na(phenotypes)] <- "x" |
- | colnames(phen) | + | |
+ | write.table(phenotypes, | ||
+ | "/ | ||
+ | row.names=F, | ||
+ | quote = F, | ||
+ | sep=" | ||
- | phenotypes <- data.frame(famid=phen$famid, | ||
- | colnames(phenotypes) <- c(" | ||
- | pcs <- read.table("/ | + | covariates |
- | colnames(pcs) [1] <- "dbGaP_Subject_ID" | + | iid = beagess$SUBJECT_ID, |
+ | patid = "x", | ||
+ | matid = "x", | ||
+ | sex = beagess$sex, | ||
+ | age = beagess$age, | ||
+ | age2 = beagess$age^2, | ||
+ | Barretts.esophagus.case_v_control = beagess$assocBEvsCO, | ||
+ | Esophageal.adenocarcinoma.case_v_control = beagess$assocEAvsCO) | ||
- | gscan.phenotypes <- merge(phenotypes, | + | covariates[is.na(covariates)] <- " |
- | gscan.phenotypes[is.na(gscan.phenotypes)] <- " | + | |
- | ### EUROPEANS - entire sample is EUR | + | write.table(covariates, |
- | phenotypes.EUR.ped <- subset(gscan.phenotypes, | + | |
- | write.table(phenotypes.EUR.ped, file=" | + | |
- | + | quote = F, | |
- | covariates.EUR.ped <- subset(gscan.phenotypes, | + | |
- | write.table(covariates.EUR.ped, | + | |