search

bernatgel / CopyNumberPlots

R package to create plots representing copy number data using karyoploteR
6 stars 2 forks source link

Problem loading CNVkit Copy number calls #18

Open moldach opened 4 years ago

moldach commented 4 years ago

I'm having a problem withloadCopyNumberCallsCNVkit(). Which output file from CNVkit is it supposed to accept?

Trying the .cns file doesn't produce an error (or output) after ~25 minutes.

s2 <- CopyNumberPlots::loadCopyNumberCallsCNVkit("data/470.sorted.dedupped.cns")
Reading data from data/470.sorted.dedupped.cns

Using the .cnr file loads with seqnames, ranges and strand; however, there is no lrr or baf in the metadata and it looks weird:

R> s2 <- CopyNumberPlots::loadCopyNumberCallsCNVkit("data/470.sorted.dedupped.cnr")
Reading data from data/470.sorted.dedupped.cnr
The column identified as Segment Value is: log2
R> s2
GRanges object with 20058 ranges and 4 metadata columns:
        seqnames            ranges strand |
           <Rle>         <IRanges>  <Rle> |
      1        X            0-4999      * |
      2        X         4999-9999      * |
      3        X        9999-14999      * |
      4        X       14999-19998      * |
      5        X       19998-24998      * |
    ...      ...               ...    ... .
  20054       IV 17483829-17488829      * |
  20055       IV 17488829-17493829      * |
  20056    MtDNA            0-4598      * |
  20057    MtDNA         4598-9196      * |
  20058    MtDNA        9196-13794      * |
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     gene
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              <character>
      1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     -,cTel7X,CHROMOSOME_X,WBGene00196667,WBGene00201059,cTel7X.2,cTel7X.3,WBGene00008349,CE7X_3.1,sjj_cTel7X.1,WBsf971962,WBGene00008351,cTel7X.1:wp184,cTel7X.1,mv_cTel7X.1,cenix:243-e8,CE7X_3,WBsf973834,WBsf1003667,WBsf042175,CE7X_3.1:wp233
      2                                                                                                                                                                                                                                                                                                                                                                                                                      -,CHROMOSOME_X,WBGene00008349,CE7X_3.1,CE7X_3,WBsf042175,CE7X_3.1:wp233,WBsf899406,WBsf898320,WBsf590792,WBsf585098,WBsf552856,WBsf576981,WBsf641471,WBsf579259,WBsf1022030,WBsf578650,WBsf580657,WBsf597709,WBsf260486,WBsf640303,WBsf260485,WBsf260487,WBsf042242,CE7X_3.2:wp240,WBGene00008350,CE7X_3.2,WBsf1003668,WBsf593198,WBsf1022031,WBsf975402,WBsf1003669,WBsf1022032,WBsf595772,WBsf593906,WBsf1003670,WBsf1022033,WBsf598391,WBsf975882,WBsf1003671,WBsf1022034
      3                                                                                                                                                                                                                                                                                                                                                                                                                                                              -,CHROMOSOME_X,CE7X_3,WBsf971563,WBsf710434,WBsf594153,WBsf600890,WBsf044263,WBsf598958,WBsf565814,WBsf583518,WBsf558211,WBsf710435,WBsf688496,WBsf1022035,WBsf971994,WBsf585542,pfd-2,WBsf1003672,WBsf586175,WBsf1003673,WBsf1022036,WBsf600695,WBsf599045,WBsf1003674,Predicted_NDNAX1_CE_5322,WBsf591715,WBsf1022037,WBsf1003675,WBsf581968,WBGene00235377,CE7X_3.4,WBGene00201656,CE7X_3.3,Predicted_PALTTAA2_CE_5323,WBsf972773
      4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                -,CHROMOSOME_X,CE7X_3,WBGene00235377,CE7X_3.4,CEXL,WBGene00021342,Y35H6.3,VY35H6BL,cenix:243-e8_1,Aff_CTEL7X.A,mv_cTel7X.1_2,sjj_cTel7X.1_2,VY35H6BL.2:wp240,WBsf973734,WBsf263870
      5                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         -,CHROMOSOME_X,WBGene00021342,Y35H6.3,VY35H6BL,VY35H6BL.2:wp240,WBsf1022038,WBsf1003676,WBsf973754,WBsf207141,WBsf209190,Predicted_PALTTAA2_CE_5331,WBsf578558,WBsf575402,WBsf560922,Predicted_PALTTAA2_CE_5332,WBsf594384,VY35H6BL.1:wp240,WBsf042145,WBsf113130,WBsf209358,WBsf116194,WBsf205132,WBsf973815,Predicted_PALTTAA2_CE_5338,WBsf973842,WBsf114466,WBsf208703
    ...                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               ...
  20054                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            CHROMOSOME_IV,-,WBsf493472,4R79,WBsf718268,WBsf485680,WBsf684809,WBsf684810,4R79.1,sjj_4R79.1,cenix:83-a2,WBsf624310,WBsf250836,WBsf627805,WBGene00003525,4R79.1a,WBsf627942,WBsf962665,4R79.1b,4R79.1:wp111,4R79.1c,Aff_4R79.1,mv_4R79.1,WBsf628090,WBsf627482,WBsf627806,WBsf627943,WBsf250837,WBsf021561,WBsf472096
  20055                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               CHROMOSOME_IV,-,WBsf493472,4R79,WBGene00003525,4R79.1a,4R79.1:wp111,4R79.1c,mv_4R79.1,WBsf472096,WBsf962666,WBGene00174183,4R79.5,WBsf484140,WBsf487193,WBsf652431,WBsf486082,WBsf684811,WBsf684812,WBsf898532,WBsf899195,cTel79B,WBGene00198969,WBGene00202392,cTel79B.1,cTel79B.2
  20056 -,WBGene00014450,MTCE.1,CHROMOSOME_MtDNA,RNAz-514325,WBGene00014451,MTCE.2,MTCE.3,WBGene00010957,mv_CAA38153,smd_ND6,WBsf979332,WBsf979461,MTCE.4,WBGene00010958,smd_NL4L,WBsf979359,WBGene00014452,MTCE.5,WBGene00014453,MTCE.6,WBGene00014454,MTCE.7,smd_S-RRNA_[18S_RNA],WBsf965597,WBsf965598,RNAz-514327,WBsf965589,WBsf965592,RNAz-514328,WBsf979280,WBsf965593,WBsf979288,RNAz-514329,WBsf979412,WBGene00014455,MTCE.8,WBGene00014456,MTCE.9,WBGene00014457,MTCE.10,WBGene00010959,MTCE.11,smd_ND1,sjj2_MTCE.11,MTCE.12,WBGene00010960,mv_CAA38154,smd_ATP6,WBsf264632,WBsf979535,WBGene00014458,MTCE.13,WBGene00014459,MTCE.14,RNAz-514331,WBGene00014460,MTCE.15,WBGene00010961,MTCE.16,mv_CAA38155,smd_ND2,sjj2_MTCE.16,WBsf979406,WBGene00014461,MTCE.17,WBGene00014462,MTCE.18,WBGene00014463,MTCE.19,WBGene00014464,MTCE.20,sjj2_MTCE.21,MTCE.21,WBGene00000829,mv_CAA38156,smd_CYTB
  20057                                                                                                                                                                                                                                                                                                                                                                                                                                   CHROMOSOME_MtDNA,-,sjj2_MTCE.21,MTCE.21,WBGene00000829,mv_CAA38156,smd_CYTB,WBsf965594,WBsf264633,WBsf264634,WBsf264635,WBsf264636,WBsf979442,WBsf979321,WBsf979271,WBsf979510,WBGene00014465,MTCE.22,MTCE.23,WBGene00010962,mv_CAA38157,smd_COIII,sjj2_MTCE.23,WBsf965595,ctc-3,WBsf979265,WBGene00014466,MTCE.24,MTCE.25,WBGene00010963,mv_CAA38158,nduo-4,sjj2_MTCE.25,smd_ND4,WBsf979516,WBGene00010964,MTCE.26,mv_CAA38159,sjj2_MTCE.26,WBsf899443,smd_COI
  20058                                                                                                                                                                                                                                               CHROMOSOME_MtDNA,-,WBGene00010964,MTCE.26,mv_CAA38159,smd_COI,WBsf965596,WBsf979259,WBGene00014467,MTCE.27,WBGene00014468,MTCE.28,sjj2_MTCE.31,WBGene00014469,MTCE.29,WBGene00014470,MTCE.30,MTCE.31,WBGene00010965,mv_CAA38160,smd_COII,WBsf979293,WBsf965582,ctc-2,WBsf965583,WBsf979353,WBsf979267,WBsf979281,WBGene00014471,MTCE.32,WBGene00014472,MTCE.33,RNAz-514332,smd_L-RRNA,WBsf965584,WBsf965585,WBsf965586,WBsf965587,RNAz-514333,WBsf965588,WBsf965590,MTCE.34,WBGene00010966,mv_CAA38161,WBsf965591,smd_ND3,WBsf979299,MTCE.35,WBGene00010967,WBsf264637,smd_ND5,sjj2_MTCE.35,WBsf899316,WBsf979298,WBGene00014473,MTCE.36,WBsf980914
            depth segment.value    weight
        <numeric>     <numeric> <numeric>
      1   87.6325      0.137598    0.9996
      2  161.3030      1.023410    0.9998
      3  165.3880      1.031640    0.9998
      4   75.2577     -0.126536    0.9996
      5  130.7140      0.792624    0.9998
    ...       ...           ...       ...
  20054   78.4524    -0.1261020  0.999800
  20055   75.5212    -0.0306035  0.999800
  20056 1412.9700     4.0109000  0.919416
  20057 1654.0200     4.2390800  0.919416
  20058 1049.0400     3.5812400  0.919416
  -------
  seqinfo: 7 sequences from an unspecified genome; no seqlengths

