1 VCFについて

• https://samtools.github.io/hts-specs/VCFv4.2.pdf

2 ライブラリ

library(readr); library(dplyr); library(tidyr)

3 VCFファイル読み込み

• readr::read_delimで *.vcf.gzをそのまま読み込める.
  • comment = "##"でヘッダーを読み込まない
• 1列目の列名は#CHROMになっているので変更する.
• サンプルデータはRNA-seqでの変異検出のデータ
  • transdecoderで予測されたCDSに対してマッピングしている(対応するアミノ酸配列がある)
  • 6サンプルに対して変異コールしている.
  • bcftools callは-mオプション(デフォルト)で実行

in_vcf <- "~/pub/sampledata/flt.vcf.gz"
vcf <- readr::read_delim(in_vcf, delim = "\t", col_names = T, comment = "##")
names(vcf)[1] <- "CHROM"

4 FILTER

• PASS:
• MinDP
• MinAB
• MinMQ
• Qual
• SnpGap
• EndDistBias
• GapWin
• StrandBias

# PASSしたサイトのみ取り出す.
fltvcf <- filter(vcf, FILTER == "PASS")

5 QUAL

• quality scoreの値

hist(vcf$QUAL, breaks = 50)

6 INFO

indel <- vcf %>% filter(grepl("^INDEL", INFO))
snp <- vcf %>% filter(!grepl("^INDEL", INFO))

# INFO値を行持ちデータに置き換えた後,列持ちデータに変換
info.dat <- strsplit(snp$INFO, ";") %>% # INFO列分割
  setNames(., 1:nrow(snp)) %>% # リストにユニークIDを付与
  stack() %>% # 行持ちデータに変換
  mutate(id = snp$CHROM[.$ind]) %>%　# 対応するCHROM列を追加
  {cbind(.$ind, .$id, do.call(rbind, strsplit(.$values, "=")))} %>% # keyとvalueに分割
  data.frame(., stringsAsFactors = F) %>% # data.frameに変換
  setNames(., c("ind","CHROM", "k", "v")) %>%
  tidyr::spread(., "k", "v") %>% # 列持ちデータに変換
  arrange(as.integer(ind))

# 各値のhist 
par(mfrow = c(3,3))
hist(as.numeric(info.dat$AF1), main = "AF1"); 
hist(as.numeric(info.dat$FQ), main = "FQ");
hist(as.numeric(info.dat$DP), main = "DP"); 
hist(as.numeric(info.dat$HWE), main = "HWE");
hist(as.numeric(info.dat$MQ), main = "MQ"); 
hist(as.numeric(info.dat$MQSB), main = "MQSB")
hist(as.numeric(info.dat$RPB), main = "RPB")
hist(as.numeric(info.dat$SGB), main = "SGB")
hist(as.numeric(info.dat$VDB), main = "VDB")

# DP4を分割する
dp4 <- do.call(rbind, strsplit(info.dat$DP4, ",")) %>% 
  data.frame() %>%
  setNames(., c("rf","rr","af","ar")) %>%
  bind_cols(info.dat[1:2], .) 

7 GT(Genotype)とPL(Phred-scaled genotype likelihoods)

• GT: 0/0(REFのホモ), 0/1(REFとALTのヘテロ), 1/1(ALTのホモ)
• PL: 191,0,255のようにカンマ区切りの値は, 各GTのもっともらしさ. 低いほど信頼度高
• REF=C, ALT=Gの場合
  • GTの値が0/0(CC), 0/1(CG), 1/1(GG)に対応する.
  • PLの値(104,0,139)はそれぞれCC,CG,GGに対応しており,CG(0/1)が採用.
• REF=T, ALT=C,G (ALTの候補が2つある)の場合
  • PL(255,255,0,255,255,255)はそれぞれ, TT,TC,CC,TG,CG,GGに対応しているこの場合GTはCC(1/1)が採用

# サンプル別の列ラベル作成
smp <- names(snp)[10:ncol(snp)]
gt.smp <- paste("GT", smp, sep = ".")
pl.smp <- paste("PL", smp, sep = ".")

# GT,PL
gtpl <- lapply(snp[10:ncol(snp)], function(x){
    data.frame(do.call(rbind, strsplit(x, ":")))
  }) %>%
  {
    cbind(do.call(cbind, lapply(., "[", 1)), do.call(cbind, lapply(., "[", 2)))
  } %>%
  setNames(., c(gt.smp, pl.smp)) %>%
  bind_cols(snp[1:5], .)

# 複数候補があるサイトを取り出す
mgt <- gtpl[sapply(strsplit(gtpl$ALT, ","), length) > 1,]

8 変異アノテーション.

