cell type specific genes - heatmap using rank based approach
The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:

&gt; expr &lt;- matrix(runif(9*24000), nrow=24000, ncol=9)
&gt; top_genes_each_celltype &lt;- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
&gt; dim(top_genes_each_celltype)
[1] 40 9
&gt; top_genes_all_celltypes &lt;- unique(as.vector(top_genes_each_celltype))
&gt; length(top_genes_all_celltypes)
[1] 358


Then you can build a biclustered heatmap with expr[top_genes_all_celltypes,], for example, using heatmap2.