WES 데이터를 통해 해상도 높은 CNV 데이터를 예측할 수 있지만, exon targeting 중 bias, artifact, 실험 노이즈로 인해 정확한 copy number 예측이 쉽지 않다. CODEX는 다수 샘플의 WES 데이터로 부터, 적절한 데이터 정규화와 함께 CNV를 탐지하는 프로그램으로, Bioconductor 모듈 가운데 하나이다.

개요 #

CNV를 탐지하는 기존의 알고리즘은 다음과 같다.

1. Read depth를 기본 데이터로 이용
2. 암 시료, 정상시료와 같은 Case/control 데이터로부터 차이나는 부분을 선별
3. SVD와 같은 행렬 분해(Matrix decomposition) 방법으로 노이즈 제거

기존 방식과 차이나는 CODEX의 특징은 다음과 같다.

1. matched normal control이 필요 없으며, 대신 같이 시퀀싱한 다수의 샘플이 필요함.
2. 이런 discrete count 데이터(coverage)는 Poisson log-linear model(Poisson regression)이 더 알맞음.
3. 이 정규화방법은 GC 함량, 엑손 길이/캡처/증폭 효율, latent systematic artifact에 의한 비뚤림(biase)을 제거함

예제 #

HapMap WES 샘플데이터로 CNV 탐지하기 예제는 다음과 같다.

설치 #

R 3.2 이상에서 구동한다. R 설치 후,

$ R
> source("http://bioconductor.org/biocLite.R")
> biocLite("CODEX")

위 명령은 의존성있는 아래 모듈들을 함께 설치한다.

• ‘lambda.r’
• ‘futile.options’
• ‘futile.logger’
• ‘snow’
• ‘BiocParallel’
• ‘XML’
• ‘RCurl’
• ‘GenomicAlignments’
• ‘rtracklayer’
• ‘S4Vectors’
• ‘IRanges’
• ‘GenomicRanges’
• ‘XVector’
• ‘Biostrings’
• ‘BiocGenerics’
• ‘zlibbioc’
• ‘bitops’
• ‘BSgenome’
• ‘Rsamtools’
• ‘GenomeInfoDb’
• ‘BSgenome.Hsapiens.UCSC.hg19’

"WES.1KG.WUGSC" 모듈에 포함된 데이터를 분석 예제 데이터로 사용한다.

> biocLite("WES.1KG.WUGSC")

위 명령으로 바로 설치 할 수 있다. .bam 파일들이 위치한 디렉토리와 엑손 타겟 위치가 저장된 .bed 파일, 샘플 이름 정보로 부터 분석을 위한 기본 설정을 수행한다.

> library(CODEX)
> library(WES.1KG.WUGSC)
> dirPath <- system.file("extdata", package="WES.1KG.WUGSC")
> bamFile <- list.files(dirPath, pattern='*.bam$')
> bamdir <- file.path(dirPath, bamFile)
> sampname <- as.matrix(read.table(file.path(dirPath, "sampname")))
> bedFile <- file.path(dirPath, "chr22_400_to_500.bed")
> chr <- 22
> bambedObj <- getbambed(bamdir=bamdir, bedFile=bedFile,
                   sampname=sampname, projectname="CODEX_demo", chr)
> bamdir <- bambedObj$bamdir; sampname <- bambedObj$sampname
> ref <- bambedObj$ref
> projectname <- bambedObj$projectname
> chr <- bambedObj$chr

.bam 파일로 부터 depth 정보를 추출한다.

> coverageObj <- getcoverage(bambedObj, mapqthres=20)
> Y <- coverageObj$Y; readlength <- coverageObj$readlength

이어서 GC 함량과 mappability를 계산한다.

> gc <- getgc(chr, ref)
> mapp <- getmapp(chr, ref)

샘플별, 엑손별로 QC 수행

