warning messages in negative binomial regression in R (forward stepwise and backward stepwise)

Question

I am doing some count data analysis. The data is in this link: [1]: https://www.dropbox.com/s/q7fwqicw3ebvwlg/stackquestion.csv?dl=0

Column A is the count data, and other columns are the independent variables. At first I used Poisson regression to analyze it:

m0<-glm(A~.,data=d,family="poisson")
summary(m0)
#We see that the residual deviance is greater than the degrees of freedom so that we have over-dispersion. 

Call:
glm(formula = A ~ ., family = "poisson", data = d)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-28.8979   -4.5110    0.0384    5.4327   20.3809  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  8.7054842  0.9100882   9.566  < 2e-16 ***
B           -0.1173783  0.0172330  -6.811 9.68e-12 ***
C            0.0864118  0.0182549   4.734 2.21e-06 ***
D            0.1169891  0.0301960   3.874 0.000107 ***
E            0.0738377  0.0098131   7.524 5.30e-14 ***
F            0.3814588  0.0093793  40.670  < 2e-16 ***
G           -0.3712263  0.0274347 -13.531  < 2e-16 ***
H           -0.0694672  0.0022137 -31.380  < 2e-16 ***
I           -0.0634488  0.0034316 -18.490  < 2e-16 ***
J           -0.0098852  0.0064538  -1.532 0.125602    
K           -0.1105270  0.0128016  -8.634  < 2e-16 ***
L           -0.3304606  0.0155454 -21.258  < 2e-16 ***
M            0.2274175  0.0259872   8.751  < 2e-16 ***
N            0.2922063  0.0174406  16.754  < 2e-16 ***
O            0.1179708  0.0119332   9.886  < 2e-16 ***
P            0.0618776  0.0260646   2.374 0.017596 *  
Q           -0.0303909  0.0060060  -5.060 4.19e-07 ***
R           -0.0018939  0.0037642  -0.503 0.614864    
S            0.0383040  0.0065841   5.818 5.97e-09 ***
T            0.0318111  0.0116611   2.728 0.006373 ** 
U            0.2421129  0.0145502  16.640  < 2e-16 ***
V            0.1782144  0.0090858  19.615  < 2e-16 ***
W           -0.5105135  0.0258136 -19.777  < 2e-16 ***
X           -0.0583590  0.0043641 -13.373  < 2e-16 ***
Y           -0.1554609  0.0042604 -36.489  < 2e-16 ***
Z            0.0064478  0.0001184  54.459  < 2e-16 ***
AA           0.3880479  0.0164929  23.528  < 2e-16 ***
AB           0.1511362  0.0050471  29.945  < 2e-16 ***
AC           0.0557880  0.0181129   3.080 0.002070 ** 
AD          -0.6569099  0.0368771 -17.813  < 2e-16 ***
AE          -0.0040679  0.0003960 -10.273  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 97109.0  on 56  degrees of freedom
Residual deviance:  5649.7  on 26  degrees of freedom
AIC: 6117.1

Number of Fisher Scoring iterations: 6

Then I think I should use negative binomial regression for the over-dispersion data. Since you can see I have many independent variables, and I wanted to select the important variables. And I decide to use stepwise regression to select the independent variable. At first, I create a full model:

full.model <- glm.nb(A~., data=d,maxit=1000)
# when not indicating maxit, or maxit=100, it shows Warning messages: 1: glm.fit: algorithm did not converge; 2: In glm.nb(A ~ ., data = d, maxit = 100) : alternation limit reached
# When indicating maxit=1000, the warning message disappear.

summary(full.model)