Any idea what's the problem?

miriammagallon commented 4 years ago

Hi @moldach ,

loadCopyNumberCallsCNVkit() is able to load all cnvkit data (.cnr and .cns).

Regarding the long time trying to load the .cns file, I think that it is because of the long list of genes annotated by region. Try to run CNVkit again without annotation and then to reload the .cns file.

Secondly, the .cnr file is correctly loaded. The baf column is given in the .cns file. The lrr column is named as segment.value. This is a mislabeling problem that we have to solve in this function.

So, to sum up, you should run again CNVkit without annotation and reload de .cns file with the function. In addition, you must remember that the lrr column is named as segment.value.

moldach commented 4 years ago

Thanks for the tip @bernatgel.

Now that it's working I have another question about labelling. From the Vignette:

library(CopyNumberPlots)
s1.calls.file <- system.file("extdata", "S1.segments.txt", package = "CopyNumberPlots", mustWork = TRUE)
s1.calls <- loadCopyNumberCalls(s1.calls.file)

kp <- plotKaryotype(chromosomes="chr1")
plotCopyNumberCalls(kp, s1.calls, r0=0, r1=0.10)

Looks like:

s1.calls
GRanges object with 13 ranges and 2 metadata columns:
     seqnames              ranges strand |        cn       loh
        <Rle>           <IRanges>  <Rle> | <integer> <integer>
   1     chr1          1-60000000      * |         1         1
   2     chr1   60000001-60000999      * |         2         0
   3     chr1   60001000-62990000      * |         0         1
   4     chr1   62990001-62999999      * |         2         0
   5     chr1  63000000-121500000      * |         1         1
  ..      ...                 ...    ... .       ...       ...
   9     chr1 189600352-220352872      * |         3         0
  10     chr1 220352873-220352971      * |         2         0
  11     chr1 220352972-234920000      * |         5         0
  12     chr1 234920001-234999999      * |         2         0
  13     chr1 235000000-249250621      * |         3         0
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

