R (프로그래밍 언어)
Bioconductor
CODEX
#
Find similar titles
- 최초 작성자
- 최근 업데이트
Structured data
- Category
- Software
WES 데이터를 통해 해상도 높은 CNV 데이터를 예측할 수 있지만, exon targeting 중 bias, artifact, 실험 노이즈로 인해 정확한 copy number 예측이 쉽지 않다. CODEX는 다수 샘플의 WES 데이터로 부터, 적절한 데이터 정규화와 함께 CNV를 탐지하는 프로그램으로, Bioconductor 모듈 가운데 하나이다.
Table of Contents
개요 #
CNV를 탐지하는 기존의 알고리즘은 다음과 같다.
- Read depth를 기본 데이터로 이용
- 암 시료, 정상시료와 같은 Case/control 데이터로부터 차이나는 부분을 선별
- SVD와 같은 행렬 분해(Matrix decomposition) 방법으로 노이즈 제거
기존 방식과 차이나는 CODEX의 특징은 다음과 같다.
- matched normal control이 필요 없으며, 대신 같이 시퀀싱한 다수의 샘플이 필요함.
- 이런 discrete count 데이터(coverage)는 Poisson log-linear model(Poisson regression)이 더 알맞음.
- 이 정규화방법은 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
Suggested Pages #
- 0.025 R
- 0.025 CSV
- 0.025 normalization
- 0.013 R (체스)
- 0.013 R (수학)
- 0.006 실수
- More suggestions...