Just remove the undesirable equations and variables:
Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 +
(T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;
NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];
inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};
eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];
var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];
(# = 400) & /@ inner;
solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;
ListPlot@s
BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.