NormalizeQuantiles

For julia 0.4, 0.5, 0.6 see: https://github.com/oheil/NormalizeQuantiles.jl/tree/backport-0.6

Package NormalizeQuantiles implements quantile normalization

qn = normalizeQuantiles(array);

and provides a function to calculate sample ranks

(r,m) = sampleRanks(array);

of a given vector or matrix.

References

• Amaratunga, D.; Cabrera, J. (2001). "Analysis of Data from Viral DNA Microchips". Journal of the American Statistical Association. 96 (456): 1161. doi:10.1198/016214501753381814
• Bolstad, B. M.; Irizarry, R. A.; Astrand, M.; Speed, T. P. (2003). "A comparison of normalization methods for high density oligonucleotide array data based on variance and bias". Bioinformatics. 19 (2): 185–193. doi:10.1093/bioinformatics/19.2.185 PMID 12538238
• Wikipedia contributors. (2018, June 12). Quantile normalization. In Wikipedia, The Free Encyclopedia. Retrieved 11:54, August 3, 2018, from https://en.wikipedia.org/w/index.php?title=Quantile_normalization

Table of Contents

• Dependencies
• Remarks
• Usage examples normalizeQuantiles
  • General usage
  • Missing Values
  • SharedArray and multicore usage examples
  • OffsetArrays
• Behaviour of function normalizeQuantiles
• Data prerequisites
• Remarks on data with missing values
• List of all exported definitions for normalizeQuantiles
• Usage examples sampleRanks
• List of all exported definitions for sampleRanks

• Julia 0.7 or above

• none

• Distributed
• SharedArrays
• Random
• Statistics

• for julia 0.4, 0.5, 0.6 see: https://github.com/oheil/NormalizeQuantiles.jl/tree/backport-0.6
• Code examples and output on this page have been used on and copied from the julia 0.7 REPL
• Last commit with julia 0.3 support: Jan 20, 2017, eb97d24ff77d470d0d121fabf83d59979ad0db36
  • git checkout eb97d24ff77d470d0d121fabf83d59979ad0db36

Pkg.add("NormalizeQuantiles");
using NormalizeQuantiles;

The following array is interpreted as a matrix with 4 rows and 3 columns:

array = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ];
qn = normalizeQuantiles(array)

	julia> qn
	4×3 Array{Float64,2}:
	 2.0  3.0  2.0
	 4.0  6.0  4.0
	 8.0  8.0  7.0
	 6.0  3.0  7.0

The columns in qn are now quantile normalized to each other.

The input array must not have dimension larger than 2.

Return type of function normalizeQuantiles is always Array{Float64,2}

If your data contain some missing values like NaN (Not a Number) or something else, they will be changed to NaN:

array = [ NaN 2.0 1.0 ; 4.0 "empty" 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ];

	julia> array
	4×3 Array{Any,2}:
	 NaN    2.0       1.0
	4.0   "empty"  6.0
	9.0  7.0       8.0
	5.0  2.0       8.0

qn = normalizeQuantiles(array)

	julia> qn
	4×3 Array{Float64,2}:
	 NaN      3.25  1.5
	   5.0  NaN     5.0
	   8.0    8.0   6.5
	   5.0    3.25  6.5

NaN is of type Float64, so there is nothing similar for Int types.

	julia> typeof(NaN)
	Float64

You can convert the result to Array{Union{Missing, Float64},2} with:

qnMissing = convert(Array{Union{Missing,Float64}},qn)

	julia> qnMissing
	4×3 Array{Union{Missing, Float64},2}:
	 NaN      3.25  1.5
	   5.0  NaN     5.0
	   8.0    8.0   6.5
	   5.0    3.25  6.5

qnMissing[isnan.(qnMissing)] = missing;

	julia> qnMissing
	4×3 Array{Union{Missing, Float64},2}:
	  missing  3.25      1.5
	 5.0        missing  5.0
	 8.0       8.0       6.5
	 5.0       3.25      6.5

Remark: restart julia now. addprocs() must be called before using NormalizeQuantiles;.

To use multiple cores on a single machine you can use the standard packages Distributed and SharedArrays:

