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/01/05 18:30] hannah_young /* Phenotypes */ |
keller_and_evans_lab:gscan_db_ga_p [2017/02/14 21:14] scott /* Phenotypes */ |
||
---|---|---|---|
Line 899: | Line 899: | ||
+ | ### Date: Feb 13 2017 | ||
+ | ### Author: Scott Vrieze | ||
- | options(stringsAs Factors=F) | + | options(stringsAsFactors=F) |
- | ### Load in dataset ### | + | ### Load in dataset ### |
- | ninds <- read.table(gzfile("/ | + | ninds <- read.table("/ |
+ | | ||
+ | ### 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(ninds)[1] <- " | ||
+ | ninds$dbGaP_Subject_ID <- row.names(ninds) | ||
+ | ninds$smokingStatus <- NULL | ||
+ | names(ninds) <- c(names(ninds)[2: | ||
+ | | ||
- | ### The file reads into R incorrectly, | + | ### subset |
- | tmp <- colnames(ninds) | + | pheno <- subset(ninds, select=c(" |
- | tmp2 <- tmp[c(2:19)] | + | |
- | tmp2 <- c(tmp2, "tmp") | + | |
- | colnames(ninds) <- tmp2 | + | |
- | ### subset the only variables needed | ||
- | phenotypes <- c(" | ||
- | pheno <- ninds2[phenotypes] | ||
### | ### | ||
Line 922: | Line 927: | ||
### NINDS variable is “smokingStatus” | ### NINDS variable is “smokingStatus” | ||
- | ### Variables are “NEVER”, | + | ### Variables are “NEVER”, |
- | ### | + | ### |
- | ### ’CURRENT' | + | ### ’CURRENT' |
+ | ### ' | ||
+ | ### but no smoking within the last 30 ### days; ’NEVER' | ||
+ | ### less than 100 cigarettes smoked in one's lifetime. | ||
### table(pheno$smokingStatus) | ### table(pheno$smokingStatus) | ||
### | ### | ||
- | ### CURRENT | + | ### CURRENT |
- | ### | + | ### |
- | pheno[smokingStatus == " | + | si <- pheno$smokingStatus |
- | pheno[pheno$smokingStatus | + | si[si == " |
+ | si[si == " | ||
+ | si[si != " | ||
+ | si <- as.numeric(si) | ||
- | ### table(pheno$V4) | + | ### table(si) |
### | ### | ||
- | ### 1 2 | + | ### 1 2 |
- | ### 2401 1864 | + | ### 2401 1864 |
- | + | ||
- | colnames(pheno)[4] <- " | + | |
Line 946: | Line 955: | ||
### NINDS variable is “smokingStatus” | ### NINDS variable is “smokingStatus” | ||
- | ### Variables are “NEVER”, | + | ### Variables are “NEVER”, |
- | ### | + | ### |
- | ### ’CURRENT' | + | ### ’CURRENT' |
+ | ### ' | ||
+ | ### no smoking within the last 30 days; ’NEVER' | ||
+ | ### 100 cigarettes smoked in one's lifetime. | ||
### table(pheno$smokingStatus) | ### table(pheno$smokingStatus) | ||
### | ### | ||
- | ### CURRENT | + | ### CURRENT |
- | ### | + | ### |
- | pheno[pheno$smokingStatus == " | + | sc <- pheno$smokingStatus |
- | pheno[pheno$smokingStatus | + | sc[sc == " |
+ | sc[sc == " | ||
+ | sc[sc != " | ||
+ | sc <- as.numeric(sc) | ||
### table(pheno$V5) | ### table(pheno$V5) | ||
### | ### | ||
- | ### 1 2 | + | ### 1 2 |
- | ### 1137 727 | + | ### 1137 727 |
- | colnames(pheno)[5] <- " | ||
### | ### | ||
Line 972: | Line 986: | ||
### table(pheno$gender) | ### table(pheno$gender) | ||
### | ### | ||
- | ### F M | + | ### F M |
- | ### 2627 1952 | + | ### 2627 1952 |
+ | sex <- pheno$gender | ||
+ | sex[sex == " | ||
+ | sex[sex == " | ||
+ | sex <- as.numeric(sex) | ||
- | pheno[gender==" | + | ###----------------### |
- | pheno[gender==" | + | ### Write to files ### |
+ | ### | ||
- | colnames(pheno)[6] | + | ### This study uses the " |
+ | ### so use that here, instead of the dbGaP_Subject_ID | ||
+ | phenotypes | ||
+ | iid = pheno$subject_id, | ||
+ | patid = "x", | ||
+ | matid = " | ||
+ | sex = sex, | ||
+ | si = si, | ||
+ | sc = sc) | ||
- | ######################### | + | write.table(phenotypes, |
- | ######################### | + | "/ |
- | ###### FORM TABLES ###### | + | |
- | ######################### | + | quote = F, |
- | ######################### | + | |
- | cbind(pheno, | ||
- | unknown <- NA | ||
- | ### PHENOTYPE TABLE ### | + | covariates |
+ | iid = pheno$subject_id, | ||
+ | patid = " | ||
+ | matid = " | ||
+ | sex = sex, | ||
+ | age = pheno$age, | ||
+ | age2 = pheno$age^2, | ||
+ | AFFECTION_STATUS = pheno$AFFECTION_STATUS) | ||
- | phenotype <- new[,c(7,8,6,4,5)] | + | write.table(covariates, |
- | colnames(phenotype)[1] <- "iid" | + | "/ |
+ | row.names=F, | ||
+ | quote = F, | ||
+ | sep="\t") ### save table | ||
- | ### generate fid, matid, patid columns | ||
- | phenotype$fid <- unknown | ||
- | phenotype$matid <- unknown | ||
- | phenotype$patid <- unkown | ||
- | phenotype <- phenotype[, | ||
- | phenotype <- replace(phenotype, | ||
- | write.table(phenotype, | ||
+ | ====== BEAGESS ====== | ||
- | ### COVARIATE TABLE ### | ||
- | covariate | ||
- | covariate$age2 <- covariate$age^2 | ||
- | colnames(covariate)[1] <- " | ||
- | ### generate fid, matid, patid columns | + | ===== Phenotypes ===== |
- | covariate$fid <- unknown | + | |
- | covariate$matid <- unknown | + | |
- | covariate$patid <- unknown | + | |
- | covariate <- covariate[, | + | |
- | + | ||
- | covariate <- replace(covariate, | + | |
- | write.table(covariate, | + | |
- | == BEAGESS == | ||
+ | options(stringsAsFactors=F) | ||
- | ===Phenotypes=== | + | beagess <- read.table("/ |
- | beagess_data <- read.table("/ | + | genos <- read.table("/ |
- | geno_data | + | |
- | 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!!! | ||
- | si <- beagess_data$BMI | ||
- | sc <- beagess_data$BMI | ||
- | age <- beagess_data$sex | ||
- | beagess_data <- cbind(beagess_data, beagess_data[,6]) | + | ### The file reads into R incorrectly because of a weird trailing tab |
- | colnames(beagess_data)[13] <- "sc" | + | ### 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 | + | ### SMOKING INITIATION |
### | ### | ||
### BEAGESS variable name is " | ### BEAGESS variable name is " | ||
Line 1048: | Line 1063: | ||
### -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 1075: | Line 1090: | ||
### -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 1115: | Line 1123: | ||
### 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 | + | ### -9 1 2 3 4 5 6 7 8 9 |
- | ### 27 | + | ### 27 |
### | ### | ||
- | ### > summary(beagess_data$BMI) | + | ### This is fine, but change the -9 to NA |
- | ### Min. 1st Qu. Median | + | beagess$agegpcat[beagess$agegpcat == -9] <- NA |
- | ### -9.000 | + | |
- | 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, | + | |
Line 1175: | Line 1177: | ||
====== Jackson Heart Study ====== | ====== Jackson Heart Study ====== | ||
- | |||
- | |||
- | ===== File Paths ===== | ||
- | |||
- | |||
- | **HMB consent phenotypes** | ||
- | |||
- | / | ||
- | HS_CARe_Subject_Phenotypes.HMB-IRB.txt.gz | ||
- | |||
- | **NPU consent phenotypes** | ||
- | |||
- | / | ||
- | c1.JHS_CARe_Subject_Phenotypes.HMB-IRB-NPU.txt.gz | ||
- | |||
- | **HMB consent genotypes** | ||
- | |||
- | / | ||
- | ls-matrixfmt.c3.HMB-IRB/ | ||
- | |||
- | **NPU consent genotypes** | ||
- | |||
- | / | ||
- | -calls-matrixfmt.c1.HMB-IRB-NPU/ | ||
- | |||
- | **Subject Sample Mapping File** | ||
- | |||
- | / | ||
- | ULTI/ | ||
Line 1233: | Line 1206: | ||
names(covariates) [2] <- " | names(covariates) [2] <- " | ||
names(covariates) [6] <- " | names(covariates) [6] <- " | ||
+ | |||
+ | ### | ||
+ | ### Smoking initiation ### | ||
+ | ### | ||
+ | ### | ||
+ | ### JHSCare variables are " | ||
+ | ### | ||
+ | ### | ||
+ | ### | ||
+ | ### Response options are | ||
+ | ### 0 - No | ||
+ | ### 1 - Yes | ||
+ | ### | ||
+ | ### table(si) | ||
+ | ### | ||
+ | ### 0 1 x | ||
+ | ### 1206 537 9 | ||
current.smoker <- subset(mapped_geno_pcs, | current.smoker <- subset(mapped_geno_pcs, | ||
Line 1253: | Line 1243: | ||
} | } | ||
mapped_geno_pcs_phen <- cbind(si, mapped_geno_pcs) | mapped_geno_pcs_phen <- cbind(si, mapped_geno_pcs) | ||
+ | |||
+ | |||
+ | ### | ||
+ | ### Smoking Cessation ### | ||
+ | ### | ||
+ | ### | ||
+ | ### JHSCare variables are " | ||
+ | ### | ||
+ | ### | ||
+ | ### | ||
+ | ### Response options are | ||
+ | ### 0 - No | ||
+ | ### 1 - Yes | ||
+ | ### | ||
+ | ### table(sc) | ||
+ | ### | ||
+ | ### 0 1 x | ||
+ | ### 1206 537 9 | ||
+ | |||
current.smoker <- subset(mapped_geno_pcs, | current.smoker <- subset(mapped_geno_pcs, | ||
Line 1277: | Line 1286: | ||
names(phenotypes)[1] <- " | names(phenotypes)[1] <- " | ||
- | write.table(covariates, | + | write.table(covariates, |
- | write.table(phenotypes, | + | write.table(phenotypes, |
- | =Genotype Processing= | ||
- | ===Pre-Phasing QC=== | + | |
+ | ======= Genotype Processing ======= | ||
+ | |||
+ | |||
+ | |||
+ | ===== Pre-Phasing QC ===== | ||
QC parameters that we chose: | QC parameters that we chose: | ||
MAF > 0.01 | MAF > 0.01 | ||
Line 1292: | Line 1306: | ||
HWE = 0.05 / number of markers but greater than 5e-8 | HWE = 0.05 / number of markers but greater than 5e-8 | ||
- | < | + | |
## To update the strand builds: http:// | ## To update the strand builds: http:// | ||
## Check strands against latest 1000G: http:// | ## Check strands against latest 1000G: http:// | ||
Line 1321: | Line 1335: | ||
| | ||
| | ||
+ | |||