As Mikhael mentioned, all ratio-based normalization methods are somewhat sensitive to low counts. The most obvious problem is that ratios involving zeroes are either undefined or yield useless scaling factors of zero. Even if we were to put those aside, we now have the problem of highly discrete ratios generated from very low counts. Taking the median or trimming by quantile for such a discrete distribution performs poorly for getting an estimate of the true scaling factor. A related problem is the fact that low counts yield highly variable ratios, which results in further degradation in performance.

If you absolutely must get normalized expression values for all genes for some reason, then one approach is to filter genes by abundance, compute normalization factors from the filtered subset, and then use the normalization factors on the entire set of genes ("transplanting" the y$samples$norm.factors back into the full set). This assumes that the same scaling biases apply across both low- and high-abundance genes, which is tolerable if you're mainly concerned about composition biases. If you do this, it is important that you do not set keep.lib.sizes=FALSE during subsetting (i.e., don't go out of your way to change the defaults). The normalization factors only make sense in the context of the library sizes with which they are calculated, so if you're doing the normalization manually, you need to make sure that the library sizes during calculation of the factors are the same as the library sizes during their application.

Basically, I have ~ 58.000 transcripts, and I just want to normalize the raw counts and transform them so that I can make comparisons (I have 2 time points and 60 samples per time point).

If you want to do differential expression comparisons, this should be done on the filtered data anyway. Because it is (i) computationally faster, (ii) improves the accuracy of the mean-variance modelling and (iii) reduces the severity of the multiple testing correction.