using Distributed
addprocs();
@everywhere using SharedArrays
@everywhere using NormalizeQuantiles

sa = SharedArray{Float64}([ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ])

	julia> sa
	4×3 SharedArray{Float64,2}:
	 3.0  2.0  1.0
	 4.0  5.0  6.0
	 9.0  7.0  8.0
	 5.0  2.0  8.0

qn = normalizeQuantiles(sa)

	julia> qn
	4×3 Array{Float64,2}:
	 2.0  3.0  2.0
	 4.0  6.0  4.0
	 8.0  8.0  7.0
	 6.0  3.0  7.0

Remark: restart julia again.

For small data sets using Distributed and SharedArrays decreases performance:

using NormalizeQuantiles
la = randn((100,100));
normalizeQuantiles(la); @time normalizeQuantiles(la);

	julia> @time normalizeQuantiles(la);
	  0.003178 seconds (8.35 k allocations: 2.813 MiB)

Remark: restart julia.

using Distributed
addprocs();
@everywhere using SharedArrays
@everywhere using NormalizeQuantiles
sa = SharedArray{Float64}( randn((100,100)) );
normalizeQuantiles(sa); @time normalizeQuantiles(sa);

	julia> @time normalizeQuantiles(sa);
	  0.013014 seconds (12.10 k allocations: 432.146 KiB)

Remark: restart julia.

For larger data sets performance increases with multicore processors:

using NormalizeQuantiles
la = randn((1000,10000));
normalizeQuantiles(la); @time normalizeQuantiles(la);

	julia> @time normalizeQuantiles(la);
	  2.959431 seconds (784.18 k allocations: 2.281 GiB, 12.13% gc time)

Remark: restart julia.

using Distributed
addprocs();
@everywhere using SharedArrays
@everywhere using NormalizeQuantiles
la = randn((1000,10000));
sa = SharedArray{Float64}(la);
normalizeQuantiles(sa); @time normalizeQuantiles(sa);

	julia> @time normalizeQuantiles(sa);
	  1.030016 seconds (83.85 k allocations: 80.754 MiB, 5.77% gc time)

Using non-SharedArrays in a multicore setup is slowest:

	julia> normalizeQuantiles(la); @time normalizeQuantiles(la);
	  5.776685 seconds (294.06 k allocations: 92.532 MiB, 0.28% gc time)

Remark: with Julia 1.3.1 OffsetArrays are not supported until JuliaLang/julia#34886 is released (expected in Julia 1.5)

using NormalizeQuantiles, OffsetArrays

array = [ 3 missing 1 ; 4 5 6 ; missing 7 8 ; 5 2 8 ];
oa = OffsetArray(array,-1,-1);

julia> oa
4×3 OffsetArray(::Array{Union{Missing, Int64},2}, 0:3, 0:2) with eltype Union{Missing, Int64} with indices 0:3×0:2:
 3          missing  1
 4         5         6
  missing  7         8
 5         2         8

qn = normalizeQuantiles(oa);

The quantile normalized result is not an OffsetArray:

julia> qn
4×3 Array{Float64,2}:
   2.0      NaN        2.0
   4.0        6.5      4.0
 NaN          6.66667  6.58333
   6.66667    4.0      6.58333

After quantile normalization the sets of values of each column have the same statistical properties. This is quantile normalization without a reference column.

The function normalizeQuantiles expects an array with dimension <= 2 and always returns a matrix of same dimensions as the input matrix and of type Array{Float64,2}.

NaN values of type Float64 and any other non-numbers, like strings, are treated as random missing values and the result value will be NaN. See "Remarks on data with missing values" below.

Equal values in a column of the input matrix will have different quantile normalized values. Those different result values can't be assigned back to the proper original positions because they are indistinguishable. The mean value of the different result values are therefor put back into original positions.

Example:

julia> array = [ 1 2 ; 2  3 ; 2 5 ]
3×2 Matrix{Int64}:
 1  2
 2  3
 2  5

julia> qn = normalizeQuantiles(array)
3×2 Matrix{Float64}:
 1.5  1.5
 3.0  2.5
 3.0  3.5

