Need R help (fdrtool, PatternCNV)

Question

[eyb]

I am writing here since I never heard back from the authors of PatternCNV. It is a library for discovering CNV variations in exome sequencing. http://bioinformaticstools.mayo.edu/research/patterncnv/. There is a user manual. Among other things there is a protocol for a CNV detection in germline samples.

The problem I am facing is that I can not calculate FDR. PatternCNV function uses fdrtool to calculate FDR and takes null model as input. Since I am only analysing ~90kb long section of 1st chromosome my null model looks like this:

> germline.null_model
$type
[1] "autosome"*

$Chr
 [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"  "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16"
[17] "chr17" "chr18" "chr19" "chr20" "chr21" "chr22"

$location
 [1]  0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

$scale
 [1] 0.5283454        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[14]        NA        NA        NA        NA        NA        NA        NA        NA        NA

$SD
 [1] 0.7972599        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[14]        NA        NA        NA        NA        NA        NA        NA        NA        NA

So when I try to calculate fdr I get the following error:

> germline.FDR_res <- patCNV.estimate.FDR(session_info=germline.session_info, cnv_res=germline.cnv_res, null_model = germline.null_model)

Error in if (!is.numeric(scale) | scale <= 0) { :
  missing value where TRUE/FALSE needed

I've tried to remove all NA values from null model class and now it looks like this:

> germline_2.null_model
$type
$type[[1]]
[1] "autosome"

$Chr
$Chr[[1]]
[1] "chr1"

$scale
$scale[[1]]
[1] 0.5283454

$SD
$SD[[1]]
[1] 0.5283454

$location
[1] 1

However running patCNV FDR function still gives me error. But another one:

germline.FDR_res <- patCNV.estimate.FDR(session_info=germline.session_info, cnv_res=germline.cnv_res, null_model = germline_2.null_model)

Error in patCNV.lap.pval(q = sel_vec, location = null_model$location[j], :
   input scale must be positive number

Can anyone suggest a fix please?

[Giovanni M Dall'Olio]

I am not familiar with patterncnv at all. Is there a way you can extract the p-values from germline.cnv_res? Note that you removed the NAs from the null model, but not from germline.cnv_res.

[eyb]

In my understanding germline.cnv_res does not contain neither NA nor p-values values. I think this class is for CNV probability per exon

[Charles Plessy]

There is something strange in your example output: you should not see things such as $Chr[[1]]. If this is not a cut-paste accident, it may reflect a problem.

> germline_2.null_model
$type
$type[[1]]
[1] "autosome"

$Chr
$Chr[[1]]
[1] "chr1"

In your case, you seem to want to remove all but the first chromosome. See in the following example one way to do it:

> test <- list(type = "test", Chr = c("chr1", "chr2"), scale = c(123, NA))
> test
$type
[1] "test"

$Chr
[1] "chr1" "chr2"

$scale
[1] 123  NA

> lapply(test, function(x) x[1])
$type
[1] "test"

$Chr
[1] "chr1"

$scale
[1] 123

(Note that this example does not apply to more complex cases where it is not just the first chromosome that is to be kept).

Answer 1

[oghzzang]

Hello. I had same problem and fixed that problem.

Go " PatternCNV_1.0/Rlib/functions/patCNV.export.CNV.tables.R "

And fix "for(j in 1:22)"

example

if(null_model$type=='autosome')

{

for (k in 1:N_sample)
{

  for(j in 22:22) # Chromosome_name (your case 1:1)**

  {

I know only manual editing .. Sorry ㅜ.ㅜ