Error while running DE in DESeq2
1
0
Entering edit mode
@6deab231
Last seen 3.3 years ago

Hello Michael Love

When I try to run the DE analysis I get the following error. Do you have any clue why this error arise?

> ddsDE <- DESeq(ddsSE)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing

warning: solve(): system seems singular; attempting approx solution
Error in fitBeta(ySEXP = ySEXP, xSEXP = xSEXP, nfSEXP = nfSEXP, alpha_hatSEXP = alpha_hatSEXP,  : 
  BLAS/LAPACK routine 'DGELSD' gave error code -4

> sessionInfo( )
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.1 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libmkl_rt.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=el_GR.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=el_GR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=el_GR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=el_GR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods  
[9] base     

other attached packages:
 [1] DESeqAnalysis_0.3.10        DESeq2_1.30.0               bcbioRNASeq_0.3.39         
 [4] basejump_0.13.4             SummarizedExperiment_1.20.0 Biobase_2.50.0             
 [7] GenomicRanges_1.42.0        GenomeInfoDb_1.26.2         IRanges_2.24.1             
[10] S4Vectors_0.28.1            BiocGenerics_0.36.0         MatrixGenerics_1.2.0       
[13] matrixStats_0.57.0         

loaded via a namespace (and not attached):
  [1] colorspace_2.0-0              ellipsis_0.3.1               
  [3] ggridges_0.5.2                XVector_0.30.0               
  [5] AcidPlots_0.3.0               rstudioapi_0.13              
  [7] ggrepel_0.9.0                 bit64_4.0.5                  
  [9] interactiveDisplayBase_1.28.0 AnnotationDbi_1.52.0         
 [11] fansi_0.4.1                   xml2_1.3.2                   
 [13] splines_4.0.3                 tximport_1.18.0              
 [15] goalie_0.4.11                 geneplotter_1.68.0           
 [17] knitr_1.30                    Rsamtools_2.6.0              
 [19] annotate_1.68.0               dbplyr_2.0.0                 
 [21] pheatmap_1.0.12               shiny_1.5.0                  
 [23] BiocManager_1.30.10           readr_1.4.0                  
 [25] compiler_4.0.3                httr_1.4.2                   
 [27] assertthat_0.2.1              Matrix_1.2-18                
 [29] fastmap_1.0.1                 lazyeval_0.2.2               
 [31] limma_3.46.0                  cli_2.2.0                    
 [33] later_1.1.0.1                 htmltools_0.5.0              
 [35] prettyunits_1.1.1             tools_4.0.3                  
 [37] gtable_0.3.0                  glue_1.4.2                   
 [39] GenomeInfoDbData_1.2.4        dplyr_1.0.2                  
 [41] rappdirs_0.3.1                tinytex_0.28                 
 [43] Rcpp_1.0.5                    vctrs_0.3.6                  
 [45] AcidGenomes_0.1.1             Biostrings_2.58.0            
 [47] AcidGenerics_0.4.1            rtracklayer_1.50.0           
 [49] xfun_0.19                     stringr_1.4.0                
 [51] syntactic_0.4.3               mime_0.9                     
 [53] lifecycle_0.2.0               ensembldb_2.14.0             
 [55] XML_3.99-0.5                  edgeR_3.32.0                 
 [57] AnnotationHub_2.22.0          zlibbioc_1.36.0              
 [59] scales_1.1.1                  vroom_1.3.2                  
 [61] hms_0.5.3                     promises_1.1.1               
 [63] ProtGenerics_1.22.0           AnnotationFilter_1.14.0      
 [65] RColorBrewer_1.1-2            SingleCellExperiment_1.12.0  
 [67] yaml_2.2.1                    curl_4.3                     
 [69] memoise_1.1.0                 gridExtra_2.3                
 [71] ggplot2_3.3.2                 UpSetR_1.4.0                 
 [73] biomaRt_2.46.0                stringi_1.5.3                
 [75] RSQLite_2.2.1                 genefilter_1.72.0            
 [77] BiocVersion_3.12.0            GenomicFeatures_1.42.1       
 [79] BiocParallel_1.24.1           rlang_0.4.9                  
 [81] pkgconfig_2.0.3               bitops_1.0-6                 
 [83] lattice_0.20-41               purrr_0.3.4                  
 [85] GenomicAlignments_1.26.0      cowplot_1.1.0                
 [87] bit_4.0.4                     tidyselect_1.1.0             
 [89] plyr_1.8.6                    magrittr_2.0.1               
 [91] AcidPlyr_0.1.4                R6_2.5.0                     
 [93] generics_0.1.0                DelayedArray_0.16.0          
 [95] DBI_1.1.0                     pillar_1.4.7                 
 [97] withr_2.3.0                   survival_3.2-7               
 [99] RCurl_1.98-1.2                bcbioBase_0.6.16             
[101] tibble_3.0.4                  pipette_0.4.22               
[103] crayon_1.3.4                  BiocFileCache_1.14.0         
[105] progress_1.2.2                locfit_1.5-9.4               
[107] grid_4.0.3                    data.table_1.13.4            
[109] blob_1.2.1                    digest_0.6.27                
[111] xtable_1.8-4                  AcidBase_0.2.6               
[113] httpuv_1.5.4                  openssl_1.4.3                
[115] munsell_0.5.0                 sessioninfo_1.1.1            
[117] askpass_1.1
DESeq2 • 2.0k views
ADD COMMENT •
0
Entering edit mode
@mikelove
Last seen 6 hours ago
United States

Can you post:

model.matrix(colData(dds), design(dds))

It says the system is (nearly) singular, which could mean you have a design that is close to linearly dependent.

Also I note the standard parametric dispersion fit didn't work, is this typical RNA-seq data?
ADD COMMENT •
0
Entering edit mode

Michael Love Yes they are typical RNASeq data

> design(ddsSE)
~disease_status_followup
> colData(ddsSE)
DataFrame with 172 rows and 18 columns
            Age_at_Diagnosis Age_at_Enrollment Age_at_last_FollowUP
                   <numeric>         <numeric>            <numeric>
patient100             49.58           49.5852                61.43
patient1056            54.99           73.0075                75.69
patient1059            63.83           67.4114                70.36
patient1062            71.74           71.7399                74.93
patient1076            33.94           33.9439                39.48
...                      ...               ...                  ...
patient884             64.93           68.5284                70.48
patient885             61.16           63.9097                65.86
patient886             83.21           91.4743                92.01
patient887             65.90           67.6715                67.67
patient9               48.78           51.5975                61.60
            Interval_of_last_FollowUP Survival_Time binet_stage_dg binet_stage_sampling
                            <numeric>     <numeric>      <numeric>            <numeric>
patient100                       4327          4327              1                    1
patient1056                      7558          7558              1                    1
patient1059                      2385          2385              2                    2
patient1062                      1165          1165              2                    2
patient1076                      2023          2023              2                    2
...                               ...           ...            ...                  ...
patient884                       2027          2027              1                    1
patient885                       1719          1719              1                    1
patient886                       3213          3213              1                   NA
patient887                        646           646              2                    2
patient9                         4684          4684              1                    2
            consensus_diagnosis disease_status_followup    gender       icd_o
                      <numeric>                <factor> <numeric> <character>
patient100                    1              Control            1      9823/3
patient1056                   3              Pathogenic         1      9823/3
patient1059                   1              Pathogenic         0      9823/3
patient1062                   1              Control            1      9823/3
patient1076                   1              Control            1      9823/3
...                         ...                     ...       ...         ...
patient884                    1              Pathogenic         0      9823/3
patient885                    1              Pathogenic         1      9823/3
patient886                    1              Pathogenic         1      9823/3
patient887                    1              Pathogenic         1      9823/3
patient9                      1              Pathogenic         1      9823/3
            last_followup_status patient_code progression     batch phenotype
                       <numeric>    <numeric>   <numeric> <logical> <logical>
patient100                     0          100           1        NA        NA
patient1056                    0         1056           0        NA        NA
patient1059                    0         1059           1        NA        NA
patient1062                    0         1062           1        NA        NA
patient1076                    0         1076           1        NA        NA
...                          ...          ...         ...       ...       ...
patient884                     0          884           0        NA        NA
patient885                     0          885           0        NA        NA
patient886                     0          886           0        NA        NA
patient887                     0          887           1        NA        NA
patient9                       1            9           1        NA        NA
            Donor_relapse_interval      sample
                         <numeric> <character>
patient100                       2  patient100
patient1056                     NA patient1056
patient1059                   1307 patient1059
patient1062                    959 patient1062
patient1076                     41 patient1076
...                            ...         ...
patient884                      NA  patient884
patient885                      NA  patient885
patient886                      NA  patient886
patient887                     616  patient887
patient9                      1176    patient9
> model.matrix(colData(ddsSE), design(ddsSE))
Error in terms.default(object, data = data) : 
  no terms component nor attribute

Am I doing somehting wrong?

> ddsSE$disease_status_followup
  [1] Control    Pathogenic Pathogenic Control    Control    Pathogenic Pathogenic
  [8] Pathogenic Pathogenic Pathogenic Pathogenic Control    Pathogenic Pathogenic
 [15] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control   
 [22] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control   
 [29] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control   
 [36] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control   
 [43] Control    Pathogenic Pathogenic Control    Pathogenic Pathogenic Pathogenic
 [50] Pathogenic Pathogenic Pathogenic Control    Pathogenic Pathogenic Pathogenic
 [57] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
 [64] Pathogenic Control    Pathogenic Pathogenic Pathogenic Pathogenic Control   
 [71] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control    Control   
 [78] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
 [85] Control    Pathogenic Pathogenic Pathogenic Control    Control    Control   
 [92] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
 [99] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Control    Pathogenic
[106] Control    Pathogenic Pathogenic Pathogenic Pathogenic Control    Control   
[113] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
[120] Control    Pathogenic Pathogenic Control    Pathogenic Control    Pathogenic
[127] Pathogenic Pathogenic Pathogenic Control    Control    Pathogenic Pathogenic
[134] Control    Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
[141] Pathogenic Control    Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
[148] Pathogenic Pathogenic Control    Pathogenic Pathogenic Pathogenic Pathogenic
[155] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
[162] Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic Pathogenic
[169] Pathogenic Pathogenic Pathogenic Pathogenic
Levels: Control Pathogenic
ADD REPLY •
0
Entering edit mode

I do not know what the problem is but it works on my Macbook just fine. I encountered the problem on my Linux workstation only.

ADD REPLY •
0
Entering edit mode

That is strange. The design looks fine.

One technique would be to pre-filter problemsome genes with:

keep <- rowSums(counts(dds) >= 10) >= X
dds <- dds[keep,]
dds <- DESeq(dds)

Where X could be chosen as the smallest group size, i.e. number of samples in the smaller group of the design.

ADD REPLY •
0
Entering edit mode

Michael Love It seems that I get this error in certain runs only. The problem occurs when I have many samples (eg 74 and 144 in my case). In other runs with fewer samples it works fine.

I tried what you proposed but I get the same error again. Do you have any other advice?

EDIT:

It somehow WORKS in my Macbook! I do not know what is going on. There is definitely a bug in UBUNTU 20.04. I made a fresh install of R but the issue is not resolved! I do not know if you can somehow fix it, but keep this issue in mind for Ubuntu 20.04.

keep <- rowSums(counts(dds) >= 10) >= 48
dds <- dds[keep, ]

> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

warning: solve(): system seems singular; attempting approx solution
Error in fitBeta(ySEXP = ySEXP, xSEXP = xSEXP, nfSEXP = nfSEXP, alpha_hatSEXP = alpha_hatSEXP,  : 
  BLAS/LAPACK routine 'DGELSD' gave error code -4
ADD REPLY •
0
Entering edit mode

Thanks for the update. I suppose if other users encounter this, they can bump the thread, but I haven't been able to reproduce on my end.

ADD REPLY •
0
Entering edit mode

Updating Deseq2 1.32.0 error was no longer an issue. May be an issue when using 1.30.1

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 824 users visited in the last hour
Content Search
Users
Tags
Badges
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