CODEX pipeline fails: Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) : subscript contains NAs
1
0
Entering edit mode
@nickjagiella-9054
Last seen 6.1 years ago

I want to use the CODEX pipeline to analyze the copy number variants (CNV) of a sample of whole exon sequencing (WES) data. 

  1. I installed the packages into R-Studio by the following commands:
    source("https://bioconductor.org/biocLite.R")
    biocLite("CODEX")
  2. I created a BED File (PDL1_exons.bed) as following:

    chr9    5357966    5358360    uc003zjd.3_exon_0_0_chr9_5357967_r    0    -
    chr9    5361077    5361187    uc003zjd.3_exon_1_0_chr9_5361078_r    0    -
    chr9    5361757    5361888    uc003zjd.3_exon_2_0_chr9_5361758_r    0    -
    chr9    5431896    5431983    uc003zjd.3_exon_3_0_chr9_5431897_r    0    -
    chr9    5436568    5436694    uc003zjd.3_exon_4_0_chr9_5436569_r    0    -
    chr9    5439142    5439266    uc003zjd.3_exon_5_0_chr9_5439143_r    0    -
    chr9    5524535    5524629    uc003zjd.3_exon_6_0_chr9_5524536_r    0    -
    chr9    5629566    5630071    uc003zjd.3_exon_7_0_chr9_5629567_r    0    -
  3. And I executed the CODEX pipeline (here) adapted to my case with 3 samples: 
library(CODEX)
dirPath <- "/Users/njagiella/Documents/Projects/PID"
bamFile <- list.files(dirPath, pattern = "*.bam$")
bamdir <- file.path(dirPath, bamFile)
sampname <- as.matrix(bamFile)
sampname
bedFile <- file.path(dirPath, "PDL1_exons.bed")
chr <- "chr9"
bambedObj <- getbambed(bamdir = bamdir, bedFile = bedFile, sampname = sampname, projectname = "TEST_WES", chr)
bamdir <- bambedObj$bamdir; sampname <- bambedObj$sampname
ref <- bambedObj$ref; projectname <- bambedObj$projectname; chr <- bambedObj$chr
coverageObj <- getcoverage(bambedObj, mapqthres = 20)
Y <- coverageObj$Y; readlength <- coverageObj$readlength
gc <- getgc(chr, ref)
mapp <- getmapp(chr, ref)
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))

The script stops with the following output:

> library(CODEX)
> dirPath <- "/Users/njagiella/Documents/Projects/PID"
> bamFile <- list.files(dirPath, pattern = "*.bam$")
> bamdir <- file.path(dirPath, bamFile)
> sampname <- as.matrix(bamFile)
> sampname
     [,1]           
[1,] "21.sorted.bam"
[2,] "24.sorted.bam"
[3,] "28.sorted.bam"
> bedFile <- file.path(dirPath, "PDL1_exons.bed")
> chr <- "chr9"
> bambedObj <- getbambed(bamdir = bamdir, bedFile = bedFile, sampname = sampname, projectname = "TEST_WES", chr)
> bamdir <- bambedObj$bamdir; sampname <- bambedObj$sampname
> ref <- bambedObj$ref; projectname <- bambedObj$projectname; chr <- bambedObj$chr
> coverageObj <- getcoverage(bambedObj, mapqthres = 20)
Getting coverage for sample 21.sorted.bam: read length 100.
Getting coverage for sample 24.sorted.bam: read length 100.
Getting coverage for sample 28.sorted.bam: read length 100.
> Y <- coverageObj$Y; readlength <- coverageObj$readlength
> gc <- getgc(chr, ref)
> mapp <- getmapp(chr, ref)
> 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))
Excluded NA exons due to extreme coverage.
Excluded 0 exons due to extreme exonic length.
Excluded 0 exons due to extreme mappability.
Excluded 0 exons due to extreme GC content.
After taking union, excluded NA out of 8 exons in QC.
Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) : 
  subscript contains NAs 
13 stop("subscript contains NAs") 
12 NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) 
11 NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) 
10 normalizeSingleBracketSubscript(i, x, as.NSBS = TRUE) 
9 extractROWS(x, i) 
8 extractROWS(x, i) 
7 .nextMethod(x, i) 
6 eval(expr, envir, enclos) 
5 eval(call, callEnv) 
4 callNextMethod(x, i) 
3 ref[binfilter] 
2 ref[binfilter] at qc.R#32
1 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)) 
bioconductor codex error • 1.2k views
ADD COMMENT
0
Entering edit mode
yuchaoj • 0
@yuchaoj-7979
Last seen 6.1 years ago
United States

Hi Nick,

Thanks for your interest! This code has been used quite extensively so there shouldn't be any huge systematic bugs. I wonder it's due to some initialization inconsistence in your code. Also have you sorted your bed file?

Can you use save.image(file='CODEX_debug.rda') after the mapp=getmapp() step and send the file to me? I will look into it. My email is yuchaoj@wharton.upenn.edu. If it is too large, you can share that with me via Dropbox: yj235@cornell.edu.

Cheers,

Yuchao

ADD COMMENT

Login before adding your answer.

Traffic: 260 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6