Editing GSCAN dbGaP (section)

===Phenotypes===

<syntaxhighlight lang="rsplus">
options(stringsAsFactors=F)

beagess <- read.table("/work/KellerLab/GSCAN/dbGaP/Barretts/PhenoGenotypeFiles/RootStudyConsentSet_phs000869.BEAGESS.v1.p1.c1.GRU-MDS/PhenotypeFiles/phs000869.v1.pht004610.v1.p1.c1.BEAGESS_Subject_Phenotypes.GRU-MDS.txt.gz", header=TRUE, sep="\t", stringsAsFactors=F)

genos <- read.table("/work/KellerLab/GSCAN/dbGaP/Barretts/PhenoGenotypeFiles/RootStudyConsentSet_phs000869.BEAGESS.v1.p1.c1.GRU-MDS/GenotypeFiles/phg000580.v1.NCI_BEAGESS.genotype-calls-matrixfmt.c1.GRU-MDS.update/BEAGESS_dbGaP_29Jan2015.fam", header=FALSE, stringsAsFactors=F)


### 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] <- "XXX"
beagess$dbGaP_Subject_ID <- row.names(beagess)
beagess$assocEABEvsCO <- NULL
names(beagess) <- c(names(beagess)[2:11], "assocEABEvsCO", "dbGaP_Subject_ID")


### SMOKING INITIATION
###
### BEAGESS variable name is "cig_smk_status" and is listed under the column "BMI".
###    "Cigarette smoking status."
###    Response options are "-99", "-9", "0", "1" and "2".
###       -99 = not consented
###       -9 = Missing
###        0 = Never
###        1 = Former
###        2 = Current
###
### Descriptives:
###
### > table(beagess_data$BMI)
###
### -99   -9    0    1    2
### 494 2051 1515 2288  575
###
### > summary(beagess_data$BMI)
### Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
### -99.000  -9.000   0.000  -9.234   1.000   2.000

si <- beagess$cig_smk_status
si[si == 1 | si == 2] <- 2
si[si == 0] <- 1
si[si != 1 & si != 2] <- NA

### SMOKING Cessation
###
### BEAGESS variable name is "cig_smk_status" and is listed under the column "BMI".
###    "Cigarette smoking status."
###    Response options are "-99", "-9", "0", "1" and "2".
###       -99 = not consented
###       -9 = Missing
###        0 = Never
###        1 = Former
###        2 = Current
###
### Descriptives:
### table(beagess$cig_smk_status)
###  -99   -9    0    1    2 
###  494 2051 1515 2288  575 
###
### Current == 2 & Former == 1 in GSCAN. This is already the case for these data.

sc <- beagess$cig_smk_status
sc[sc == 0 | sc == -9 | sc == -99] <- NA

### BEAGESS variable name is "agegpcat"
###   "Age in 5 year categories"
###   Response options are integers 1 through 14 or -9.
###       -9 = Missing
###       1 = 15-29 years of age
###       2 = 30-34 years of age
###       3 = 35-39 years of age
###       4 = 40-44 years of age
###       5 = 45-49 years of age
###       6 = 50-54 years of age
###       7 = 55-59 years of age
###       8 = 60-64 years of age
###       9 = 65-69 years of age
###       10 = 70-74 years of age
###       11 = 75-79 years of age
###       12 = 80-84 years of age
###       13 = 85-89 years of age
###       14 = 90-100 years of age
###
### Descriptives:
###
### > table(beagess_data$agegpcat)
###
### -9    1    2    3    4    5    6    7    8    9   10   11   12   13   14
### 27   19   48  112  190  378  657  958 1134 1238 1018  794  246   87   17
###
### This is fine, but change the -9 to NA
beagess$agegpcat[beagess$agegpcat == -9] <- NA



### looks like BEAGESS uses their own internal "SUBJECT_ID" in the
### genotype files, so we'll use that in our pedigree files
phenotypes <- data.frame(fid = beagess$SUBJECT_ID,
                         iid = beagess$SUBJECT_ID,
                         patid = "x",
                         matid = "x",
                         sex = beagess$sex,
                         si = si,
                         sc = sc)

phenotypes[is.na(phenotypes)] <- "x"

write.table(phenotypes,
            "/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/BEAGESS/GSCAN_BEAGESS_phenotypes.ped",
            row.names=F,
            quote = F,
            sep="\t")


covariates  <- data.frame(fid = beagess$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)

covariates[is.na(covariates)] <- "x"

write.table(covariates,
                        "/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/BEAGESS/GSCAN_BEAGESS_covariates.ped",
                        row.names=F,
                        quote = F,
                        sep="\t") ### save table                 

</syntaxhighlight>