We have applied density functional calculations to investigate the formation of boron nitride islands in the triangular and hexagonal graphene nanoflakes (GNFs). Based on the calculated reaction energies, BN doping of GNFs is found to be endothermic. The variation of reaction energies indicates that BN doping of the GNFs follow previously proposed 'hexagon filling" rule. Based on the calculated NICS indexes, aromaticity inside these molecules can be explained within the framework of the Clar's Sextet Theory. In the case of BN-doped hexagonal and triangular GNFs, NICS(1) index shows aromaticity and nonaromaticity for the borazine like substructures, respectively. Moreover, the substitution of one BN unit in the 1,2 azaborine like rings leads to the increase of aromaticity and nonaromaticity in hexagonal and triangular GNFs, respectively. Based on our results, the BN doped hexagonal GNFs not only lose electrons somewhat more easily to form positive ions, but also obtain electrons more easily to form anions.