Get["FeynCalc`"]

fAmp=-(EL^2*Spinor[Momentum[p],0,1].(-((EL*CKM[1,1]*DiracGamma[7])/(Sqrt[2]*SW))).(MD-DiracGamma[Momentum[q1]]).((I*EL*Conjugate[CKM[1,1]]*DiracGamma[LorentzIndex[Lor1]].DiracGamma[7])/(Sqrt[2]*SW)).Spinor[Momentum[p],MU,1]*FeynAmpDenominator[PropagatorDenominator[-p-q1,MW],PropagatorDenominator[-p-q1,MW*Sqrt[GaugeXi[W]]],PropagatorDenominator[-p-q1,MW*Sqrt[GaugeXi[W]]],PropagatorDenominator[-q1,MD],PropagatorDenominator[q2,MH]]*Pair[LorentzIndex[Lor1],LorentzIndex[Lor2]]*Pair[LorentzIndex[Lor2],LorentzIndex[Lor4]]*Pair[LorentzIndex[Lor4],Momentum[-p-q1]]*SumOver[SUNFIndex[Col1],3,External]*SumOver[SUNFIndex[Col2],3,External])/(2048*Pi^8*SW^2)-(EL^2*Spinor[Momentum[p],0,1].(-((EL*CKM[1,1]*DiracGamma[7])/(Sqrt[2]*SW))).(MD-DiracGamma[Momentum[q1]]).((I*EL*Conjugate[CKM[1,1]]*DiracGamma[LorentzIndex[Lor1]].DiracGamma[7])/(Sqrt[2]*SW)).Spinor[Momentum[p],MU,1]*FeynAmpDenominator[PropagatorDenominator[-p-q1,MW],PropagatorDenominator[-p-q1,MW*Sqrt[GaugeXi[W]]],PropagatorDenominator[-p-q1,MW*Sqrt[GaugeXi[W]]],PropagatorDenominator[-p-q1,MW*Sqrt[GaugeXi[W]]],PropagatorDenominator[-q1,MD],PropagatorDenominator[q2,MH]]*(1-GaugeXi[W])*Pair[LorentzIndex[Lor1],Momentum[-p-q1]]*Pair[LorentzIndex[Lor2],LorentzIndex[Lor4]]*Pair[LorentzIndex[Lor2],Momentum[p+q1]]*Pair[LorentzIndex[Lor4],Momentum[-p-q1]]*SumOver[SUNFIndex[Col1],3,External]*SumOver[SUNFIndex[Col2],3,External])/(2048*Pi^8*SW^2)

ToTFI[Contract[fAmp],q1,q2,p]