When I load my data:

s2.calls <- CopyNumberPlots::loadCopyNumberCallsCNVkit("data/470.sorted.dedupped.cns")

I see the following:

s2.calls
GRanges object with 12 ranges and 7 metadata columns:
     seqnames           ranges strand |        gene segment.value     depth    probes     weight
        <Rle>        <IRanges>  <Rle> | <character>     <numeric> <numeric> <integer>  <numeric>
   1        X       0-17718942      * |           -            NA   65.0631      3500  3426.3700
   2        I       0-15072434      * |           -            NA   69.0837      2771  2808.6200
   3        V        0-2349907      * |           -            NA   69.1714       453   454.3200
   4        V  2349907-2564899      * |           -      -1.01244   30.6558        43    41.8325
   5        V 2564899-15299400      * |           -            NA   67.0854      2510  2458.3900
  ..      ...              ...    ... .         ...           ...       ...       ...        ...
   8       II       0-15279421      * |           -            NA   65.2653      2899 2897.14000
   9      III        0-7404355      * |           -            NA   69.0932      1364 1390.20000
  10      III  7404355-7454351      * |           -       2.06886  290.4790        10    9.72848
  11      III 7454351-13783801      * |           -            NA   65.6806      1193 1186.87000
  12       IV       0-17493829      * |           -            NA   69.0006      3037 3074.20000
         ci_lo     ci_hi
     <numeric> <numeric>
   1        NA        NA
   2        NA        NA
   3        NA        NA
   4  -1.03255 -0.993146
   5        NA        NA
  ..       ...       ...
   8        NA        NA
   9        NA        NA
  10   1.92209   2.14312
  11        NA        NA
  12        NA        NA
  -------
  seqinfo: 6 sequences from an unspecified genome; no seqlengths

