WGCNA correcting p value for varying sample numbers in each of the trait columns
You can use the function corAndPvalue. Simply replace the lines

moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)


by something like

cp = corAndPvalue(MEs, status)
moduleTraitCor = cp$cor
moduleTraitPvalue = cp$p


corAndPvalue takes automatically into account the number of valid samples for each element of the resulting correlation matrix.