Getting Started With Biological Data In R
R is the default tool for most people working with RNA-seq, microarray data, or ecological count data. Not because it is the prettiest tool, but because the infrastructure around it is thick enough to walk on without falling through. The Bioconductor project has been maintaining specialized packages for over two decades, and that means if someone else encountered a problem you are about to face, a package already exists for it. I still find myself installing DESeq2, edgeR, limma, and ggplot2 in roughly that order for most projects, and that order is not arbitrary. When you open a raw count matrix from an RNA-seq experiment, it looks like a spreadsheet with genes in rows and samples in columns. It is not. Those raw counts carry library size differences, composition bias, and mean-variance relationships that ordinary statistical tests cannot handle. The first thing most people get wrong is running a t-test on raw counts. The second thing is then wondering why their volcano plots look like garbage. You need to normalize first. With DESeq2 that means calling the median-of-ratios method, which estimates size factors from the geometric mean across samples. It sounds straightforward until you hit a sample with extreme composition bias — like one sample where half the reads map to a few highly expressed genes. The size factors will be distorted, and your downstream results will be too. I ran into this with a dataset where one condition had a contamination issue that only showed up after PCA, not during QC. The workaround was to filter out the most variable genes before estimating size factors, then re-run estimateSizeFactors() on the filtered object. It fixed the clustering. Visualization is where most beginners lose time. ggplot2 is the standard, but the learning curve is steep because the grammar of graphics is not intuitive until it is. You do not build plots by adding elements in sequence. You build them by mapping variables to aesthetic attributes and then layering geoms. A volcano plot is not a scatter plot with labels. It is a layer of points mapped to log2 fold change and negative log10 p-value, with a second layer of text labels applied conditionally to significant genes, over a background grid. The code takes about twelve lines the first time you write it. After that, it takes about four because you have a template.
Heatmaps are another trap. The default output from pheatmap() or even ComplexHeatmap can look impressive, but without proper row scaling and a meaningful clustering method, you are just displaying noise with colors. I once spent three days debugging a heatmap that showed what looked like clear condition-specific patterns. The patterns disappeared when I realized the clustering was driven by a single outlier sample that had a much higher sequencing depth than the rest. Re-running with scale="row" and removing that sample from the dendrogram calculation revealed the actual signal. The lesson was not about heatmaps. It was about checking your data before you decorate it.
The Practical Workflow
A typical pipeline starts with quality control using fastqc or MultiQC, then alignment with STAR or HISAT2, then quantification with featureCounts or Salmon. The output is a count matrix that you bring into R. From there, the workflow splits depending on your question. Differential expression uses DESeq2 or edgeR. Trajectory inference uses Monocle3 or slingshot. Single-cell analysis uses Seurat or scater. Each of these has its own preprocessing expectations, and none of them accept raw data without some preparation. With DESeq2, the pipeline is:
Get the Full Details
- Create a
colDatadataframe with your experimental design - Run
DESeqDataSetFromMatrix() - Filter low-count genes before dispersion estimation
- Run
DESeq() - Extract results with
results()and applylfcShrink()for more accurate log2 fold changes
The last step is the one most tutorials skip. The regularized log transformation from rlog() or the variance stabilizing transformation from vst() is essential for visualization. Raw log2(counts + 1) compresses high counts and inflates low counts in a way that distorts clustering. vst() corrects for that. I use it for PCA plots, heatmaps, and any distance-based analysis. For visualization, ggplot2 combined with ggrepel for label placement and patchwork for arranging multiple plots is the standard combination. The pals package provides good color schemes, but diverging palettes like RdBu or viridis work better for heatmaps than sequential ones. Colorblind-friendly palettes matter more than people admit. I stopped using rainbow() after a reviewer pointed out that my heatmap was unreadable to someone with deuteranopia. The fix was replacing it with viridis from viridisLite.
Where Things Break
R is not a universal solution. Memory usage becomes a real problem with single-cell datasets. A typical 10x Genomics dataset with 50,000 cells and 20,000 genes will consume several gigabytes in memory when loaded as a dense matrix. Seurat handles this with sparse matrices, but even then, operations like UMAP embedding can push a 16GB machine to its limits. If you are working with large datasets, you need either more RAM or a subset of cells for initial exploration. There is no way around this constraint. Batch effects are another area where people underestimate the damage. ComBat from the sva package works well for microarray data, but applying it to RNA-seq count data is wrong. You should correct batch effects on normalized, transformed data, not on raw counts. And even then, batch correction can remove biological signal if the batch is confounded with your condition. I worked on a project where the treatment group was processed on one day and the control on another. Batch correction would have erased the effect I was trying to detect. The only valid approach was to include batch as a covariate in the DESeq2 design formula: ~ batch + condition. This is a limitation of the data collection, not the analysis, but it is easy to miss if you just follow a tutorial blindly. Reproducibility is easier now than it was five years ago, but it is not automatic. Saving your session state with save.image() captures package versions, which is useful. But it does not capture the exact environment if you switched between conda and system packages mid-project. Using renv or conda to lock dependencies at the project level is the only reliable method. I lost two weeks of work to a package update that changed the default behavior of lfcShrink(). The new version shrank toward zero differently, and my effect sizes shifted enough to change the significance calls on borderline genes. Pinning versions with renv prevents this.
What To Learn First
Don't start with single-cell. Start with bulk RNA-seq differential expression. The concepts transfer, and the datasets are smaller, so you can debug faster. Learn how to read a DESeq2 results table, how to interpret the dispersion plot, and how to filter on adjusted p-value and log2 fold change simultaneously. These are the foundations. Everything else builds on them. For visualization, learn ggplot2 separately from your biological analysis. They are different skills. Knowing how to transform data for a plot is useful, but understanding layers, scales, and coordinates takes practice that is unrelated to statistics. I kept trying to combine both at the same time and made slow progress on both. Separating them cut my learning time roughly in half. The documentation for Bioconductor packages is generally better than CRAN packages, but it assumes you already know what you are looking for. The ?DESeq2 help page is thorough but dense. The vignettes are more useful for learning. Run them in order. Do not skip the one on independent filtering, which explains why DESeq2 filters low-count genes automatically and how that improves power. Skipping it means you will not understand why your manual filtering produced different results than the automatic pipeline.
If you need a starting point, the Bioconductor workflow site has a collection of end-to-end pipelines that cover everything from FASTQ to publication-ready figures. They are maintained by the package authors, so they tend to reflect current best practices rather than outdated conventions. I still use them as a reference even though I have been doing this for years. The documentation for OrganismDbi and AnnotationDbi packages is also worth reading if you plan to do any gene annotation work. Most people ignore it until they need it, and by then they are stuck trying to map gene IDs across species-specific databases. The field moves fast. New packages appear every release cycle. The core tools — DESeq2, edgeR, limma, ggplot2, Seurat — have been stable enough that investing time in them pays off for a long time. Don't chase every new tool. Learn the ones that solve the problems you actually encounter.
