Gamma-rays emitted by extragalactic sources propagating over cosmological distances has attracted growing interests and motivated a large variety of theoretical models as well as observational studies at GeV and TeV energies. However, the mean free path of these gamma-rays is limited above a certain energy (e.g. ∼ 100 GeV) since they are attenuated due to pair-production process with extragalactic background light (EBL) photons. The resulted opacity of Universe can be reduced with gamma-ray photon mixing with a proposed hypothetical axion-like particle (ALP) in the presence of external magnetic field. This is known as the anomalous transparency of Universe towards TeV gamma-rays. Whereas at GeV energies, the photon-ALP mixing can also give rise to significant observable effects on the spectral energy distributions (SEDs) of gamma-rays, where the significance of such effects depends on the strength of photons coupling to ALPs (gγγ ) as well as the ambient magnetic fields. Therefore, additional modulations shown on the observed SEDs in comparison to conventional physics are likely to result from the photon-ALP mixing effect. The significance of the irregularities on the SEDs can then be used to probe the existence of this hypothetical particle. In this work, we use the energy spectra of 20 extragalactic gamma-ray sources of AGN recorded by Fermi -LAT in a 10-year observation. The propagation of the gamma-rays is divided approximately into three distinct regions, where different magnetic field profiles are configured and used for photon-ALP conversions. For the region in the vicinity of gamma-ray sources, we model the magnetic field with a large-scale homogeneous field (characterized with free parameters B for its field strength and s for its spatial scale) as little is known about its true nature. When photon/ALP beam continues to propagate in the intergalactic space, we take into account EBL absorption effect and neglect photon-ALP conversion here since the field strength is too weak to induce any prominent modulation on the final observed SEDs. Finally, as the beam enters Milky Way, the present Galactic magnetic field contributes another possible conversion region. After fully propagating the photon/ALP beam, we fit the observational spectral data points of each source to two models with and without the introduction of photon-ALP mixing effect respectively, namely, the H0 (without photon-ALP mixing) and H1 hypotheses (with photon-ALP mixing). Motivated from a possible ALP signal in a previous Galactic pulsar analysis we fix the ALP parameters at values of gaγγ = 2.3 × 10−10 GeV−1 with ma = 3.6 neV, while searching for best-fitting parameters of (B, s) on a grid map over large ranges of B (1 nG-100 μG) and s (100 pc-10 Mpc). We find for 18 of 20 sources a favorable fit, particularly for Markarian 421 and NGC 1275, under H1 hypothesis in a likelihood fitting analysis. A significance level of 5.3 σ is obtained from a combinatorial analysis of all 20 sources in hypotheses testing using bootstrapping method, where likelihood results of local maxima for 20 sources are chosen in the B-s planes. The local maximum of each source refers to the minimization of energy requirement that are used to sustain a magnetic field with energy density ∝ B2s3. In a similar way, a 6 σ is also achieved when choosing all best-fitting parameters of (B, s) from the global maxima, instead of the local maxima. It is important to note that the best-fitting parameters of (B, s) would still remain valid for different combinations of ALP parameters, due to the internal degeneracy of the parameters (gaγγ , ma, B and s) that are used to describe spectral modulations on SEDs. The best-fitting values of (B, s) for 20 sources fall into ranges that are expected for large-scale magnetic fields present in relevant astrophysical environments in vicinity of the sources. A χ2 analysis is also performed for the purposes of consistency check and possibly extending the analyzed LAT energy (100 MeV-500 GeV) to TeV range. We obtain consistent results from χ2 analysis in terms of best-fitting parameter values and hypotheses testing.