What corresponds to cn and loh in s1.calls?

Trying this on my data produces an error:

s2.calls <- CopyNumberPlots::loadCopyNumberCallsCNVkit("data/470.sorted.dedupped.cns")

custom.genome <- toGRanges(data.frame(chr=c("I", "II", "III", "IV", "V", "X", "MtDNA"), start=c(1, 1, 1, 1, 1, 1, 1), end=c(15072434, 15279421, 13783801, 17493829, 20924180, 17718942, 13794)))

kp <- plotKaryotype(genome = custom.genome)
plotCopyNumberCalls(kp, s2.calls, r0=0, r1=0.10)
Error in plotCopyNumberCalls(kp, s2.calls, r0 = 0, r1 = 0.1) : 
  The cn.calls object does not have a column cn. No copy number data is available

Not sure if there is a lack of copy number data (if so why?) or a mislabelling issue.

miriammagallon commented 4 years ago

Hi @moldach

The .cns file you have loaded doesn't contain the cn data.

Firstly, you have to take into account CNVkit gives the results in different files. In .cnr file, you will find the lrr information. Here is where you will find the mislabelling after loading this type of file with loadCopyNumberCallsCNVkit(). We are considering developing a function to specifically load the CNVkit's .cnr files.

Another important thing you have to know concerning CNVkit is that you obtain two different .cns files. The first .cns file is the one you get after running the segment function. In this file, you won't have the cn information. You will have the segment.value information. This means that when you load .cns files with loadCopyNumberCallsCNVkit() the column named as segment.value will be correct.

The last thing you have to consider is that if you want to obtain the cn information in your .cns file you have to run the function call of CNVkit. That's why you don't have this information loaded now.

So, you can load all type of CNVkit data with loadCopyNumberCallsCNVkit(). The issue here is that you have to know what is contained in each CNVkit file.

moldach commented 4 years ago

Hi @miriammagallon.

I ran the following in cnvkit: cnvkit.py batch N2_trim_bwaMEM_sort_dedupped.bam -n -m wgs -f /scratch/moldach/data/references/c_elegans.PRJNA13758.WS265.genomic.fa

which produced the following files

base) mtg@mtg-ThinkPad-P53:~/projects/data/celegans$ ll
total 7668364
drwxrwxr-x 6 mtg mtg       4096 May 20 16:33 ./
drwxrwxr-x 3 mtg mtg       4096 May  1 15:43 ../
-rw-rw-r-- 1 mtg mtg         31 May 20 16:33 470.sorted.dedupped.antitargetcoverage.cnn
-rw-r----- 1 mtg mtg 5352369367 May 19 23:47 470.sorted.dedupped.bam
-rw-rw-r-- 1 mtg mtg     307672 May 20 08:47 470.sorted.dedupped.bam.bai
-rw-rw-r-- 1 mtg mtg         54 May 20 16:33 470.sorted.dedupped.bintest.cns
-rw-rw-r-- 1 mtg mtg        595 May 20 16:33 470.sorted.dedupped.call.cns
-rw-rw-r-- 1 mtg mtg     922610 May 20 16:33 470.sorted.dedupped.cnr
-rw-rw-r-- 1 mtg mtg        609 May 20 16:33 470.sorted.dedupped.cns
-rw-rw-r-- 1 mtg mtg     741769 May 20 16:33 470.sorted.dedupped.targetcoverage.cnn
-rw-rw-r-- 1 mtg mtg          0 May 20 08:46 c_elegans.PRJEB28388.WS274.genomic.antitarget.bed
-rw-rw-r-- 1 mtg mtg        136 May 20 08:46 c_elegans.PRJEB28388.WS274.genomic.bed
-rw-rw-r-- 1 mtg mtg     618515 May 20 08:46 c_elegans.PRJEB28388.WS274.genomic.target.bed
-rw-rw-r-- 1 mtg mtg          0 May 20 16:30 c_elegans.PRJNA13758.WS265.genomic.antitarget.bed
-rw-rw-r-- 1 mtg mtg         96 May 20 16:30 c_elegans.PRJNA13758.WS265.genomic.bed
-rw-rw-r-- 1 mtg mtg  101957874 May 20 08:55 c_elegans.PRJNA13758.WS265.genomic.fa
-rw-rw-r-- 1 mtg mtg        181 May 20 08:57 c_elegans.PRJNA13758.WS265.genomic.fa.fai
-rw-rw-r-- 1 mtg mtg     426633 May 20 16:30 c_elegans.PRJNA13758.WS265.genomic.target.bed
drwxrwxr-x 2 mtg mtg       4096 May 19 18:58 MADDOG/
-rwxr-x--- 1 mtg mtg 2391634258 May 19 22:52 N2_trim_bwaMEM_sort_dedupped.bam*
drwxrws--- 2 mtg mtg       4096 May  1 16:43 PRJEB28388/
-rw-rw-r-- 1 mtg mtg     939345 May 20 13:46 pyenv.log
-rw-rw-r-- 1 mtg mtg     738490 May 20 16:30 reference.cnn
-rw-rw-r-- 1 mtg mtg     738490 May 20 08:57 reference.cnn.1
-rw-rw-r-- 1 mtg mtg     738490 May 20 16:25 reference.cnn.2
drwxrwxr-x 6 mtg mtg       4096 May 20 16:18 venv/
drwxrwx--- 3 mtg mtg       4096 May  1 19:47 WS265_wormbase/