In row 2 and 3 instead if 2.5 and 3.5 the mean 3.0 is the result in both rows.

To use quantile normalization your data should have the following properties:

• the distribution of values in each column should be similar
• number of values for each column should be large
• number of missing values in the data should be small and of random nature

Currently there seems to be no general agreement on how to deal with missing values during quantile normalization. Here we put any given missing value back into the sorted column at the original position before calculating the mean of the rows.

sampleRanks of a given vector calculates for each element the rank, which is the position of the element in the sorted vector.

using NormalizeQuantiles
a = [ 5.0 2.0 4.0 3.0 1.0 ];
(r,m) = sampleRanks(a);   # here only return value r is relevant, for m see below
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	 5
	 2
	 4
	 3
	 1

If you provide a matrix like

array = [ 1.0 2.0 3.0 ; 4.0 5.0 6.0 ; 7.0 8.0 9.0 ; 10.0 11.0 12.0 ]

	julia> array
	4×3 Array{Float64,2}:
	  1.0   2.0   3.0
	  4.0   5.0   6.0
	  7.0   8.0   9.0
	 10.0  11.0  12.0

ranks are calculated column wise, or in other words, array is treated as array[:]:

(r,m) = sampleRanks(array);
r

	julia> r
	12-element Array{Union{Missing, Int64},1}:
	  1
	  4
	  7
	 10
	  2
	  5
	  8
	 11
	  3
	  6
	  9
	 12

There are three optional keyword parameters tiesMethod, naIncreasesRank and resultMatrix:

(r,m) = sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true);
(r,m) = sampleRanks(a,resultMatrix=true);

Equal values in the vector are called ties. There are several methods available on how to treat ties:

• tmMin : the smallest rank for all ties (default)
• tmMax : the largest rank
• tmOrder : increasing ranks
• tmReverse : decreasing ranks
• tmRandom : the ranks are randomly distributed
• tmAverage : the average rounded to the next integer

These methods are defined and exported as

	@enum qnTiesMethods tmMin tmMax tmOrder tmReverse tmRandom tmAverage

Internally ties have increasing ranks. On these the chosen method is applied.

The next rank for the successive values after the ties is the so far highest used rank plus 1.

Examples:

a = [ 7.0 2.0 4.0 2.0 1.0 ];
(r,m) = sampleRanks(a); #which is the same as (r,m)=sampleRanks(a,tiesMethod=tmMin)
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	 4
	 2
	 3
	 2
	 1

(r,m) = sampleRanks(a,tiesMethod=tmMax);
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	 5
	 3
	 4
	 3
	 1

(r,m) = sampleRanks(a,tiesMethod=tmReverse);
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	 5
	 3
	 4
	 2
	 1

One or more missing values in the vector are never equal and remain on there position after sorting. The rank of each missing value is always missing::Missing. The default is that a missing value does not increase the rank for successive values. Giving true keyword parameter naIncreasesRank changes that behavior to increasing the rank by 1 for successive values:

a = [ 7.0 2.0 4.0 2.0 1.0 ];
a[1] = NaN;
(r,m) = sampleRanks(a);
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	  missing
	 2
	 3
	 2
	 1

(r,m) = sampleRanks(a,naIncreasesRank=true);
r

	julia> r
	5-element Array{Union{Missing, Int64},1}:
	  missing
	 3
	 4
	 3
	 2

The keyword parameter resultMatrix lets you generate a dictionary of rank indices to allow direct access to all values with a given rank. For large vectors this may have a large memory consumption therefor the default is to return an empty dictionary of type Dict{Int64,Array{Int64,N}}:

a = [ 7.0 2.0 4.0 2.0 1.0 ];
(r,m) = sampleRanks(a,resultMatrix=true);
m

	julia> m
	Dict{Int64,Array{Int64,N} where N} with 4 entries:
	  4 => [1]
	  2 => [2,4]
	  3 => [3]
	  1 => [5]

haskey(m,2)   #does rank 2 exist?

	julia> haskey(m,2)
	true

a[m[2]]   #all values of rank 2

	julia> a[m[2]]
	2-element Array{Float64,1}:
	 2.0
	 2.0