> qcObj <- qc(Y, sampname, chr, ref, mapp, gc, cov_thresh = c(20, 4000), length_thresh = c(20, 2000), mapp_thresh = 0.9, gc_thresh = c(20, 80))
> Y_qc <- qcObj$Y_qc; sampname_qc <- qcObj$sampname_qc; gc_qc <- qcObj$gc_qc
> mapp_qc <- qcObj$mapp_qc; ref_qc <- qcObj$ref_qc; qcmat <- qcObj$qcmat
> write.table(qcmat, file = paste(projectname, '_', chr, '_qcmat', '.txt', sep=''), sep='\t', quote=FALSE, row.names=FALSE)

QC결과는 디렉토리에 CODEX_demo_22_qcmat.txt 이름으로 저장된다. 파일의 내용은 대략 다음과 같다.

$ head CODEX_demo_22_qcmat.txt
chr start_bp    end_bp  pass    median_depth    pass_depth  length_kb   pass_length mapp    pass_mapp   gc  pass_gc
22  21345867    21346168    0   9   0   0.302   1   1   1   64.24   1
22  21346453    21346708    1   53.5    1   0.256   1   1   1   60.55   1
22  21347033    21347243    0   12  0   0.211   1   1   1   64.93   1
22  21347901    21348093    0   20  0   0.193   1   1   1   63.73   1
22  21348163    21348608    1   27  1   0.446   1   1   1   67.94   1
22  21348797    21349066    1   63.5    1   0.27    1   1   1   64.07   1
22  21349109    21349365    1   95.5    1   0.257   1   1   1   61.87   1
22  21349985    21350451    1   151.5   1   0.467   1   1   1   62.74   1
22  21350935    21351305    1   94.5    1   0.371   1   1   1   61.19   1

데이터 정규화를 진행한다. Poisson latent factor를 1~9 사이의 값으로 각각 피팅한다.

> normObj <- normalize(Y_qc, gc_qc, K = 1:9)
> Yhat <- normObj$Yhat; AIC <- normObj$AIC; BIC <- normObj$BIC
> RSS <- normObj$RSS; K <- normObj$K

만일 case-control 실험이라면 별도로 normal 데이터를 지정해서 정규화할 수 있다.

> normObj <- normalize2(Y_qc, gc_qc, K = 1:9, normal_index=seq(1,45,2))
> Yhat <- normObj$Yhat; AIC <- normObj$AIC; BIC <- normObj$BIC
> RSS <- normObj$RSS; K <- normObj$K

계산된 latent factor, AIC, BIC에 대한 deviance reduction plot을 PDF 형식으로 저장할 수 있다.

> choiceofK(AIC, BIC, RSS, K, filename = paste(projectname, "_", chr, "_choiceofK", ".pdf", sep = ""))

이상의 정보를 바탕으로 CNV 영역을 결정할 수 있다. germline의 경우 "integer" 모드를 사용하며, somatic의 경우, "fraction" 모드를 사용한다.

> optK = K[which.max(BIC)] 
> finalcall <- segment(Y_qc, Yhat, optK = optK, K = K, sampname_qc, ref_qc, chr, lmax = 200, mode = "integer")
> finalcall

sample_name chr  cnv   st_bp      ed_bp      length_kb st_exon ed_exon raw_cov
"NA18990"   "22" "dup" "22312814" "22326373" "13.56"   "60"    "72"    "1382"
norm_cov copy_no lratio  mBIC     pvalue
"997"    "3"     "61.83" "51.228" "1e-28"

> write.table(finalcall, file = paste(projectname, '_', chr, '_', optK, '_CODEX_frac.txt', sep=''), sep='\t', quote=FALSE, row.names=FALSE) 
> save.image(file = paste(projectname, '_', chr, '_image', '.rda', sep=''), compress='xz')

위 결과 내용의 copy_no 부분이 해당 영역의 copy number를 의미한다.

참고문헌 #

1 CODEX: a normalization and copy number variation detection method for whole exome sequencing Nucleic Acids Research 2015-01-23

Suggested Pages #

• 0.025 R
• 0.025 CSV
• 0.025 normalization
• 0.013 R (체스)
• 0.013 R (수학)
• 0.006 실수
• More suggestions...