More information about the files used for UKBiobank are here. In brief, we used the UK10K + 1kgp3 imputed vcfs provided by UKBionank and added in dosages w/ this python script:
import gzip, argparse, re, os, datetime from subprocess import Popen, PIPE
def add_dosage(pair):
a, b = pair
probs = b.split(b',')
dose = float(probs[1]) + (float(probs[2]) * 2)
return a + b':' + str(dose).encode('ascii') + b':' + b
def gziplines(fname):
f = Popen(['zcat', fname], stdout=PIPE)
for line in f.stdout:
yield line
parser = argparse.ArgumentParser() parser.add_argument('inputVCF', help = 'The path to the VCF') args = parser.parse_args()
flag = False
for line in gziplines(args.inputVCF):
if line.startswith(b'#'):
os.write(1, line.rstrip() + b'\n')
if not flag:
os.write(1, b'##FORMAT=<ID=DS,Number=1,Type=Float,Description="Genotype Dosages">\n')
os.write(1, b'##Dosages added using the script add.dosages.subprocess.py at ' +
str(datetime.datetime.now()).encode('ascii') + b'\n')
flag = True
else:
elements = re.split(b'\t|:', line.rstrip())
first8 = elements[:8]
genotypes = elements[10:]
form = b'GT:DS:GP'
genotypes_split = zip(genotypes[::2], genotypes[1::2])
try:
dose_genos = [add_dosage(pair) for pair in genotypes_split]
except (ValueError, IndexError) as e:
os.write(2, "\n" + line)
os.write(2, line + "\n" + args.inputVCF + "\n\n")
raise e
os.write(1, b'\t'.join(first8) + b'\t' + form + b'\t' + b'\t'.join(dose_genos) + b'\n')
options(stringsAsFactors=F)
### Offspring cohort off ← read.table(gzfile(“/work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/RootStudyConsentSet_phs000007.Framingham.v28.p10.c1.HMB-IRB-MDS/PhenotypeFiles/phs000007.v28.pht000030.v7.p10.c1.ex1_1s.HMB-IRB-MDS.txt.gz”), header=T, sep=“\t”)
### Generation 3 cohort g3 ← read.table(gzfile(“/work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/RootStudyConsentSet_phs000007.Framingham.v28.p10.c1.HMB-IRB-MDS/PhenotypeFiles/phs000007.v28.pht000074.v9.p10.c1.ex3_1s.HMB-IRB-MDS.txt.gz”), header=T, skip=10, sep=“\t”, fill=T)
###———————### ### Cigarettes per Day ### ###———————###
### Offspring Cohort Exam 1 ### Framingham variable name is “A102”. ### “Usual number of cigarettes smoked (now or formerly)” ### Response option is an integer value ### ### table(off$A102) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 1208 50 29 44 26 47 31 16 30 3 175 1 20 1 9 94 ### 16 17 18 19 20 22 24 25 26 27 28 30 34 35 40 42 ### 5 6 50 1 556 10 3 39 1 1 2 214 1 15 147 2 ### 45 50 60 80 88 90 ### 4 21 19 1 14 2 ### ### summary(off$A102) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.00 0.00 7.00 12.11 20.00 90.00 8
cpd.off ← off$A102 cpd.off[cpd.off==0] ← NA ### Set those who smoke “0” cpd to NA cpd.off[cpd.off > 60] ← 60 ### Set those who report smoking >40 cpd to 40
cpd.off[cpd.off ⇐ 5 & cpd.off >= 1] ← 1 cpd.off[cpd.off ⇐ 15 & cpd.off >= 6] ← 2 cpd.off[cpd.off ⇐ 25 & cpd.off >= 16] ← 3 cpd.off[cpd.off ⇐ 35 & cpd.off >= 26] ← 4 cpd.off[cpd.off >= 36 & cpd.off ⇐ 60] ← 5 cpd.off[cpd.off > 60 | is.na(cpd.off)] ← NA
### G3 cohort Exam 1 ### Framingham variable name is “G3A074 ### “IF EVER SMOKED CIGS REGULARLY: ON THE AVERAGE ### OF THE ENTIRE TIME YOU SMOKED, HOW MANY CIGARETTES ### DID YOU SMOKE PER DAY?” ### Response option is an integer value ### ### table(g3$G3A074) ### ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 3180 57 14 15 20 30 8 13 7 2 122 2 16 3 2 60 ### 16 17 18 19 20 21 22 23 24 25 28 30 35 40 45 50 ### 3 3 10 1 287 2 1 1 1 12 1 53 1 47 1 12 ### 55 60 80 ### 1 6 1 ### ### summary(g3$G3A074) ### Min. 1st Qu. Median Mean 3rd Qu. Max. ### 0.000 0.000 0.000 3.528 0.000 80.000
cpd.g3 ← g3$G3A074 cpd.g3[cpd.g3==0] ← NA ### Set those who smoke “0” cpd to NA cpd.g3[cpd.g3 > 60] ← 60 ### Set those who report smoking >40 cpd to 40
cpd.g3[cpd.g3 ⇐ 5 & cpd.g3 >= 1] ← 1 cpd.g3[cpd.g3 ⇐ 15 & cpd.g3 >= 6] ← 2 cpd.g3[cpd.g3 ⇐ 25 & cpd.g3 >= 16] ← 3 cpd.g3[cpd.g3 ⇐ 35 & cpd.g3 >= 26] ← 4 cpd.g3[cpd.g3 >= 36 & cpd.g3 ⇐ 60] ← 5 cpd.g3[cpd.g3 > 60 | is.na(cpd.g3)] ← NA
###——————–### ### Smoking initiation ### ###——————–### ### Offspring Cohort Exam 1 ### Framingham variable is “A99” ### “Smokes cigarettes: Yes(now)/No/Former” ### Response option is fill-in ### 0 is no, 1 is yes (now), 2 is former ### ### table(off$A99) ### ### 0 1 2 ### 1203 1126 566
si.off ← off$A99 si.off[si.off==1] ← 2 ### Set current smoker “1” to smoker “2” si.off[si.off==0] ← 1 ### Set never smoker “0” to nonsmoker “1”
### Former smoker "2" remains as a "2" as smoker
### G3 Cohort Exam 1 ### Framingham variables is “G3A070” ### G3A070 - HAVE YOU EVER SMOKED CIGARETTES REGULARLY? (NO MEANS ### LESS THAN 20 PACKS OF CIGARETTES OR 12 OZ OF TOBACCO IN A ### LIFETIME OR LESS THAN 1 CIGARETTE A DAY FOR A YEAR.)” ### Response options: ### “0” is no, “1” is yes, “0” is unknown ### ### table(g3$G3A070) ### 0 1 20 ### 2137 1629 1 si.g3 ← g3$G3A070 si.g3[si.g3==20] ← NA si.g3[si.g3==1] ← 2 si.g3[si.g3==0] ← 1
###———————### ### Smoking cessation ### ###———————###
### Offspring Cohort Exam 1 ### Framingham variable is “A99” ### Smokes cigarettes: Yes(now)/No/Former“ ### Response option is fill-in ### 0 is no, 1 is yes (now), 2 is former ### ### table(off$A99) ### ### 0 1 2 ### 1203 1126 566
sc.off ← off$A99 sc.off[sc.off==2] ← 3 ### temporarily set former smoker “2” as “3” sc.off[sc.off==1] ← 2 ### set current smoker “1” as current smoker “2” sc.off[sc.off==3] ← 1 ### set former smoker “3” as former smoker “1” sc.off[sc.off==0] ← NA ### code all non smokers “0” as NA
### G3 Cohort Exam 1 ### Framingham variables are “G3A070” and “G3A072” ### Response option is fill-in ### G3A070 - HAVE YOU EVER SMOKED CIGARETTES REGULARLY? (NO MEANS LESS THAN 20 PACKS OF CIGARETTES OR 12 OZ OF ###TOBACCO IN A LIFETIME OR LESS THAN 1 CIGARETTE A DAY FOR A YEAR.) ###“0” is no, “1” is yes, “20” is unknown ### G3A072 - IF EVER SMOKED CIGARETTES REGULARLY: DO YOU NOW SMOKE CIGARETTES (AS OF 1 MONTH AGO)? ###“0” is no or never smoked, “1” is yes ### ### table(g3$G3A070) ### 0 1 20 ### 2137 1629 1 ### ### table(g3$G3A072) ### 0 1 ### 3179 585 sc.g3 ← g3$G3A072 sc.g3[sc.g3==1] ← 2 ## These are our current smokers sc.g3[sc.g3==0] ← 1 sc.g3 ← ifelse(sc.g3==1 & (g3$G3A070 == 0 | is.na(g3$G3A070)), NA, sc.g3)
###———————### ### Age of initiation ### ###———————###
### Offspring Cohort Exam 1 ### Framingham variable “A100” ### Age started smoking cigarettes regularly\ ### Response option is an integer value, 88 is doesn't smoke regularly ### table(off$A100) ### ### 0 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 ### 6 2 2 2 8 4 24 37 61 147 266 233 330 166 164 85 42 26 13 25 ### 26 27 28 29 30 31 33 34 35 37 38 39 40 44 45 46 51 88 ### 2 8 4 3 7 1 1 1 3 1 2 1 3 1 1 2 1 105 ###
ai.off ← off$A100 ai.off[ai.off < 10] ← NA ### Set AI less than 10 as missing\ ai.off[ai.off > 35] ← NA ### Set AI greater than 35 as missing \
### G3 Cohort Exam 1 ### Framingham variable “G3A075” ### IF EVER SMOKED CIGARETTES REGULARLY: HOW OLD WERE YOU WHEN YOU FIRST STARTED REGULAR CIGARETTE SMOKING? ### Response option is an integer value, 0 is never smoked ### table(g3$G3A075) ### ### 0 5 7 8 9 10 11 12 13 14 15 16 17 18 19 20 ### 2138 1 2 2 3 5 14 52 108 159 152 277 193 256 98 107 ### 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 ### 51 31 19 13 31 6 4 6 2 13 2 1 1 4 1 4 ### 37 38 40 42 46 50 ### 2 2 1 1 1 2 ###
ai.g3 ← g3$G3A075 ai.g3[ai.g3 < 10] ← NA ### Set AI less than 10 as missing ai.g3[ai.g3 > 35] ← NA ### Set AI greater than 35 as missing
###——————-### ### Drinks per week ### ###——————-###
### Offspring Cohort Exam 1 ### Framingham variables “A111”, “A112”, “A113” ### A111 - Beer (bottles, cans or glasses per week) ### A112 - Wine (glasses per week) ### A113 - Cocktails, highballs, straight drinks (# per week) ### Response options are integer values. ### table(off$A111) ### ### 0 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 ### 1452 524 171 131 78 49 125 34 33 4 46 68 21 21 5 13 ### 20 21 24 25 26 28 30 35 36 40 42 48 49 50 70 90 ### 18 7 48 5 1 6 6 2 2 4 4 6 1 2 1 1 ### ### table(off$A112) ### ### 0 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 ### 1329 919 218 105 69 41 22 64 25 3 25 11 30 2 2 2 ### 20 21 24 28 30 32 50 ### 7 4 2 3 1 2 1 ### ### table(off$A113) ### ### 0 1 2 3 4 5 6 7 8 9 10 12 14 15 16 18 ### 861 1094 267 157 96 64 59 65 31 9 49 16 57 13 2 1 ### 19 20 21 24 25 28 30 35 38 42 50 55 60 61 ### 1 12 10 4 2 2 6 2 1 1 1 1 1 1 ### dpw.off ← rowSums(subset(off, select = c(“A111”, “A112”, “A113”)), na.rm=T) ### The above assumes all drinks are of equal ETOH content dpw.off[dpw.off == 0] ← NA dpw.off ← log(dpw.off)
### G3 Cohort Exam 1\ ### Framingham variables are “G3A115”, “G3A116”, “G3A119”, “G3A120”, ### “G3A123”, “G3A124”, “G3A127”, “G3A128”, ### “G3A131”, “G3A132” ### G3A115 - BEER: NUMBER OF BEER (12 OZ. BOTTLE, GLASS, CAN) YOU DRINK PER WEEK OVER THE PAST YEAR ### G3A116 - BEER: NUMBER OF BEER (12 OZ. BOTTLE, GLASS, CAN) YOU DRINK PER MONTH OVER THE PAST YEAR ### G3A119 - WHITE WINE: NUMBER OF WHITE WINE (4 OZ GLASS) YOU DRINK PER WEEK OVER THE PAST YEAR ### G3A120 - WHITE WINE: NUMBER OF WHITE WINE (4 OZ GLASS) YOU DRINK PER MONTH OVER THE PAST YEAR ### G3A123 - RED WINE: NUMBER OF RED WINE (4 OZ GLASS) YOU DRINK PER WEEK OVER THE PAST YEAR ### G3A124 - RED WINE: NUMBER OF RED WINE (4 OZ GLASS) YOU DRINK PER MONTH OVER THE PAST YEAR ### G3A127 - LIQUOR/SPIRITS: AVERAGE NUMBER OF LIQUOR/SPIRITS (1 1/4 OZ JIGGER) YOU DRINK PER WEEK OVER THE PAST YEAR ### G3A128 - LIQUOR/SPIRITS: AVERAGE NUMBER OF LIQUOR/SPIRITS (1 1/4 OZ JIGGER) YOU DRINK PER MONTH OVERTHE PAST YEAR ### G3A131 - OTHER BEVERAGE: AVERAGE NUMBER OF OTHER BEVERAGE YOU DRINK PER WEEK OVER THE PAST YEAR ### G3A132 - OTHER BEVERAGE: AVERAGE NUMBER OF OTHER BEVERAGE YOU DRINK PER MONTH OVER THE PAST YEAR ### ### *Participant was allowed to report alcohol consumption in either drinks per week or drinks per ### month. Therefore, to calculate total alcohol consumption, you must use both number of drinks per ### week and number of drinks per month (E.G. DRINKS PER MONTH = SUM(OF (DRINKS PER WEEK)(DRINKS PER ### MONTH/4)). ### ### Response options are integer values. ### all ← subset(g3, select = c(“G3A115”, “G3A116”, “G3A119”, “G3A120”, “G3A123”, “G3A124”, “G3A127”, “G3A128”,
"G3A131", "G3A132"))
names(all) ← c(“beerpw”, “beerpm”, “wwinepw”, “wwinepm”, “rwinepw”, “rwinepm”, “liqpw”, “liqpm”, “othpw”, “othpm”) w ← all[,c(1,3,5,7,9)] m ← all[,c(2,4,6,8,10)]/4
dpw.tmp ← data.frame(beer = rowSums(cbind(w$beerpw, m$beerpm), na.rm=T),
wwine = rowSums(cbind(w$wwinepw, m$wwinepm), na.rm=T),
rwine = rowSums(cbind(w$rwinepw, m$rwinepm), na.rm=T),
liq = rowSums(cbind(w$liqpw, m$liqpm), na.rm=T),
oth = rowSums(cbind(w$othpw, m$othpm), na.rm=T))
### WHITE WINE - convert from 4oz to 5oz drink dpw.tmp$wwine ← dpw.tmp$wwine*5/4 ### RED WINE - convert from 4oz to 5oz drink dpw.tmp$rwine ← dpw.tmp$rwine*5/4 ### LIQUOR - convert from 1.25 to 1.50 drink dpw.tmp$liq ← dpw.tmp$liq*1.5/1.25 ### OTHER DRINK is an unknown. leave as-is. dpw.g3 ← rowSums(dpw.tmp, na.rm=T)
### We have a few outliers, draw cuttoff at 70 drinks / week dpw.g3[dpw.g3 > 70] ← NA dpw.g3[dpw.g3 == 0] ← NA dpw.g3 ← log(dpw.g3 + .75)
###—————————### ### Drinker vs. non-drinker ### ###—————————### ### ### Use the variables already computed for dpw above. Assume ### anyone who is drinking less than 1 drink / month is a non-drinker x ← subset(off, select = c(“A111”, “A112”, “A113”)) index ← rep(NA, nrow(x)) for(i in 1:nrow(x)) {
if(is.na(x[i,1]) & is.na(x[i,2]) & is.na(x[i,3])) {
index[i] <- 1
} else {
index[i] <- 0
}
}
dnd.off ← ifelse(rowSums(subset(off, select = c(“A111”, “A112”, “A113”)), na.rm=T) == 0, 1, 2) dnd.off[index==1] ← NA
x ← all index ← rep(NA, nrow(x)) for(i in 1:nrow(x)) {
if(is.na(x[i,1]) & is.na(x[i,2]) & is.na(x[i,3])) {
index[i] <- 1
} else {
index[i] <- 0
}
}
dnd.g3 ← ifelse(rowSums(all, na.rm=T) == 0, 1, 2) dnd.g3[index==1] ← NA
###—————————-### ### Binge drinking in everyone ### ### (Only available for G3) ### ###—————————-### ### G3 Cohort Exam 1 ### Framingham variable is “G3A137” ### IF EVER CONSUMED ALCOHOL: WHAT WAS THE MAXIMUM NUMBER OF ### DRINKS YOU HAD IN A 24 HOUR PERIOD DURING THE PAST MONTH? ### Response option is an integer value ### ### summary(g3$G3A137) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 1.000 3.000 3.566 5.000 50.000 3
bde ← g3$G3A137 sex ← g3$G3A440 ### bde is defined differently for males & females ### sex == 1 is male, sex == 2 is female
tmp ← rep(NA, length(bde)) ### temporary vector to hold results and avoid conflicts
tmp[sex == 1 & bde < 5] ← 1 tmp[sex == 1 & bde >=5] ← 2 tmp[sex == 2 & bde < 4] ← 1 tmp[sex == 2 & bde >=4] ← 2
bde.g3 ← tmp
###—————————————–### ### Binge-drinking (in lifetime drinkers) ### ### (Only available for G3) ### ###—————————————–###
### G3 Cohort Exam 1 ### Framingham variable is “G3A112” ### HAVE YOU EVER CONSUMED ALCOHOLIC BEVERAGES (BEER, WINE, LIQUOR/SPIRITS)? \ ### 0 is no, 1 is yes
lifednd ← g3$G3A112 ## No missings, so next command has no problem bdl.g3 ← bde.g3 bdl.g3 ← ifelse(lifednd == 0, NA, bde.g3)
################## ################## ### ### ### Covariates ### ### ### ################## ##################
### HEIGHT ### Offspring Cohort variable is “A51” ### Generation 3 Cohort variable is “G3A446”
height.off ← off$A51 height.g3 ← g3$G3A446
### WEIGHT ### Offspring Cohort variable is “A50” ### Generation 3 Cohort variable is “G3A444”
weight.off ← off$A50 weight.g3 ← g3$G3A444
###————-### ### Age at exam ### ###————-### ### Get BIRTHDATE birth ← read.table(gzfile(”/work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/RootStudyConsentSet_phs000007.Framingham.v28.p10.c1.HMB-IRB-MDS/PhenotypeFiles/phs000007.v28.pht000740.v7.p10.c1.birthyr_alls.HMB-IRB-MDS.txt.gz“), header=T, sep=“\t”)
off.IDs ← off$dbGaP_Subject_ID g3.IDs ← g3$dbGaP_Subject_ID
birth.off ← subset(birth, dbGaP_Subject_ID %in% off.IDs) birth.g3 ← subset(birth, dbGaP_Subject_ID %in% g3.IDs)
### True exam dates spanned 1971 - 1975 for offspring ### True exam dates spanned 2002 - 2005 for generation 3 ### A good-enough approximation is to take the middle ### year and subtract birth year to get age at exam
### Offspring Cohort off.age ← data.frame(dbGaP_Subject_ID = birth.off$dbGaP_Subject_ID,
Age = 1973 - birth.off$birthyr)
### Framingham Cohort g3.age ← data.frame(dbGaP_Subject_ID = birth.g3$dbGaP_Subject_ID,
Age = 2004 - birth.g3$birthyr)
age ← rbind(off.age, g3.age) #################################### #################################### ### ### ### Create Phenotype Data Frames ### ### ### #################################### ####################################
offspring ← data.frame(dbGaP_Subject_ID = off$dbGaP_Subject_ID,
shareid = off$shareid,
cpd = cpd.off,
si = si.off,
sc = sc.off,
ai = ai.off,
dpw = dpw.off,
dnd = dnd.off,
bde = rep(NA, nrow(off)),
bdl = rep(NA, nrow(off)),
height = height.off,
weight = weight.off,
cohort = rep(1, nrow(off)))
generation3 ← data.frame(dbGaP_Subject_ID = g3$dbGaP_Subject_ID,
shareid = g3$shareid,
cpd = cpd.g3,
si = si.g3,
sc = sc.g3,
ai = ai.g3,
dpw = dpw.g3,
dnd = dnd.g3,
bde = bde.g3,
bdl = bdl.g3,
height = height.g3,
weight = weight.g3,
cohort = rep(2, nrow(g3)))
tmp ← rbind(offspring, generation3)
phenotypes ← merge(tmp, age, by=“dbGaP_Subject_ID”) phenotypes$age2 ← phenotypes$Age^2
################## ### ID MAPPING ### ################## ### ID mapping file from dbGaP_Subject_ID to SAMPID ID.map ← read.table(gzfile(”/work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/RootStudyConsentSet_phs000007.Framingham.v28.p10.c1.HMB-IRB-MDS/PhenotypeFiles/phs000007.v28.pht001415.v16.p10.Framingham_Sample.MULTI.txt.gz“), header=T, sep=“\t”, stringsAsFactors=F)
### Genotype files genotype.IDs ← read.table(”/work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/ChildStudyConsentSet_phs000342.Framingham.v16.p10.c1.HMB-IRB-MDS/GenotypeFiles/phg000006.v9.FHS_SHARe_Affy500K.genotype-calls-matrixfmt.c1/subject_level_PLINK_sets/FHS_SHARe_Affy500K_subjects_c1.fam“, header=F) names(genotype.IDs) ← c(“famid”, “SAMPID”, “patid”, “matid”, “sex”, “phenotype”)
length(which(genotype.IDs$SAMPID %in% ID.map$SAMPID)) ## 6954
x ← merge(genotype.IDs, ID.map, by=“SAMPID”, all.x=T) x ← x[,c(1:5,7)]
### There are many duplicates, which I can remove because I only want ### a single entry for every SAMPID-SUBJID x ← x[which(!duplicated(x$dbGaP_Subject_ID)),]
almost.final ← merge(phenotypes, x, by=“dbGaP_Subject_ID”, all.x=T)
#################################################### ### Bring in genetic PCs and ancestry categories ### ####################################################
PCs ← read.table(”/work/KellerLab/Zhen/FRAMINGHAM/PCA/FRAMINGHAM_pcs_and_ancestries.txt“, header=T) names(PCs)[2] ← c(“SAMPID”)
final ← merge(almost.final, PCs, by=“SAMPID”, all.x=T)
phenotypes.ped ← subset(final, ancestry==“EUR”,
select=c("famid", "SAMPID", "patid",
"matid", "sex", "cpd", "ai", "si",
"sc", "dpw", "dnd", "bde", "bdl"))
phenotypes.ped[is.na(phenotypes.ped)] ← “x”
covariates.ped ← subset(final, ancestry==“EUR”,
select=c("famid", "SAMPID", "patid",
"matid", "sex", "Age", "age2", "sc",
"height", "weight", "cohort",
"PC1", "PC2", "PC3", "PC4", "PC5",
"PC6", "PC7", "PC8", "PC9", "PC10"))
covariates.ped[is.na(covariates.ped)] ← “x”
write.table(phenotypes.ped, file=“Framingham.EUR.phenotypes.ped”, quote=F,
col.names=T, row.names=F)
write.table(covariates.ped, file=“Framingham.EUR.covariates.ped”, quote=F,
col.names=T, row.names=F)
We used the Affy 500K genotypes found here: /work/KellerLab/GSCAN/dbGaP/Framingham/PhenoGenotypeFiles/ChildStudyConsentSet_phs000342.Framingham.v16.p10.c1.HMB-IRB-MDS/GenotypeFiles/phg000006.v9.FHS_SHARe_Affy500K.genotype-calls-matrixfmt.c1/subject_level_PLINK_sets/FHS_SHARe_Affy500K_subjects_c1.[bed|bim|fam]
phenotypes ← read.table(”/work/KellerLab/GSCAN/dbGaP/ARIC/PhenoGenotypeFiles/ChildStudyConsentSet_phs000090.ARIC_RootStudy.v3.p1.c1.HMB-IRB/PhenotypeFiles/phs000090.v3.pht000114.v2.p1.c1.GENEVA_ARIC_Subject_Phenotypes.HMB-IRB.txt.gz“,header=T,sep=“\t”,stringsAsFactors=F)
phenotypes ← subset(phenotypes, select=c(“geneva_id”, “racegrp”, “gender”,”v1age01“, “anta01”,”anta04“, “drnkr01”, “hom29”, 'hom35', “hom32”, 'cigt01','evrsmk01', 'dtia90','dtia96', 'dtia97','dtia98', 'cursmk01','forsmk01'))
### rename phenotypes to be readable names(phenotypes)[c(1,2,3,4,5,6)] ← c(“geneva_id”, “race”, “sex”, “age”, “height”, “weight”)
### To connect sample ids to geneva ids, take SAMPID (the ID used in the genotype fam file, and SUBJID (aka geneva_id) ) id_map ← read.table(gzfile(”/work/KellerLab/GSCAN/dbGaP/ARIC/PhenoGenotypeFiles/ChildStudyConsentSet_phs000090.ARIC_RootStudy.v3.p1.c1.HMB-IRB/GenotypeFiles/phg000035.v1.ARIC_GEI.genotype-qc.MULTI/geno-qc/samp-subj-mapping.csv.gz“),header=T,sep=”,“,stringsAsFactors=F)[,c(2,1)] names(id_map) ← c(“SAMPID”, “geneva_id”) phenotypes ← merge(phenotypes, id_map, by=“geneva_id”, all.x=TRUE)
### import genotype data to get family info fam_data ← read.table(”/work/KellerLab/GSCAN/dbGaP/ARIC/PhenoGenotypeFiles/ChildStudyConsentSet_phs000090.ARIC_RootStudy.v3.p1.c1.HMB-IRB/GenotypeFiles/phg000035.v1.ARIC_GEI.genotype-calls-matrixfmt.c1.GRU.update1/Genotypes_with_flagged_chromosomal_abnormalities_zeroed_out/ARIC_PLINK_flagged_chromosomal_abnormalities_zeroed_out.fam“, col.names = c(“fam_id”, “SAMPID”, “patid”, “matid”, “sex”, “dummy”))
### Replace 0's with “x” for rvTest preferred formatting fam_data$patid ← fam_data$matid ← fam_data$fam_id[fam_data$fam_id == 0] ← “x”
######################################## ###—- Derive GSCAN phenotypes —–### ########################################
### DRINKER VERSUS NON-DRINKER ### ARIC variable name is “drnkr01”. ### Combination of “Do you presently drink alcoholic beverages?” and “Have you ever consumed alcoholic beverages?” ### Response option for both questions are “yes” or “no”, which are turned into the options below. ### Response options: ### 1 = Current Drinker ### 2 = Former Drinker ### 3 = Never Drinker ### 4 = Unknown ### ### Descriptives: ### table(phenotypes$drnkr01) ### 1 2 3 4 ### 7257 2309 3153 6 ### ### To obtain GSCAN “DND” collapse across Former and Never Drinkers ### and make “Non-Drinkers”. Current Drinkers will be made “Drinkers”
dnd ← phenotypes$drnkr01 dnd[dnd == 1] ← “Current Drinker” dnd[dnd == 2 | dnd == 3] ← 1 dnd[dnd == “Current Drinker”] ← 2 dnd[dnd == 4 | is.na(dnd)] ← “x”
### AGE OF INITIATION OF SMOKING ### ### ARIC variable name is “hom29”. ### “How old were you when you first started regular cigarette smoking?” ### Response option is an integer value. ### ### Descriptives: ### ### > table(phenotypes$hom29) ### 0 1 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ### 19 1 2 10 11 15 22 32 65 44 154 187 302 659 941 715 ### 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 ### 1219 567 703 447 275 129 88 247 56 45 59 22 100 8 34 8 ### 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 ### 17 51 9 9 10 8 35 3 11 5 5 16 3 4 2 3 ### 50 51 52 55 57 59 60 62 ### 7 1 2 1 1 2 2 1 ### ### > summary(phenotypes$hom29) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.00 16.00 18.00 18.77 20.00 62.00 5377 ### ai ← phenotypes$hom29 ### remove ages older than 35 and younger than 10 ai[ai > 35 | ai < 10 | is.na(ai)] ← “x”
### CIGARETTES PER DAY ### ARIC variable name is “hom35” ### “On the average of the entire time you smoked, how many cigarettes did you usually smoke per day?” ### Response option is integer, or “0” for <1 cigarette per day ### ### > table(phenotypes$hom35) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 46 54 118 187 133 254 146 84 89 17 990 23 106 30 12 520 ### 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 ### 31 30 62 5 2559 1 5 10 5 193 9 3 7 6 771 3 ### 32 33 34 35 36 37 38 40 42 43 45 50 51 54 55 58 ### 1 1 1 44 3 1 1 572 1 3 14 72 1 1 3 1 ### 60 65 70 75 80 86 90 99 ### 100 1 4 2 10 1 1 3 ### > summary(phenotypes$hom35) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.00 10.00 20.00 19.67 24.00 99.00 5420 ### Responses are binned in accordance with the GSCAN Analysis Plan. cpd ← phenotypes$hom35 cpd[cpd ⇐ 5 & cpd >= 1] ← 1 cpd[cpd ⇐ 15 & cpd >= 6] ← 2 cpd[cpd ⇐ 25 & cpd >= 16] ← 3 cpd[cpd ⇐ 35 & cpd >= 26] ← 4 cpd[cpd >= 36 & cpd ⇐ 60] ← 5 cpd[cpd > 60 | is.na(cpd)] ← “x”
### DRINKS PER WEEK ### ARIC variable names are “dtia96”, “dtia97”, and “dtia98” ### “dtia96” - “How many glasses of wine do you usualy have per week? (4oz. glasses; round down).” ### “dtia97” - “How many bottles of cans or beer do you usualy have per week? (12oz. bottles or cans; round down).” ### “dtia98” - “How many drinks of hard liquor do you usualy have per week? (4oz. glasses; round down).” ### Response option for all three is integer. ### ### Descriptives: ### ### >table(phenotypes$dtia96) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 14 15 16 ### 5226 844 461 255 147 90 75 50 27 3 34 1 15 28 9 2 ### 17 18 20 21 25 28 30 32 33 35 40 ### 1 3 7 5 1 2 3 1 1 1 1 ### ### >summary(phenotypes$dtia96) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 0.000 0.000 0.868 1.000 40.000 5478 ### ### >table(phenotypes$dtia97) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 4356 674 528 312 214 107 297 93 76 10 97 3 186 4 40 28 ### 16 18 19 20 21 22 23 24 25 28 30 32 33 35 36 40 ### 5 28 1 36 19 1 1 89 8 8 12 2 1 10 8 6 ### 42 45 48 49 50 56 60 63 70 72 80 92 ### 13 2 6 1 3 2 4 1 1 2 1 1 ### ### >summary(phenotypes$dtia97) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 0.000 0.000 2.609 2.000 92.000 5474 ### ### >table(phenotypes$dtia98) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 14 15 16 ### 5226 844 461 255 147 90 75 50 27 3 34 1 15 28 9 2 ### 17 18 20 21 25 28 30 32 33 35 40 ### 1 3 7 5 1 2 3 1 1 1 1 ### ### >summary(phenotypes$dtia96) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 0.000 0.000 0.868 1.000 40.000 5478 ### ### >table(phenotypes$dtia97) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 4356 674 528 312 214 107 297 93 76 10 97 3 186 4 40 28 ### 16 18 19 20 21 22 23 24 25 28 30 32 33 35 36 40 ### 5 28 1 36 19 1 1 89 8 8 12 2 1 10 8 6 ### 42 45 48 49 50 56 60 63 70 72 80 92 ### 13 2 6 1 3 2 4 1 1 2 1 1 ### ### >summary(phenotypes$dtia97) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 0.000 0.000 2.609 2.000 92.000 5474 ### ### >table(phenotypes$dtia98) ### 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ### 4387 735 551 295 239 190 148 138 70 11 151 14 57 3 103 33 ### 16 17 18 20 21 24 25 26 27 28 30 32 33 34 35 36 ### 9 9 7 36 30 1 10 1 1 9 6 2 1 2 4 1 ### 39 40 44 45 47 48 50 51 52 54 55 56 63 64 75 77 ### 1 7 1 1 1 2 5 1 1 1 1 3 1 2 1 2 ### 90 99 ### 1 2 ### ### >summary(phenotypes$dtia98) ### Min. 1st Qu. Median Mean 3rd Qu. Max. NA's ### 0.000 0.000 0.000 2.227 2.000 99.000 5483
wine ← phenotypes$dtia96 # 1 drink = 4oz beer ← phenotypes$dtia97 # 1 drink = 12oz spirits ← phenotypes$dtia98 # 1 drink = 1.5oz
### muliply wine by 4 and divide wine by 5 to normalize to standard drink of 5 oz for wine ### Combine all alcohol types, left-anchor at 1, and log dpw ← log1) former.smoker ← subset(phenotypes, select=c(“forsmk01”)) N ← nrow(phenotypes) sc ← rep(NA, N) for(i in 1:N){
if(is.na(current.smoker[i,1]) | is.na(former.smoker[i,1])){
sc[i] <- NA
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 0){
sc[i] <- NA
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 1){
sc[i] <- 1 ### former smokers are coded as 1
}
else if (current.smoker[i,1] == 1 & former.smoker[i,1] == 0){
sc[i] <- 2 ### current smokers are coded as 2 }
} sc[is.na(sc)] ← “x”
### Create dataframe with our new GSCAN variables N ← nrow(phenotypes) NAs ← rep(“x”, N) gscan.phenotypes ← data.frame(SAMPID = phenotypes$SAMPID, famid = NAs,geneva_id = phenotypes$geneva_id,patid = NAs,matid = NAs,sex = ifelse(phenotypes$sex == “M”, 1, 2),cpd = cpd,ai = ai,si = si,sc = sc,dnd = dnd,dpw = dpw,age = phenotypes$age,age2 = phenotypes$age^2,height = phenotypes$height,weight = phenotypes$weight,currentformersmoker = sc)
### Merge in the SAMPID, which is used in the genotype files gscan.phenotypes ← merge(gscan.phenotypes, id_map, by=“geneva_id”, all.x=TRUE)
### Reorder phenotype file to make pedigree file consistent with genotype IDs gscan.phenotypes ← gscan.phenotypes[c(17,2,16,3:15)] colnames(gscan.phenotypes) [2] ← “SAMPID”
### Read in PCs and add to pedigree file, then write out to a phenotype and covariate file ### [ here read in PCs and merge into phenotype file (probably by the SAMPID) ] pcs ← read.table(”/rc_scratch/hayo0753/aric/aric_ancestry_and_pcs“, head=TRUE, stringsAsFactors=F) colnames(pcs) [1] ← “SAMPID” pcs ← merge(pcs, gscan.phenotypes, by=“SAMPID”, all.x=TRUE)
############# PRELIMINARY #####################
### Write to file [NOTE TO HANNAH: will have to be changed once we ### have PCs and ancestry groups identified. PCs will have to be read ### in like with read.table() and we'll have to subset the dataset ### into European and African ancestry, and then write out one ### phenotype and covariate file per ancestry group.
### EUROPEANS phenotypes.EUR.ped ← subset(pcs, ancestry == “EUR”,select=c(“famid”,”SAMPID“,”patid“,”matid“, “sex”,”cpd“, “ai”,”si“, “sc”, “dnd”,”dpw“)) write.table(phenotypes.EUR.ped, file=“ARIC.EUR.phenotypes1.ped”, quote=F, col.names=T, row.names=F, sep=“\t”)
covariates.EUR.ped ← subset(pcs, ancestry == “EUR”, select=c(“famid”,”SAMPID“,”patid“,”matid“, “sex”,”age“, “age2”, “height”, “weight”, “currentformersmoker”,”PC1“, “PC2”, “PC3”, “PC4”, “PC5”, “PC6”, “PC7”, “PC8”,”PC9“, “PC10”)) write.table(covariates.EUR.ped, file=“ARIC.EUR.covariates1.ped”, quote=F, col.names=T, row.names=F, sep=“\t”)
phenotypes.AFR.ped ← subset(pcs, ancestry == “AFR”,select=c(“famid”,”SAMPID“,”patid“,”matid“, “sex”,”cpd“, “ai”,”si“, “sc”, “dnd”,”dpw“)) write.table(phenotypes.AFR.ped, file=“ARIC.AFR.phenotypes1.ped”, quote=F, col.names=T, row.names=F, sep=“\t”)
covariates.AFR.ped ← subset(pcs, ancestry == “AFR”, select=c(“famid”,”SAMPID“,”patid“,”matid“, “sex”,”age“, “age2”, “height”, “weight”, “currentformersmoker”, “PC1”, “PC2”, “PC3”, “PC4”, “PC5”, “PC6”, “PC7”, “PC8”,”PC9“, “PC10”)) write.table(covariates.AFR.ped, file=“ARIC.AFR.covariates1.ped”, quote=F, col.names=T, row.names=F, sep=“\t”)
### Must remove duplicates from all files, use UNIX ### sort ARIC.EUR.phenotypes1.ped | uniq > ARIC.EUR.phenotypes.ped ### sort ARIC.EUR.covariates.ped | uniq > ARIC.EUR.covariates.ped ### sort ARIC.AFR.phenotypes1.ped | uniq > ARIC.AFR.phenotypes.ped ### sort ARIC.AFR.covariates.ped | uniq > ARIC.AFR.covariates.ped
(Hannah/Joyce to update following Framingham as a guide)
Description of phenotypes can be found here: file_mesa_phenotypes_-_final.pdf
options(stringsAsFactors=F)
### eMERGE is broken into different consent classes. We can conduct analyses on hmb, hmb-gso-nic, and
emerge.hmb ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c1.HMB/PhenotypeFiles/phs000360.v3.pht003255.v2.p1.c1.MergedSet_Subject_Phenotypes.HMB.txt.gz“, header=TRUE, sep=“\t”, stringsAsFactors=F) emerge.hmb.genos ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c1.HMB/GenotypeFiles/matrix/c1.HMB/eMerge_660_11212012_c1.fam“, header=FALSE, sep=“\t”, stringsAsFactors=F)
emerge.hmb.gso.nic ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c3.HM-B-GSO-NIC/PhenotypeFiles/phs000360.v3.pht003255.v2.p1.c3.MergedSet_Subject_Phenotypes.HM-B-GSO-NIC.txt.gz“, header=TRUE, sep=“\t”, stringsAsFactors=F) emerge.hmb.gso.nic.genos ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c3.HM-B-GSO-NIC/GenotypeFiles/matrix/c3.HM-B-GSO-NIC/eMerge_660_11212012_c3.fam“, header=FALSE, sep=“\t”, stringsAsFactors=F)
emerge.hmb.gso ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c4.HMB-GSO/PhenotypeFiles/phs000360.v3.pht003255.v2.p1.c4.MergedSet_Subject_Phenotypes.HMB-GSO.txt.gz“, header=TRUE, sep=“\t”, stringsAsFactors=F) emerge.hmb.gso.genos ← read.table(”/work/KellerLab/dbGaP/eMERGE-MergedSet/PhenoGenotypeFiles/RootStudyConsentSet_phs000360.eMERGE_MergedSet.v3.p1.c4.HMB-GSO/GenotypeFiles/matrix/c4.HMB-GSO/eMerge_660_11212012_c4.fam“, header=FALSE, sep=“\t”, stringsAsFactors=F)
### Merge all files above according to SUBJID, which is used in the ### genotype files.
emerge ← merge(emerge.hmb, emerge.hmb.gso, all=T) emerge ← merge(d, emerge.hmb.gso.nic, all=T)
### 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 ### 2217 1736 3457 9635
si ← emerge$SMOKING_STATUS si[si == “C67147” | si == “C67148”] ← 2 si[si == “C65108”] ← 1 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 == “C67147”] ← 2 sc[sc == “C67148”] ← 1 sc[sc != 1 & sc != 2] ← NA
### eMERGE age variable is tricky because there is no obvious age at ### assessment. We will use their “DECADE_BIRTH” as a terrible ### approximation. ### 1=1900-1919; 2=1920-1929, 3=1930-1939; 4=1940-1949; 5=1950-1959; 6=Unknown ### ### Descriptives: ### ### table(emerge$DECADE_BIRTH) ### ### . 1 2 3 4 5 6 7 8 9 99 ### 6 612 2667 3533 4439 3127 1291 761 490 10 109 birthyear ← emerge$DECADE_BIRTH birthyear[birthyear == “99”] ← NA birthyear[birthyear == ”.“] ← NA
### SEX sex ← emerge$SEX sex[sex == “C46109”] ← 1 sex[sex == “C46110”] ← 2
### 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 = "x",
matid = "x",
sex = sex,
si = si,
sc = sc)
phenotypes[is.na(phenotypes)] ← “x”
write.table(phenotypes,
"/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/eMERGE/GSCAN_eMERGE_phenotypes.ped",
row.names=F,
quote = F,
sep="\t")
covariates ← data.frame(fid = emerge$SUBJID,
iid = emerge$SUBJID,
patid = "x",
matid = "x",
sex = sex,
birthyear = birthyear)
covariates[is.na(covariates)] ← “x”
write.table(covariates,
"/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/eMERGE/GSCAN_eMERGE_covariates.ped",
row.names=F,
quote = F,
sep="\t")
### Date: Feb 13 2017 ### Author: Scott Vrieze
options(stringsAsFactors=F)
### Load in dataset ###
ninds ← read.table(”/work/KellerLab/GSCAN/dbGaP/CIDR_StrokeGenetics/PhenoGenotypeFiles/RootStudyConsentSet_phs000615.CIDR_International_StrokeGenetics.v1.p1.c3.GRU-NPU/PhenotypeFiles/phs000615.v1.pht003307.v1.p1.c3.IschemicStroke_Subject_Phenotypes.GRU-NPU.txt.gz“,
header = T, sep="\t")
### 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] ← “XXX” ninds$dbGaP_Subject_ID ← row.names(ninds) ninds$smokingStatus ← NULL names(ninds) ← c(names(ninds)[2:17], “smokingStatus”, “dbGaP_Subject_ID”)
### subset the only variables needed pheno ← subset(ninds, select=c(“subject_id”, “smokingStatus”, “age”, “gender”, “AFFECTION_STATUS”))
###————————————————————### ### SMOKING INITIATION ### ###————————————————————###
### NINDS variable is “smokingStatus” ### Variables are “NEVER”, “FORMER”, “CURRENT”, “NO”, “UNKNOWN”, and blank. ### ### ’CURRENT' defined as cigarette smoking within last 30 days, ### 'FORMER' defined as ### more than 100 cigarettes in one's lifetime ### but no smoking within the last 30 ### days; ’NEVER' defined as ### less than 100 cigarettes smoked in one's lifetime. ### table(pheno$smokingStatus) ### ### CURRENT FORMER NEVER NO UNKNOWN ### 63 727 1137 2397 4 251
si ← pheno$smokingStatus si[si == “CURRENT” | si == “FORMER”] ← 2 si[si == “NEVER” | si == “NO”] ← 1 si[si != “1” & si != “2”] ← NA si ← as.numeric(si)
### table(si) ### ### 1 2 ### 2401 1864
###————————————————————-### ### SMOKING CESSATION ### ###————————————————————-###
### NINDS variable is “smokingStatus” ### Variables are “NEVER”, “FORMER”, “CURRENT”, “NO”, “UNKNOWN”, and blank. ### ### ’CURRENT' defined as cigarette smoking within last 30 days, ### 'FORMER' defined as more than 100 cigarettes in one's lifetime but ### no smoking within the last 30 days; ’NEVER' defined as less than ### 100 cigarettes smoked in one's lifetime. ### table(pheno$smokingStatus) ### ### CURRENT FORMER NEVER NO UNKNOWN ### 63 727 1137 2397 4 251
sc ← pheno$smokingStatus sc[sc == “CURRENT”] ← 2 sc[sc == “FORMER”] ← 1 sc[sc != “1” & sc != “2”] ← NA sc ← as.numeric(sc)
### table(pheno$V5) ### ### 1 2 ### 1137 727
###——————————### ### GENDER ### ###——————————###
### NINDS variable is “gender” ### Variables are “F” and “M” ### table(pheno$gender) ### ### F M ### 2627 1952
sex ← pheno$gender sex[sex == “M”] ← 1 sex[sex == “F”] ← 2 sex ← as.numeric(sex)
###—————-### ### Write to files ### ###—————-###
### This study uses the “subject_id” ID field in the genotype files, ### so use that here, instead of the dbGaP_Subject_ID phenotypes ← data.frame(fid = pheno$subject_id,
iid = pheno$subject_id,
patid = "x",
matid = "x",
sex = sex,
si = si,
sc = sc)
write.table(phenotypes,
"/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/Stroke/GSCAN_STROKE_phenotypes.ped",
row.names=F,
quote = F,
sep="\t")
covariates ← data.frame(fid = pheno$subject_id,
iid = pheno$subject_id,
patid = "x",
matid = "x",
sex = sex,
age = pheno$age,
age2 = pheno$age^2,
AFFECTION_STATUS = pheno$AFFECTION_STATUS)
write.table(covariates,
"/work/KellerLab/vrieze/GSCAN/GWAS/summary_stats_generated_internally/Stroke/GSCAN_STROKE_covariates.ped",
row.names=F,
quote = F,
sep="\t") ### save table
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
hmb ← read.table(”/work/KellerLab/GSCAN/dbGaP/JHS/PhenoGenotypeFiles/ChildStudyConsentSet_phs000499.JHS.v3.p1.c3.HMB-IRB/PhenotypeFiles/phs000499.v3.pht003209.v1.p1.c3.JHS_CARe_Subject_Phenotypes.HMB-IRB.txt.gz“, header=T, sep=“\t”, stringsAsFactors=F) npu ← read.table(”/work/KellerLab/GSCAN/dbGaP/JHS/PhenoGenotypeFiles/ChildStudyConsentSet_phs000499.JHS.v3.p1.c1.HMB-IRB-NPU/PhenotypeFiles/phs000499.v3.pht003209.v1.p1.c1.JHS_CARe_Subject_Phenotypes.HMB-IRB-NPU.txt.gz“, header=T, sep=“\t”, stringsAsFactors=F) phens ← rbind(hmb, npu)
hmb_geno ← read.table(”/work/KellerLab/GSCAN/dbGaP/JHS/PhenoGenotypeFiles/ChildStudyConsentSet_phs000499.JHS.v3.p1.c3.HMB-IRB/GenotypeFiles/phg000103.v1.JHS.genotype-calls-matrixfmt.c3.HMB-IRB/JHS_PLINK_illu_GRU.fam“, header=F, stringsAsFactors=F) npu_geno ← read.table(”/work/KellerLab/GSCAN/dbGaP/JHS/PhenoGenotypeFiles/ChildStudyConsentSet_phs000499.JHS.v3.p1.c1.HMB-IRB-NPU/GenotypeFiles/phg000103.v1.JHS.genotype-calls-matrixfmt.c1.HMB-IRB-NPU/JHS_PLINK_illu_NCU.fam“, header=F, stringsAsFactor=F) geno ← rbind(hmb_geno, npu_geno)
jhs_map ← read.table(”/work/KellerLab/GSCAN/dbGaP/JHS/PhenoGenotypeFiles/ChildStudyConsentSet_phs000499.JHS.v3.p1.c3.HMB-IRB/GenotypeFiles/phg000104.v1.JHS.sample-info.MULTI/JHS_subject_sample_map_8_10_2010.txt“, header=T, sep=“\t”, stringsAsFactors=F) jhs_pcs ← read.table(”/rc_scratch/hayo0753/JHS_files/JHS_ancestry_and_pcs.txt“, header=T, stringsAsFactors=F)
names(jhs_pcs) [2] ← “SUBJID” names(geno) ← c(“SUBJID”, “iid”, “patid”, “matid”, “sex”, “phen”) filtered_geno ← geno[(geno$SUBJID %in% phens$SUBJID),] filtered_phens ← phens[(phens$SUBJID %in% geno$SUBJID),]
mapped_geno ← merge(filtered_phens, filtered_geno, by=“SUBJID”, all.x=T) mapped_geno_pcs ← merge(mapped_geno, jhs_pcs, by=“SUBJID”, all.x=T)
covariates ← mapped_geno_pcs[c(1,69,70:72,5,75:84)] names(covariates) [1] ← “fid” names(covariates) [2] ← “iid” names(covariates) [6] ← “age”
###——————–### ### Smoking initiation ### ###——————–### ### ### JHSCare variables are “current_smoker_baseline” and “former_smoking_baseline” ### current_smoker_baseline = Current smoking status at first participant visit ### former_smoker_baseline = Former smoking status at first participant visit ### ### Response options are ### 0 - No ### 1 - Yes ### ### table(si) ### ### 0 1 x ### 1206 537 9
current.smoker ← subset(mapped_geno_pcs, select=c(“current_smoker_baseline”, “SUBJID”)) former.smoker ← subset(mapped_geno_pcs, select=c(“former_smoker_baseline”, “SUBJID”)) N ← nrow(mapped_geno_pcs) si ← rep(NA, N) for(i in 1:N){
if(is.na(current.smoker[i,1]) | is.na(former.smoker[i,1])){
si[i] <- "x"
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 0){
si[i] <- "1"
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 1){
si[i] <- 2 ### former smokers are coded as 2
}
else if (current.smoker[i,1] == 1 & former.smoker[i,1] == 0){
si[i] <- 2 ### current smokers are coded as 2 }
} mapped_geno_pcs_phen ← cbind(si, mapped_geno_pcs)
###——————–### ### Smoking Cessation ### ###——————–### ### ### JHSCare variables are “current_smoker_baseline” and “former_smoking_baseline” ### current_smoker_baseline = Current smoking status at first participant visit ### former_smoker_baseline = Former smoking status at first participant visit ### ### Response options are ### 0 - No ### 1 - Yes ### ### table(sc) ### ### 0 1 x ### 1206 537 9
current.smoker ← subset(mapped_geno_pcs, select=c(“current_smoker_baseline”, “SUBJID”)) former.smoker ← subset(mapped_geno_pcs, select=c(“former_smoker_baseline”,”SUBJID“)) N ← nrow(mapped_geno_pcs) sc ← rep(NA, N) for(i in 1:N){
if(is.na(current.smoker[i,1]) | is.na(former.smoker[i,1])){
sc[i] <- "x"
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 0){
sc[i] <- "x"
}
else if (current.smoker[i,1] == 0 & former.smoker[i,1] == 1){
sc[i] <- 1 ### former smokers are coded as 1
}
else if (current.smoker[i,1] == 1 & former.smoker[i,1] == 0){
sc[i] <- 2 ### current smokers are coded as 2 }
} phenotypes ← cbind(sc, mapped_geno_pcs_phen)
phenotypes ← phenotypes[c(3,71,72:74,1:2)] names(phenotypes)[1] ← “fid”
write.table(covariates, “AFR.JHS.covariates.ped”, row.names=F, quote=F) write.table(phenotypes, “AFR.JHS.phenotypes.ped”, row.names=F, quote=F)
QC parameters that we chose: MAF > 0.01
SNP callrate > 0.95
Missingness per individual > 0.95
HWE = 0.05 / number of markers but greater than 5e-8
## To update the strand builds: http://www.well.ox.ac.uk/~wrayner/strand/ ## Check strands against latest 1000G: http://www.well.ox.ac.uk/~wrayner/tools/ perl HRC-1000G-check-bim.pl \
## Phasing using shapeit (ARIC used as an example) shapeit -B ARIC_b37_filtered-updated-chr${1} -M /rc_scratch/meli7712/dbGAP/1000GP_Phase3/genetic_map_chr${1}_combined_b37.txt -O phased/ARIC_b37_filtered-updated-chr${1}.phased -T 48
## To convert the shapeit output into vcf shapeit -convert –input-haps mesa-chr${1}.phased \
## Imputation Minimac3-omp3 –haps mesa-chr${1}.phased.vcf \