Loading the .cns sample with the copy numbers looks like this:

> s2.calls <- CopyNumberPlots::loadCopyNumberCallsCNVkit("data/470.sorted.dedupped.call.cns")
> s2.calls
GRanges object with 10 ranges and 7 metadata columns:
     seqnames           ranges strand |        gene segment.value
        <Rle>        <IRanges>  <Rle> | <character>     <numeric>
   1        X       0-17718942      * |           -            NA
   2        I       0-15072434      * |           -            NA
   3        V        0-2349907      * |           -            NA
   4        V  2349907-2564899      * |           -      -1.54065
   5        V 2564899-20924180      * |           -            NA
   6       II       0-15279421      * |           -            NA
   7      III        0-7404355      * |           -            NA
   8      III  7404355-7454351      * |           -       1.54065
   9      III 7454351-13783801      * |           -            NA
  10       IV       0-17493829      * |           -            NA
            cn     depth     p_ttest    probes     weight
     <integer> <numeric>   <numeric> <integer>  <numeric>
   1         2   65.0631 5.05544e-24      3500 3426.37000
   2         2   69.0837 9.63092e-10      2771 2808.62000
   3         2   69.1714 9.22086e-01       453  454.32000
   4         0   30.6558 2.05389e-34        43   41.83250
   5         2   70.3195 1.32293e-03      3551 3506.15000
   6         2   65.2653 3.02552e-01      2899 2897.14000
   7         2   69.0932 8.75159e-07      1364 1390.20000
   8         6  290.4790 1.18285e-06        10    9.72848
   9         2   65.6806 2.89004e-01      1193 1186.87000
  10         2   69.0006 5.70572e-02      3037 3074.20000
  -------
  seqinfo: 6 sequences from an unspecified genome; no seqlengths

Okay so I can see that this file has cn and segment.value but [according to section 5.4 of the vignette] (https://bioconductor.org/packages/release/bioc/vignettes/CopyNumberPlots/inst/doc/CopyNumberPlots.html) you need

a GRanges object with a at least one column of: “cn” for integer copy number calls “segment.value” for non-integer segment regional values * “loh” a logical for loss-of-heterozygosity

Where is loh found? If it's in another file how am I supposed to merge them? For example, the following doesn't work on Granges:

dplyr::inner_join(s2.calls,s3.calls) Error in UseMethod("inner_join") : no applicable method for 'inner_join' applied to an object of class "c('GRanges', 'GenomicRanges', 'GRanges_OR_NULL', 'Ranges', 'GenomicRanges_OR_missing', 'GenomicRanges_OR_GenomicRangesList', 'GenomicRanges_OR_GRangesList', 'List', 'Vector', 'list_OR_List', 'Annotated', 'vector_OR_Vector')"

miriammagallon commented 4 years ago

Hi @moldach,

CNVkit has not generated the LOH column but you don't need this information to use CopyNumberPlots. You can already plot your data with CopyNumberPlots functions. If you want to get the LOH information you should take a look at CNVkit documentation or write to the developers.