Data Analysis For Life Sciences

Getting Your First Single-Cell Dataset Through a Viable Pipeline

I spent three days debugging a sequencing run that turned out to be completely fine. The problem was my own quality control filter set too aggressively on the ambient RNA estimation. I had removed cells that were actually viable because my tool flagged them as outliers. This kind of thing happens constantly when you first start doing Data Analysis For Life Sciences, especially if you are working with single-cell RNA sequencing data. You follow a published pipeline exactly, get poor results, and spend a week questioning your assumptions before realizing the preprocessing parameters were tuned for a different tissue type. The first real step is not running the pipeline. It is understanding what your sequencing platform generated. Raw FASTQ files from Illumina machines contain adapter sequences, low-quality bases at the read ends, and sometimes index hopping artifacts if you multiplexed samples heavily. Trimming with tools like Trim Galore or fastp before alignment usually recovers a meaningful fraction of otherwise discarded reads. Skipping this step costs you roughly 5 to 12 percent of usable reads on a standard 10x Genomics run, which translates directly into lower gene detection per cell.

Data Analysis For Life Sciences

The field covers a wide range of analytical work, but most projects you encounter will fall into one of several categories. Bulk RNA-seq differential expression, single-cell transcriptomics, variant calling from whole exome or genome sequencing, proteomics mass spectrometry quantification, and microbiome 16S or shotgun metagenomics. Each category has its own standard tools and its own common failure modes. Beginners often try to apply the same workflow to all of them. That does not work well. I will walk through a bulk RNA-seq differential expression project because it is the most common entry point and it teaches the core principles you need for everything else. The general workflow involves raw read QC, alignment or pseudoalignment, quantification, normalization, differential expression testing, and then functional interpretation. The order matters less than you might think, but the dependency between steps is strict. You cannot normalize accurately until you have reliable counts, and you cannot get reliable counts until your reads are not full of adapter contamination. Here is the part that is not in most tutorials. Most people learn to use DESeq2 or edgeR and then stop. They produce a volcano plot and call it done. The thing nobody tells you is that batch effects will destroy your results if you ignore them during the experimental design phase, not during the analysis phase. If you collected all your treated samples on Tuesday and all your controls on Wednesday, the day variable is a confounder that no amount of batch correction afterward will cleanly separate from the biological signal. The fix is randomization during sample processing, but if you are stuck with a messy dataset, you include the batch variable in your design matrix. In DESeq2 that looks like adding a column to your colData and specifying ~batch + condition in your model formula. It is not a perfect solution, but it is usually better than ignoring the problem.

Another practical detail that saves serious time is using tximport instead of importing raw counts directly into DESeq2 when your quantification came from Salmon or kallisto. These tools produce transcript-level estimates, and tximport collapses them to gene-level counts while also carrying the effective length information into the DESeq2 model. Using raw gene counts from Salmon without tximport means you lose that length correction, and your dispersion estimates become slightly biased. The difference is small on well-annotated genomes but noticeable when you are working with a non-model organism or a poorly annotated transcriptome. I remember a project where we were analyzing mouse brain tissue with ~15 million reads per sample. A colleague insisted on using a 0.1 CPM filtering threshold before differential expression, which removed about 40 percent of the genes. The final result had fewer significant hits than a run with a much looser filter, and the effect sizes were distorted for lowly expressed transcription factors that happened to fall below that arbitrary cutoff. We ended up using the default DESeq2 independent filtering, which is built into the function and automatically optimizes the threshold based on the relationship between mean expression and raw p-value distribution. It usually keeps more genes without inflating the false discovery rate. For normalization, the median-of-ratios method in DESeq2 and the TMM method in edgeR are both solid choices for standard RNA-seq experiments. They perform similarly on well-controlled datasets. The situation changes when you have extreme composition bias, meaning a small number of genes dominate the read counts in one condition. This happens frequently in experiments involving strong transcriptional activation or knockdown studies. TMM handles this better than median-of-ratios in those cases. If you have a dataset where the top 10 genes account for more than 30 percent of total reads in one group, switch to TMM or consider using a specialized normalization method like upper-quartile normalization.

