QuackenbushLab / yarn

YARN: Robust Multi-Tissue RNA-Seq Preprocessing and Normalization
13 stars 2 forks source link

Quantile normalization is not working #3

Open dn070017 opened 5 years ago

dn070017 commented 5 years ago

Hi, I was trying to compare the normalization result from "qsmooth" and "quantile". "qsmooth" works perfectly; however, I cannot get the correct result from quantile normalization.

Here's the code that I used:

data(skin)
skin <- normalizeTissueAware(skin, "SMTSD", "quantile")
assayData(skin)$normalizedMatrix   # Return a matrix full of NA

and the session info:

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] yarn_1.10.0         Biobase_2.44.0      BiocGenerics_0.30.0 BiocManager_1.30.4  optparse_1.6.2     

loaded via a namespace (and not attached):
  [1] colorspace_1.4-1            siggenes_1.58.0             mclust_5.4.3                XVector_0.24.0             
  [5] GenomicRanges_1.36.0        quantro_1.18.0              base64_2.0                  getopt_1.20.3              
  [9] bit64_0.9-7                 AnnotationDbi_1.46.0        xml2_1.2.0                  codetools_0.2-16           
 [13] splines_3.6.0               doParallel_1.0.14           scrime_1.3.5                Rsamtools_2.0.0            
 [17] annotate_1.62.0             HDF5Array_1.12.1            readr_1.3.1                 compiler_3.6.0             
 [21] httr_1.4.0                  assertthat_0.2.1            Matrix_1.2-17               lazyeval_0.2.2             
 [25] limma_3.40.2                prettyunits_1.0.2           tools_3.6.0                 gtable_0.3.0               
 [29] glue_1.3.1                  GenomeInfoDbData_1.2.1      dplyr_0.8.1                 doRNG_1.7.1                
 [33] Rcpp_1.0.1                  bumphunter_1.26.0           Biostrings_2.52.0           multtest_2.40.0            
 [37] gdata_2.18.0                preprocessCore_1.46.0       nlme_3.1-140                rtracklayer_1.44.0         
 [41] iterators_1.0.10            DelayedMatrixStats_1.6.0    stringr_1.4.0               rngtools_1.3.1.1           
 [45] gtools_3.8.1                XML_3.98-1.19               beanplot_1.2                edgeR_3.26.4               
 [49] zlibbioc_1.30.0             MASS_7.3-51.4               scales_1.0.0                hms_0.4.2                  
 [53] SummarizedExperiment_1.14.0 minfi_1.30.0                rhdf5_2.28.0                GEOquery_2.52.0            
 [57] RColorBrewer_1.1-2          memoise_1.1.0               ggplot2_3.1.1               downloader_0.4             
 [61] pkgmaker_0.27               biomaRt_2.40.0              reshape_0.8.8               stringi_1.4.3              
 [65] RSQLite_2.1.1               genefilter_1.66.0           S4Vectors_0.22.0            foreach_1.4.4              
 [69] GenomicFeatures_1.36.1      caTools_1.17.1.2            BiocParallel_1.18.0         bibtex_0.4.2               
 [73] GenomeInfoDb_1.20.0         rlang_0.3.4                 pkgconfig_2.0.2             matrixStats_0.54.0         
 [77] bitops_1.0-6                nor1mix_1.2-3               lattice_0.20-38             purrr_0.3.2                
 [81] Rhdf5lib_1.6.0              GenomicAlignments_1.20.0    bit_1.1-14                  tidyselect_0.2.5           
 [85] plyr_1.8.4                  magrittr_1.5                R6_2.4.0                    IRanges_2.18.0             
 [89] gplots_3.0.1.1              DelayedArray_0.10.0         DBI_1.0.0                   pillar_1.4.1               
 [93] withr_2.1.2                 survival_2.44-1.1           RCurl_1.95-4.12             tibble_2.1.2               
 [97] crayon_1.3.4                KernSmooth_2.23-15          progress_1.2.2              locfit_1.5-9.1             
[101] grid_3.6.0                  data.table_1.12.2           blob_1.1.1                  digest_0.6.19              
[105] xtable_1.8-4                tidyr_0.8.3                 illuminaio_0.26.0           openssl_1.3                
[109] stats4_3.6.0                munsell_0.5.0               registry_0.5-1              askpass_1.1                
[113] quadprog_1.5-7

I quickly go through the source code, I think maybe the problems lies in line 46 of normalizeTissueAware.R

cnts <- exprs(obj[, which(pData(obj)$our %in% i)])

Maybe it should be change to something like:

cnts <- exprs(obj[, which(select(pData(obj), groups) %in% i)])

Could you please help me take a look at this issue? Thank you very much for your time.

Best, Ping-Han