Here are a few thoughts:

• "the Pearson correlation coefficient of mRNA abundance". I don't think this would be a good idea because mRNA abondance (that is, I guess, read count per gene) follows a long tail distribution that would skew the calculation the correlation. Instead, I would either use Spearman correlation on the read counts per gene, or Pearson correlation on (r)log-transformed read counts per genes.

• "Can I use deeptools?" I think that it is possible to use deeptools for that purpose with the bamCorrelate and plotCorrelation functions. Another approach, perhaps more common, is to calculate the number of read per genes, feed that into DESeq/EdgeR, normalize and transform the data (with rlog transformation or instance), then calculate correlation.

• "How to analyze the Pearson correlation coefficient of modification sites of nucleotide in mRNA?" Well, in that case I would not use a correlation because you do not have a continuous variable (it is just presence/absence from what I understand). Instead I would make a Venn diagram + FIsher's exact test to see if there is significant overlap between the modified sites in sample A vs B.

How to analyze the Pearson correlation coefficient of mRNA abundance between two biological replicates? Can I use deeptools?

Yes