CYP2C19.Rmd

```{r}
library(dplyr)

c19.pos <- read.table('pos_CYP2C19.txt', header=F) %>% rename(POS=V1)
c19.rsid <- read.table('rsid_CYP2C19.txt', header=F) %>% rename(RSID=V1) %>% cbind(.,c19.pos)

# load vcf files containing CYP2C19 region
vcf<-read.table(gzfile('CYP2C19.vcf.gz'), stringsAsFactors = FALSE)[1:3] %>% rename(CHROM=V1,POS=V2,ID=V3)

# check how many CYP2C19 SNPs were covered in phased vcf files
c19.vcf <- vcf %>% merge(.,c19.rsid, by='POS') 
```


# check which SNPs belong to one star allele in CYP2C19 using CPIC reference table
```{r}
ref <- readxl::read_excel('CYP2C19_allele_definition_table.xlsx', sheet=1) %>% .[,1:ncol(.)-1]  

c19.pos <- rbind(c('POS'), c19.pos)

ref.pos <- ref %>% 
    rbind(c19.pos$POS, .) %>% 
    `colnames<-`(.[1, ]) %>% 
    .[8:nrow(.),] %>% 
    rename(allele=POS)

cand.pos <- intersect(as.character(c19.vcf$POS), names(ref.pos))  

# include alleles that only contain SNPs in the candidate SNP list
allele1 <- ref.pos[apply(ref.pos, 1, function(x) {
  non_na <- names(ref.pos)[which(!is.na(x))][-1]
  allele_with_other <- names(ref.pos)[which(x=='M'|x=='Y'|x=='R')]
  allele_names <- setdiff(non_na, allele_with_other)
  all_present <- all(allele_names %in% cand.pos)
  return(all_present)
}), ] %>% as.data.frame

# select SNPs shown in the vcf file from the reference, because all included SNPs have rsid, change this to colnames
rsid <- ref %>% 
    rbind(c19.pos$POS, .) %>% 
    .[, c(names(.)[1], names(intersect(c19.vcf$POS, .[1,])))] %>% 
    `colnames<-`(.[6, ]) %>%
    .[8:nrow(.),] %>% 
    rename(allele=rsID)

# remove alleles that have all NAs in the listed SNPs
allele2 <- rsid %>% .[rowSums(is.na(.)) != ncol(.)-1, 1]  
allele.sel <- intersect(allele1$allele, allele2$allele)  

```

# manually call CYP2C19 from genotype data
```{r}
# read files
vcf<-read.table(gzfile('CYP2C19.vcf.gz'), stringsAsFactors = FALSE) %>% 
    filter(V2 %in% c19.vcf$POS) %>% 
    merge(c19.vcf[,c('POS','RSID')], ., by.x='POS', by.y='V2') %>% 
    t %>% 
    as.data.frame %>% 
    .[c(2,5,6,11:nrow(.)),] %>% 
    `colnames<-`(.[1, ]) %>% 
    .[-1, ]
 
# only leave genotype from vcf
vcf <- apply(vcf, 2, function(x) {do.call(rbind, strsplit(x,':'))[,1]}) %>% as.data.frame

# separate diplotype to haplotype, the star alleles are ordered from the lowest activity score to the highest
vcf <- vcf %>% mutate(hap_1='*1', hap_2='*1')

# *2 rs12769205 rs4244285 rs3758581
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs12769205, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs4244285, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*2'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs12769205, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit(vcf$rs4244285, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,2] =='1' & 
    vcf$hap_2 =='*1')]  <- '*2'

# *3 rs4986893 rs3758581 
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs4986893, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*3'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs4986893, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit  (vcf$rs3758581, '\\|'))[,2] =='0' 
    & vcf$hap_2 =='*1')] <- '*3'

# *8 rs41291556 
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs41291556, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*8'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs41291556, '\\|'))[,2] =='1' & 
    vcf$hap_2 =='*1')] <- '*8'

# *35 rs12769205 rs3758581
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs12769205, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*35'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs12769205, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,2] =='1' & 
    vcf$hap_2 =='*1')] <- '*35'

# *17 rs12248560 rs3758581
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs12248560, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*17'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs12248560, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,2] =='1' & 
    vcf$hap_2 =='*1')] <- '*17'

# *11 rs58973490 rs3758581
vcf$hap_1[which(do.call(rbind, strsplit(vcf$rs58973490, '\\|'))[,1] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,1] =='1' & 
    vcf$hap_1 =='*1')] <- '*11'
vcf$hap_2[which(do.call(rbind, strsplit(vcf$rs58973490, '\\|'))[,2] =='1' & 
    do.call(rbind, strsplit(vcf$rs3758581, '\\|'))[,2] =='1'& 
    vcf$hap_2 =='*1')] <- '*11'

# activity value of each allele 
vcf <- vcf%>% mutate(as_hap1 = ifelse(
    grepl('^\\*2$|^\\*4$|^\\*8$|^\\*35$', hap_1), 0, ifelse(
        grepl('^\\*17$', hap_1), 1.5, 1
    )), as_hap2 = ifelse(
    grepl('^\\*2$|^\\*4$|^\\*8$|^\\*35$', hap_2), 0, ifelse(
        grepl('^\\*17$', hap_2), 1.5, 1
    )),
    # activity score
    as = as_hap1 + as_hap2, diplotype = paste0(hap_1, '/', hap_2
    ))

# phenotype
vcf <- vcf %>% mutate(phenotype = ifelse(
    as>=2.5, 'Ultrarapid', ifelse(
        as==0, 'Poor', ifelse(
            as<=1.5 & as>0, 'Intermediate', 'Normal'
        )
    )))
vcf <- vcf %>% .[3:nrow(vcf),]

```