I could not find any good source for this type of calculation. I'd be very appreciated if someone who is experienced with that to show it step by step
I assume that your structure is fully relaxed.
(I) Self-consistent field calculation at DFT-PBE level with spin-orbit coupling (SOC) to obtain the converged charge density.
#SOC related tags: LSORBIT=.TRUE. SAXIS= 0 0 1 ISYM=0
(II) Usual band calculation at DFT-PBE level with SOC by reading the converged charge density at step I. Note that the k-sampling assigned in KPOINTS can be found in OUTCAR and used in a subsequent step.)
grep "k-points in reciprocal lattice and weights" OUTCAR -A n (n is the total k points in your band path.)
(III) Self-consistent field calculation at HSE06 level with SOC by reading the converged charge density at step I.
#HSE06 related tags: LHFCALC = .TRUE. HFSCREEN = 0.2 ALGO = Damped TIME = 0.4 PRECFOCK=F
(IV) Non-self-consistent field calculation at HSE06 level with SOC by reading the converged wave function at step III. ( One special point is we first need to move the IBZKPT at step III to the KPOINTS file and also add the k-points read from OUTCAR at step II to this KPOINTS file with zero weight.)
#HSE06 related tags: LHFCALC = .TRUE. HFSCREEN = 0.2 ALGO = Normal TIME = 0.4 PRECFOCK=F
(V) Read the EIGENVALUE to plot bandstructure (You can make a comparison with the result of step II). Note that the EIGENVALUE from SCF IBZKPT at step III needs to be excluded. For this step, you may take a look at this post: How to remove the fake k-points from vasprun.xml file in the calculation of HSE06 band structure?