When you move into multiple testing correction, the Benjamini-Hochberg procedure is standard, but it assumes a reasonable proportion of true null hypotheses. In practice, with RNA-seq you often have thousands of truly differentially expressed genes, which means BH can be overly conservative. Some researchers use the storey q-value approach via the qvalue R package as an alternative. It estimates pi0, the proportion of true nulls, and adjusts accordingly. The difference is usually modest, but it can matter when you are working with a small sample size and every significant hit is precious. Functional interpretation is where many projects stall. You have your list of differentially expressed genes and now you need to make sense of them. GO enrichment with clusterProfiler is the most common route. It is straightforward but has limitations. Gene ontology terms are hierarchical, and standard enrichment tests do not account for term-to-term dependencies. This means you get a long list of overlapping terms that say basically the same thing in slightly different wording. The g:Profiler web tool or the enrichR package handle this somewhat better by reordering and deduplicating results. For pathway analysis, MSigDB collections are more comprehensive than traditional KEGG alone. I usually pull from the Hallmark gene sets and the C2 curated pathways collection. The Hallmark sets are pre-filtered to remove redundancy, which makes downstream interpretation cleaner. Volume and speed become real issues when you scale up. A standard bulk RNA-seq differential expression analysis on a modern laptop with 16 GB of RAM takes roughly 10 to 20 minutes from raw FASTQ to volcano plot if you use a pseudoalignment approach with Salmon followed by tximport and DESeq2. A full alignment with STAR followed by featureCounts takes longer, usually 45 to 90 minutes depending on sample count and genome complexity. The speed difference is meaningful when you are processing 50 or 100 samples. The accuracy difference is negligible for standard gene-level quantification.

Here is a scenario that catches people off guard. You run a batch of samples through a pipeline, everything looks fine, and then you discover that two of your samples have a dramatically different GC content profile compared to the rest. This often happens when there is a library preparation issue or when the RNA integrity number dropped below 7 for those particular samples. GC bias skews quantification in a way that standard normalization does not fully correct. The workaround is to check GC content distributions early, ideally right after the alignment or pseudoalignment step, using tools like samtools idxstats or Salmon's own bias correction output. If you detect a strong outlier, you can either exclude the sample or use a GC-bias correction method like cqn in R before downstream analysis. Including a biased sample in a differential expression test usually produces false positives in genes that correlate with the GC content artifact rather than with your actual experimental condition. For variant calling, the pipeline is more rigid. You align with BWA-MEM or minimap2, mark duplicates with Picard or sambamba, realign around indels if you are following GATK best practices, and then call variants with HaplotypeCaller or FreeBayes. The GATK pipeline is the industry standard but it is also slow and memory hungry. A whole exome sequencing run processed through the full GATK pipeline on a single node can take 3 to 6 hours. freebayes is faster but less rigorous in its Bayesian model for complex regions. If you need speed and your samples are relatively clean, freebayes or varscan2 can work. For clinical or publication-grade variant calls, GATK is still the safest choice. A common mistake in variant calling is not hard-filtering the raw VCF output. Raw calls from HaplotypeCaller contain a lot of noise, especially in repetitive regions and near indels. Applying the Variant Quality Score Recalibration or at minimum a hard filter based on QD, FS, MQ, and MQRankSum values removes a significant portion of false positives. Skipping this step on a WES dataset can leave you with thousands of spurious variants that look plausible until you cross-reference them with dbSNP and gnomAD and realize most of them do not exist in any population database.