# fastaファイル読み込み.
in_orf <- "~/pub/sampledata/transdecoder/longest_orfs.cds"
orf <- Biostrings::readDNAStringSet(in_cds)
names(orf) <- sapply(strsplit(names(orf), " "), "[", 1)

snv.anno <- function(vcf, cds){
  
  # VCFに記述されているCHROMに対応するCDSを取り出し, 翻訳する.
  cds <- cds[match(vcf$CHROM, names(cds))]
  trcds <- Biostrings::translate(cds[match(vcf$CHROM, names(cds))], no.init.codon = T)
  # length(cds); length(trcds)
  # all.equal(Biostrings::width(cds), Biostrings::width(trcds)*3)
  
  # mismatch positionを見つける関数 
  seqdiff <- function(seq1, seq2){
    seq <- strsplit(c(seq1, seq2), split = '')
    mismatches <- which(seq[[1]] != seq[[2]])
    return(mismatches)
  }
  
  # annotation 付加 vcfを一行づつ処理
  pos_alt <- vcf[c(1,2,4,5)] 
  snv_anno <- unlist(sapply(1:nrow(pos_alt), function(i){
    # extract vcf data ----
    idx <- as.character(unlist(pos_alt[i,]))
    id <- idx[1]; pos <- as.integer(idx[2]); ref <- idx[3]; alt <- idx[4]
    # id; pos; ref; alt;pos_alt[i,]
    
    # REF cds, ALT cds ----
    ref_cds <- toString(cds[[i]])
    alt_cds <- ref_cds
    if (substr(ref_cds, pos, pos) == ref) {
      substr(alt_cds, pos, pos) <- alt
    } else {
      cat(paste0('The "REF" nuc is different from "', ref, '".'))
    }
    
    # REF aa, ALT aa ----
    ref_pep <- toString(trcds[[i]])
    alt_pep <- toString(Biostrings::translate(Biostrings::DNAString(alt_cds),
                                             no.init.codon = T))
    mismatch_pos <- seqdiff(ref_pep, alt_pep)

    # compare REF aa, ALT aa -> synonymous or non-synonymous ----
    if (length(mismatch_pos)) {
      ref_aa <- substr(ref_pep, mismatch_pos, mismatch_pos)
      alt_aa <- substr(alt_pep, mismatch_pos, mismatch_pos)
      paste0(mismatch_pos,":", ref_aa, "->", alt_aa)
  
    } else {
      "synonymous"
    }
  }))
  
  # substitution data and aa seq for provean
  prove <- sapply(strsplit(snv_anno, ":|->"), function(x){
    if (length(x) != 1) { 
      paste(x[c(2,1,3)], collapse = "")
    } else {
      "synonymous"
    }
  })
  vtrcds <- sapply(trcds, function(x) toString(x))
  prove.seq <- ifelse(snv_anno != "synonymous", vtrcds, NA)
  
  snv_flt <- pos_alt %>% 
    mutate(aa = snv_anno, 
           provean = prove, 
           prove.seq) %>%
    mutate(nuc = paste0(POS,":",REF,"->",ALT)) %>%
    select(1:4,8,5,6,7)
  return(snv_flt)
} 

res <- snv.anno(snp, cds)

9 環境.

sessionInfo()

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] tidyr_0.8.3   dplyr_0.8.0.1 readr_1.3.1  
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.1          highr_0.7           pillar_1.3.1       
##  [4] compiler_3.5.2      plyr_1.8.4          XVector_0.22.0     
##  [7] tools_3.5.2         zlibbioc_1.28.0     digest_0.6.18      
## [10] rskodat_0.1.0       viridisLite_0.3.0   evaluate_0.13      
## [13] tibble_2.1.1        pkgconfig_2.0.2     rlang_0.3.3        
## [16] rstudioapi_0.9.0    yaml_2.2.0          parallel_3.5.2     
## [19] xfun_0.5            rsko_0.1.0          kableExtra_1.0.1   
## [22] stringr_1.4.0       httr_1.4.0          knitr_1.21         
## [25] xml2_1.2.0          Biostrings_2.50.2   S4Vectors_0.20.1   
## [28] hms_0.4.2           IRanges_2.16.0      webshot_0.5.1      
## [31] stats4_3.5.2        tidyselect_0.2.5    glue_1.3.1         
## [34] R6_2.4.0            rmarkdown_1.11      purrr_0.3.2        
## [37] magrittr_1.5        scales_1.0.0        htmltools_0.3.6    
## [40] BiocGenerics_0.28.0 assertthat_0.2.1    rvest_0.3.2        
## [43] colorspace_1.4-1    stringi_1.4.3       munsell_0.5.0      
## [46] crayon_1.3.4