Question: Need R help (fdrtool, PatternCNV)
0
gravatar for eyb
4.2 years ago by
eyb180
Russian Federation
eyb180 wrote:

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?

 

ADD COMMENTlink modified 22 months ago by oghzzang30 • written 4.2 years ago by eyb180

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.

ADD REPLYlink written 4.2 years ago by Giovanni M Dall'Olio26k

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

ADD REPLYlink modified 4.2 years ago • written 4.2 years ago by eyb180

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).

ADD REPLYlink written 22 months ago by Charles Plessy2.6k
0
gravatar for oghzzang
22 months ago by
oghzzang30
oghzzang30 wrote:

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 ㅜ.ㅜ

ADD COMMENTlink modified 22 months ago • written 22 months ago by oghzzang30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1100 users visited in the last hour