Call:
glm.nb(formula = A ~ ., data = d, maxit = 1000, init.theta = 2.730327193, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.5816  -0.8893  -0.3177   0.4882   1.9073  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)   
(Intercept) 11.8228596  8.3004322   1.424  0.15434   
B           -0.2592324  0.1732782  -1.496  0.13464   
C            0.2890696  0.1928685   1.499  0.13393   
D            0.3136262  0.3331182   0.941  0.34646   
E            0.3764257  0.1313142   2.867  0.00415 **
F            0.3257785  0.1448082   2.250  0.02447 * 
G           -0.7585881  0.2343529  -3.237  0.00121 **
H           -0.0714660  0.0343683  -2.079  0.03758 * 
I           -0.1050681  0.0357237  -2.941  0.00327 **
J            0.0810292  0.0566905   1.429  0.15291   
K            0.2582978  0.1574582   1.640  0.10092   
L           -0.2009784  0.1543773  -1.302  0.19296   
M           -0.2359658  0.3216941  -0.734  0.46325   
N           -0.0689036  0.1910518  -0.361  0.71836   
O            0.0514983  0.1383610   0.372  0.70974   
P            0.1843138  0.3253483   0.567  0.57105   
Q            0.0198326  0.0509651   0.389  0.69717   
R            0.0892239  0.0459729   1.941  0.05228 . 
S           -0.0430981  0.0856391  -0.503  0.61479   
T            0.2205653  0.1408009   1.567  0.11723   
U            0.2450243  0.1838056   1.333  0.18251   
V            0.1253683  0.0888411   1.411  0.15820   
W           -0.4636739  0.2348172  -1.975  0.04831 * 
X           -0.0623290  0.0508299  -1.226  0.22011   
Y           -0.0939878  0.0606831  -1.549  0.12142   
Z            0.0019530  0.0015143   1.290  0.19716   
AA          -0.2888123  0.2449085  -1.179  0.23829   
AB           0.1185890  0.0696343   1.703  0.08856 . 
AC          -0.3401963  0.2047698  -1.661  0.09664 . 
AD          -1.3409002  0.4858741  -2.760  0.00578 **
AE          -0.0006299  0.0051338  -0.123  0.90234   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(2.7303) family taken to be 1)

    Null deviance: 516.494  on 56  degrees of freedom
Residual deviance:  61.426  on 26  degrees of freedom
AIC: 790.8

Number of Fisher Scoring iterations: 1


              Theta:  2.730 
          Std. Err.:  0.537 

 2 x log-likelihood:  -726.803 

When not indicating maxit, or maxit=100, it shows Warning messages: 1: glm.fit: algorithm did not converge; 2: In glm.nb(A ~ ., data = d, maxit = 100) : alternation limit reached. When indicating maxit=1000, the warning message disappear.

Then I create a first model:

first.model <- glm.nb(A ~ 1, data = d)

Then I tried the forward stepwise regression:

step.model <- step(first.model, direction="forward", scope=formula(full.model))
#Error in glm.fit(X, y, wt, offset = offset, family = object$family, control = object$control) : 
#NA/NaN/Inf in 'x'
#In addition: Warning message:
#  step size truncated due to divergence 
#What is the problem?

It gives me error message: Error in glm.fit(X, y, wt, offset = offset, family = object$family, control = object$control) : NA/NaN/Inf in 'x' In addition: Warning message: step size truncated due to divergence

I also tried the backward regression:

step.model2 <- step(full.model,direction="backward")

#the final step
Step:  AIC=770.45
A ~ B + C + E + F + G + H + I + K + L + R + T + V + W + Y + AA + 
    AB + AD

       Df Deviance    AIC
<none>      62.375 770.45
- AB    1   64.859 770.93
- H     1   65.227 771.30
- V     1   65.240 771.31
- L     1   65.291 771.36
- Y     1   65.831 771.90
- B     1   66.051 772.12
- C     1   67.941 774.01
- AA    1   69.877 775.95
- K     1   70.411 776.48
- W     1   71.526 777.60
- I     1   71.863 777.94
- E     1   72.338 778.41
- G     1   73.344 779.42
- F     1   73.510 779.58
- AD    1   79.620 785.69
- R     1   80.358 786.43
- T     1   95.725 801.80
Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: algorithm did not converge 
3: glm.fit: algorithm did not converge 
4: glm.fit: algorithm did not converge 

My question is: Why it is different in using forward and backward stepwise regression? And why do I get the error message when performing forward selection? Also, what exactly do these warning messages mean? And how should I deal with it?

I am not a stats person but need to conduct statical analysis for my research data. So I am struggling in learning how to do different regression analyses using real data. I searched online for similar questions but I still could understand ... And please let me know if I did anything wrong in my regression analysis. I would really appreciate it if you could help me with these questions!