GSCAN dbGaP
Studies
Framingham
ID Mapping
The following snppit describes how the PLINK genotype file IDs (which use "sampid") were mapped to the dbGaP_Subject_ID contained in the available phenotype files. The snippet only contains one of Joyce Ling's phenotype files for illustration only.
This script and the results it produces can be found in /work/KellerLab/GSCAN/dbGaP/Framingham/ID_mapping.r
### 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)]
### One of Joyce's phenotype ID files (eventually all phenotypes for all available participants were merged in, but this is a good example.)
phenotypes <- read.table("OffC_Exam_1.txt", header=T, sep="\t")
xx <- merge(x, phenotypes, by="dbGaP_Subject_ID", all=TRUE)
xx <- xx[,c(2:ncol(xx),1)] ## Reorder so the dbGaP_Subject_ID is the last column
write.table(xx, file="framingham_GSCAN_phenotypesCovariates.ped", quote=FALSE, sep="\t", row.names=F)
Phenotypes
(Joyce will update this section)
Genotypes
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]
ARIC
(Hannah/Joyce to update this section following Framingham as a guide)
ID Mapping
Phenotypes
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("gen_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")
### 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".
### [Hannah will fill in question text]
### 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
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 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
wine <- phenotypes$dtia96 # 1 drink = 4oz
beer <- phenotypes$dtia97 # 1 drink = 12oz
spirits <- phenotypes$dtia98 # 1 drink = 1.5oz
### divide wine by 5 to normalize to standard drink
dpw <- log((wine*4/5 + beer + spirits) + 1) ## Combine all alcohol types, left-anchor at 1, and log
dpw[is.na(dpw)] <- "x"
### SMOKING INITIATION
### ARIC variable name is "evrsmk01"
### [hannah fill in details]
si <- phenotypes$evrsmk01
si[si == 1] <- 2
si[si == 0] <- 1
si[is.na(si)] <- "x"
### SMOKING CESSATION
### ARIC variable names are "cursmk01" and "forsmk01"
current.smoker <- subset(phenotypes, select=c("cursmk01"))
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(famid = NAs,
geneva_id = phenotypes$gen_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(2,16,3:15)]
### 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) ]
############# 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.
PCs <- read.table(xxxFile from Zhenxxx)
ancestry.IDs <- read.table(xxxFile from Zhenxxx)
### EUROPEANS
phenotypes.EUR.ped <- subset(gscan.phenotypes, ancestry == "EUR",
select=c("famid","SAMPID","patid","matid", "sex",
"cpd", "ai","si", "sc", "dnd","dpw"))
write.table(phenotypes.EUR.ped, file="ARIC.EUR.phenotypes.ped", quote=F, col.names=T, row.names=F, sep="\t")
covariates.EUR.ped <- subset(gscan.phenotypes, 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.covariates.ped", quote=F, col.names=T, row.names=F, sep="\t")
phenotypes.AFR.ped <- subset(gscan.phenotypes, ancestry == "AFR",
select=c("famid","SAMPID","patid","matid", "sex",
"cpd", "ai","si", "sc", "dnd","dpw"))
write.table(phenotypes.AFR.ped, file="ARIC.AFR.phenotypes.ped", quote=F, col.names=T, row.names=F, sep="\t")
covariates.AFR.ped <- subset(gscan.phenotypes, 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.covariates.ped", quote=F, col.names=T, row.names=F, sep="\t")
Genotypes
MESA
(Hannah/Joyce to update following Framingham as a guide)
Phenotypes
Description of phenotypes can be found here: Media:MESA phenotypes - FINAL.pdf
eMERGE
(Hannah/Joyce to update following Framingham as a guide)
Phenotypes
Description of phenotypes can be found here: Media:EMERGE.pdf
Stroke
(Hannah/Joyce to update following Framingham as a guide)
Genotype Processing
Pre-Phasing QC
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/ #!/bin/bash #SBATCH --qos=blanca-ibg #SBATCH --mem=40gb perl HRC-1000G-check-bim.pl -b ARIC_b37_filtered.bim -f ARIC_b37_filtered.frq -r 1000GP_Phase3_combined.legend -g -p EUR
## Phasing using shapeit
#!/bin/bash
#SBATCH --mem=20gb
#SBATCH --time=48:00:00
#SBATCH -o shapeit_aric_%j.out
#SBATCH -e shapeit_aric_%j.err
#SBATCH --qos blanca-ibgc1
#SBATCH --ntasks-per-node 48
#SBATCH -J shapeit_aric
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
#!/bin/bash
#SBATCH --mem=20gb
#SBATCH --time=24:00:00
#SBATCH -o shapeit_mesa_%j.out
#SBATCH -e shapeit_mesa_%j.err
#SBATCH --qos janus
#SBATCH --ntasks-per-node 12
#SBATCH -J shapeit_mesa
shapeit -convert --input-haps mesa-chr${1}.phased --output-vcf mesa-chr${1}.phased.vcf -T 12
## Imputation
#!/bin/bash
#SBATCH --mem=30gb
#SBATCH --time=72:00:00
#SBATCH -o impute_mesa_%j.out
#SBATCH -e impute_mesa_%j.err
#SBATCH --qos blanca-ibgc1
#SBATCH --ntasks-per-node 48
#SBATCH -J impute_mesa
/work/KellerLab/Zhen/bin/Minimac3/bin/Minimac3 --haps mesa-chr${1}.phased.vcf --cpus 48 --refHaps /rc_scratch/meli7712/dbGAP/references/${1}.1000g.Phase3.v5.With.Parameter.Estimates.m3vcf.gz --chr ${1} --noPhoneHome --format GT,DS,GP --allTypedSites --prefix mesa-chr${1}.phased.imputed