I would recommend using the raw counts per gene and look into DeSeq2 or edgeR to do the testing. These are packages used by a large number of people and give robust results.

The formula as it is written above does mot make much sense to me.

To my knowledge, "quantile normalization" works differently: if you have a matrix of genes (in rows) and samples (in columns), then each column is ordered, and the values in the rows are substituted by the row-means. Then the values in the columns are brought back to their original position. As a result, the empirical distribution of the values in all columns are identical (and, thus, all moments like mean, variance, kurtosis, ...).

This normalization works very well for whole transcriptome data from comparable biological sources (it might not be suited when the global profiles can be drastically different, like for single-cell analyses or the comparisons of different cell types).

Please check GFOLD: http://dx.doi.org/10.1093/bioinformatics/bts515

I would strongly recommend to not use the log2((X1+median1)/(X2+median2)) formula you mentioned, because this is not a statistical test, does not take into account the data properties, and does not take into account variance at all.

1. Get the raw counts of the sequencing run

2. Read the well-documented DESeq2 vignette, it is easy to use and the vignette will guide you through all the necessary steps (https://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.pdf)

Litte explanation about the quantile normalization you mentioned: This usually means that the library sized were adjusted. In RNA-seq, you have to adjust for library sizes to compensate for different technical properties of each sequencing run (some runs result in more reads than others). Usually the library size normalization is not the difficulty in RNA-seq, it is all about the noise in the low count region. 

DESeq2 normalization (sometimes referred to as RLE), edgeR normalization (trimmed mean of m-values, TMM) and quantile normalization result in practically only little differences. What every normalization does, it basically says let's adjust for the majority of the data set which is more robust than adjusting for the whole data set (outliers!).

In DESeq2, instead of altering the raw counts, size factors are calculated. When performing the analysis you can extract them using sizeFactors(dds), but you can also overwrite those sizeFactors if for example you want to adjust for certain sample properties, also with "sizeFactors(dds)