Multi-omics integration is another area where expectations often exceed reality. People want to combine transcriptomics, proteomics, and methylation data into a single unified model. The most practical approach right now is still sequential integration. You analyze each data type separately, extract features of interest from each, and then use overlap analysis or pathway-level correlation to find convergence points. Joint dimensionality reduction methods like MOFA or Seurat's CCA-based integration exist and can work, but they require careful parameter tuning and they often introduce artifacts that are difficult to distinguish from biological signal. I have seen projects where the integration itself created clusters that matched the batch structure rather than the biological condition. Cross-validation on held-out samples is essential if you go this route. Proteomics data from mass spectrometry introduces its own set of challenges. Missing values are a major issue. Unlike RNA-seq where a zero count usually means the gene was not expressed or not captured, a missing value in proteomics can mean the peptide was not detected due to ionization efficiency, dynamic range limits, or stochastic sampling in data-dependent acquisition mode. Imputing missing values is risky. Random imputation from a narrow distribution works for some downstream analyses but it can create artificial patterns. A more defensible approach is to use left-censored imputation, where you replace missing values with values drawn from a normal distribution centered slightly below the observed minimum. The minprob function in the-impute R package implements something close to this. It is not perfect, but it is better than filling missing values with zeros or with the column mean. Microbiome analysis follows a different logic entirely. You are usually working with relative abundance data, which means all the compositional constraints apply. A change in one taxon affects the apparent abundance of every other taxon simply because the total is normalized to 100 percent. This makes standard correlation analysis misleading. Tools like SparCC or proportionally scaled log-ratio transformations address this to some extent. alpha diversity metrics like Shannon and Simpson are straightforward to calculate with QIIME2 or the phyloseq R package. Beta diversity with Bray-Curtis or UniFrac distances is where most people spend their time, and PERMANOVA with the adonis function in vegan is the standard test for group differences. The caveat is that PERMANOVA is sensitive to differences in dispersion between groups, not just differences in centroids. If your groups have different within-group variability, a significant PERMANOVA result might reflect dispersion rather than a location shift. Checking dispersion with betadisper and running the associated permutation test is a quick way to verify what your PERMANOVA is actually detecting.

Storage and data management are not glamorous but they determine whether your project survives past the first month. Raw sequencing data for a moderate RNA-seq study with 60 samples at 30 million reads each generates roughly 3 to 4 terabytes of FASTQ data. Processed files add another 500 gigabytes. If you are archiving for reproducibility, you need a strategy for versioning your code and your environment. Docker or Singularity containers solve the dependency problem but they add complexity to the workflow. R environments managed with packrat or renv are simpler for R-based pipelines. Python projects benefit from conda or pixi environments. The key is documenting exactly what you used. A pipeline that works today and cannot be reproduced in six months because an R package updated and broke a function is worse than no pipeline at all.

Tools and Resources Worth Knowing About

DESeq2 and edgeR are the workhorses for bulk RNA-seq. Seurat and scanpy dominate single-cell work. Salmon and kallisto are fast and accurate for transcript quantification. STAR and HISAT2 are the alignment standards. GATK is the variant calling reference. QIIME2 and mothur cover microbiome workflows. For visualization, ggplot2 in R and matplotlib and seaborn in Python remain the most reliable options. There are newer tools that promise more, but these established packages have the documentation, community support, and track record that matter when your analysis is under review. Data repositories are another practical concern. GEO and ArrayExpress hold most public RNA-seq datasets. The SRA is the primary source for raw reads, but downloading from SRA directly can be slow. Using the sra-toolkit with fasterq-dump is generally faster than downloading individual FASTQ files through the web interface. For clinical or protected datasets, dbGaP and the EU’sEGA require controlled access applications. Processing times vary from a few days to several weeks depending on the review board. Cloud computing has changed the economics of large-scale analysis. Terra from Broad Institute, AWS Partnered Solutions for Genomics, and Google Cloud Life Sciences all offer pre-configured environments for common pipelines. Costs can add up quickly. A single WGS analysis on cloud infrastructure with optimal settings runs roughly $2 to $5 per sample in compute costs, not including storage. For a one-off project it is manageable. For ongoing work, the costs accumulate. The main advantage is that you do not need to maintain local infrastructure and you can scale up instantly when a large batch arrives.

The hardest part of this work is usually not the technical execution. It is the decision-making at each step. Which filter threshold? Which normalization method? Which batch correction approach? Which version of the genome annotation? There is rarely a single correct answer. The best approach is to document every choice, run sensitivity analyses on the key decisions, and report the uncertainty in your final results. Readers and reviewers increasingly expect this level of transparency. It also protects you when you come back to the project six months later and cannot remember why you made a particular parameter choice.