This is all done in Jupyter (formerly IPython) and IRkernel.
Visit http://159.203.72.130:8888 to use this same notebook in your browser.
(...unless you're reading this later, of course. Go fire up your own docker container with "docker run -d -p 8888:8888 jupyter/r-notebook" or something.)

(thanks Creighton)
Do you need plotting or visualization? Use ggplot2. Completely ignore built-in plotting.
Do you need something else? Search CRAN.
Does no CRAN package solve your problems? Do you need to write "real"(tm) software for production? Strongly consider giving up.
Datasetfilter, slice - Select rows (filter is by predicate, slice is by position)arrange - Reorder rowsselect, rename - Select columnsdistinct - Choose only distinct rowsmutate, transmute - Make new columns from existing onessummarise - 'Peel off' one grouping level (or collapse frame to one row) with aggregate functionssample_n, sample_frac - Randomly sample (by count or by percentage)group_by - Group observations (most of above worked on grouped observations)x %>% f(y) = f(x,y)data.table via dtplyr)

(though actually M >> N is useful too and is the basis for kernel methods such as SVMs)
Then for $X$ containing samples as row vectors, $Y=X\mathcal{P}_x$: $$ \mathcal{P}_x= \begin{bmatrix} 2 & 0 \\ -1 & 0 \\ 0 & 1 \\ 0 & 1 \end{bmatrix} $$
Consider data expressed as an $n\times m$ matrix with each column representing one feature (of $m$) and each row one sample (of $n$):
$$ X= \begin{bmatrix} a_1 & b_1 & c_1 & \cdots\\ a_2 & b_2 & c_2 & \cdots\\ a_3 & b_3 & c_3 & \cdots\\ \vdots & \vdots & \vdots & \ddots\\ a_n & b_n & c_n & \cdots\\ \end{bmatrix} $$(if you want to know why it is $\frac{1}{n-1}$ and not $\frac{1}{n}$, ask a statistics PhD or something)
Treating $A$ and $B$ as vectors:
$$\sigma_{AB}^2=\frac{A\cdot B}{n-1}$$Recalling our data matrix:
$$ X= \begin{bmatrix} a_1 & b_1 & c_1 & \cdots\\ a_2 & b_2 & c_2 & \cdots\\ a_3 & b_3 & c_3 & \cdots\\ \vdots & \vdots & \vdots & \ddots\\ a_n & b_n & c_n & \cdots\\ \end{bmatrix} $$It can be rewritten as column vectors:
$$ X= \begin{bmatrix} a_1 & b_1 & c_1 & \cdots\\ a_2 & b_2 & c_2 & \cdots\\ a_3 & b_3 & c_3 & \cdots\\ \vdots & \vdots & \vdots & \ddots\\ a_n & b_n & c_n & \cdots\\ \end{bmatrix} =\begin{bmatrix} A & B & C & \cdots \end{bmatrix} $$Then covariance matrix is:
$$\mathbf{S}_X=\frac{X^\top X}{n-1}= \begin{bmatrix} \sigma_{A}^2 & \sigma_{AB}^2 & \sigma_{AC}^2 & \sigma_{AD}^2 & \cdots \\ \sigma_{AB}^2 & \sigma_{B}^2 & \sigma_{BC}^2 & \sigma_{BD}^2 & \cdots \\ \sigma_{AC}^2 & \sigma_{BC}^2 & \sigma_{C}^2 & \sigma_{CD}^2 & \cdots \\ \sigma_{AD}^2 & \sigma_{BD}^2 & \sigma_{CD}^2 & \sigma_{D}^2 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} $$