How this document works: every step includes (1) a plain-English explanation, (2) the underlying algorithm/statistics in depth, (3) the exact file format with an annotated example, (4) an example command as an R Markdown code chunk (eval=FALSE — illustrative, not meant to knit standalone), and (5) a citation. Commands reflect standard documented usage — always check a tool’s current --help/docs before production use, since flags and versions change.

How to use this document: every step includes (1) a plain-English explanation, (2) the underlying algorithm/statistics in more depth, (3) the exact file format with an annotated example, (4) the real command a bioinformatician would run, and (5) a citation. Command syntax is illustrative of standard usage as documented in each tool’s official manual — always check --help or the tool’s current docs before running in production, since flags/versions change.

0. Sequencing Instruments

Platform	Mechanism	Read length	Raw output	Use case
Illumina NovaSeq X / 6000	Bridge amplification + sequencing-by-synthesis (SBS): each cluster of identical DNA fragments on a flow cell emits a color-coded light pulse per base added	50–300 bp, paired-end	`.bcl` → converted to `.fastq.gz`	WGS, WES, RNA-seq, scRNA-seq
Oxford Nanopore	Single DNA/RNA strand threaded through a protein nanopore in a membrane; ionic current disruption pattern is decoded (“basecalled”) into bases in real time	Kb–Mb, no fixed limit	`.fast5`/`.pod5` (raw signal) → basecalled to `.fastq.gz`	Structural variants, de novo assembly, direct RNA sequencing
PacBio HiFi	Circular Consensus Sequencing (CCS): reads the same circularized fragment many times, then computes a consensus, cancelling out random errors	10–25 kb, ~99.9% accuracy	`.bam` (subreads) → CCS `.bam`/`.fastq.gz`	Long-read WGS, full-length isoform sequencing (Iso-Seq)
10x Genomics Chromium	Gel-bead-in-emulsion (GEM): each droplet contains one cell + one bead coated in millions of copies of a unique barcode oligo, which tags every RNA molecule from that cell before pooled Illumina sequencing	Short, barcoded reads	`.fastq.gz` (R1=barcode+UMI, R2=cDNA)	scRNA-seq, scATAC-seq, spatial

References: Illumina, “An Introduction to Next-Generation Sequencing Technology,” Illumina.com white paper; Wang Y., Zhao Y., Bollas A., Wang Y., Au K.F. (2021), “Nanopore sequencing technology, bioinformatics and applications,” Nature Biotechnology 39:1348–1365; Wenger A.M. et al. (2019), “Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome,” Nature Biotechnology 37:1155–1162 (defines PacBio HiFi); Zheng G.X.Y. et al. (2017), “Massively parallel digital transcriptional profiling of single cells,” Nature Communications 8:14049 (10x GEM method).

Key fact every bioinformatician should know: Illumina base calls are accompanied by a Phred quality score Q, where Q = -10 × log10(P_error). Q30 means a 1-in-1000 chance that base is wrong; Q20 means 1-in-100. This scoring system originates from the original Phred basecalling software (Ewing B. & Green P., 1998, Genome Research 8:186–194) and is universal across nearly all downstream tools.

1. FASTQ Format — anatomy in full detail

@A00123:45:HG7NKDSXX:1:1101:5678:1000 1:N:0:CGATGT
NAGCTGACGTTTGCAAGGCTAGCATGCATGCATGCATG
+
#FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:

Line 1 (@...): read identifier — instrument, run, flow cell, tile, x/y cluster coordinates, and (after the space) mate pair number (1 or 2), filter status, and the sample index barcode.
Line 2: the called sequence (N = base could not be determined).
Line 3 (+): separator, optionally repeats the ID (legacy).
Line 4: ASCII-encoded Phred quality, one character per base of line 2. The encoding offset is 33 (! = Q0), called “Phred+33,” the modern universal standard (older Illumina used Phred+64 — a common source of bugs when handling legacy data).

Paired-end sequencing produces two files, sample_R1.fastq.gz and sample_R2.fastq.gz, where read n in R1 and read n in R2 are the two ends of the same physical DNA fragment, sequenced from opposite directions (“forward” and “reverse”).

Reference: Cock P.J.A. et al. (2010), “The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants,” Nucleic Acids Research 38(6):1767–1771 — the formal specification.

1a. Quality Control

FastQC scans every read and reports, among other metrics: per-base quality distribution, per-sequence GC content (should roughly match expected genome/transcriptome GC%), overrepresented sequences (often adapters or PCR artifacts), and sequence duplication levels.

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc/

Reference: Andrews S. (2010). FastQC: A Quality Control tool for High Throughput Sequence Data. Babraham Bioinformatics. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

MultiQC parses the output logs of many tools (FastQC, aligners, dedup tools, etc.) across many samples into one interactive report — a standard final step of nearly every pipeline run.

multiqc qc/ -o multiqc_report/

Reference: Ewels P., Magnusson M., Lundin S., Käller M. (2016). “MultiQC: summarize analysis results for multiple tools and samples in a single report.” Bioinformatics 32(19):3047–3048.

1b. Trimming

fastp is a single, fast, combined QC+trimming tool that removes adapter sequences (auto-detected by overlap analysis between R1/R2), trims low-quality tails via a sliding window, and filters out reads below a minimum length or mean quality.

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o trim_R1.fastq.gz -O trim_R2.fastq.gz \
      --detect_adapter_for_pe -j fastp.json -h fastp.html

Reference: Chen S., Zhou Y., Chen Y., Gu J. (2018). “fastp: an ultra-fast all-in-one FASTQ preprocessor.” Bioinformatics 34(17):i884–i890.

What every bioinformatician should know: always inspect the duplication level and adapter content plots before trusting downstream results — high duplication can indicate low library complexity (too little input DNA/RNA) rather than a processing error, and won’t be “fixed” by trimming; it needs to be handled at the mark-duplicates step (Section 2) or flagged as a QC failure.

2. Alignment / Mapping (DNA)

2a. Reference preparation

bwa index GRCh38.fa            # builds the FM-index / suffix array (.amb .ann .bwt .pac .sa)
samtools faidx GRCh38.fa       # builds .fai — byte offsets for random access
gatk CreateSequenceDictionary -R GRCh38.fa   # builds .dict — required by GATK/Picard

The human reference genome (GRCh38/hg38, or the newer telomere-to-telomere CHM13) is distributed as .fasta: a header line >chr1 followed by the sequence, wrapped at ~60–80 characters per line. Reference genome release: Genome Reference Consortium, GRCh38.p14, ncbi.nlm.nih.gov/grc; Nurk S. et al. (2022), “The complete sequence of a human genome (T2T-CHM13),” Science 376:44–53.

2b. Alignment algorithm

BWA-MEM uses the Burrows-Wheeler Transform to build a compressed, searchable index of the genome, then finds “seed” exact matches for each read and extends them with Smith-Waterman local alignment to tolerate mismatches/small indels — the industry standard for short-read DNA.

bwa mem -t 8 -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA" \
  GRCh38.fa trim_R1.fastq.gz trim_R2.fastq.gz > aligned.sam

The -R read-group tag is mandatory for GATK downstream — it records sample name (SM), platform (PL), and library, which GATK uses to correctly handle multi-sample or multi-lane data.

For long reads, minimap2 uses minimizer-based seeding (a subsampling of k-mers) suited to the much higher per-base error rate of Nanopore/PacBio reads.

minimap2 -ax map-hifi GRCh38.fa reads.fastq.gz > aligned.sam

References: Li H. (2013). “Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM.” arXiv:1303.3997; Li H. (2018). “Minimap2: pairwise alignment for nucleotide sequences.” Bioinformatics 34(18):3094–3100.

2c. SAM/BAM format — full anatomy

SAM (Sequence Alignment/Map) is tab-delimited text; BAM is its binary, compressed, indexable equivalent. Every alignment line has 11 mandatory fields:

#	Field	Meaning	Example
1	QNAME	read name	`A00123:45:HG7NKDSXX:1:1101:5678:1000`
2	FLAG	bitwise code: paired, mapped, reverse strand, secondary, duplicate, etc.	`99`
3	RNAME	chromosome	`chr7`
4	POS	1-based leftmost mapping position	`140453136`
5	MAPQ	mapping quality, Phred-scaled confidence this position is correct	`60`
6	CIGAR	alignment operations: M=match/mismatch, I=insertion, D=deletion, S=soft-clip, N=skipped (splice)	`76M` or `40M2I34M`
7	RNEXT	chromosome of mate	`=`
8	PNEXT	position of mate	`140453310`
9	TLEN	inferred fragment/insert size	`250`
10	SEQ	read sequence	`AGCTGACG...`
11	QUAL	Phred quality string	`FFFFFFFF...`

The FLAG field is a bitmask — e.g., FLAG 99 = 1 (paired) + 2 (proper pair) + 32 (mate reverse strand) + 64 (first in pair). Every bioinformatician should be comfortable decoding these bits (e.g., via samtools flags or the Broad’s flag explainer) since filtering on FLAG (e.g., excluding secondary/duplicate alignments) is routine.

Reference: Li H. et al. (2009). “The Sequence Alignment/Map format and SAMtools.” Bioinformatics 25(16):2078–2079; SAM/BAM format specification maintained by the samtools/hts-specs GitHub repository (current spec at github.com/samtools/hts-specs).

samtools sort -@8 -o sorted.bam aligned.sam
samtools index sorted.bam        # produces sorted.bam.bai

2d. Mark duplicates

PCR amplification during library prep can copy the same original DNA fragment many times before sequencing; because these copies all map to the identical start/end position, they’d otherwise look like independent confirming evidence for a variant. MarkDuplicates identifies read pairs sharing identical 5′ mapping coordinates and flags all but one as duplicates (FLAG bit 1024) — it does not delete them, so the evidence is still auditable.

gatk MarkDuplicatesSpark -I sorted.bam -O dedup.bam -M dedup_metrics.txt

Reference: Picard toolkit, Broad Institute, broadinstitute.github.io/picard/; Van der Auwera G.A. & O’Connor B.D. (2020). Genomics in the Cloud. O’Reilly Media (Chapter on Data Pre-processing).

2e. Base Quality Score Recalibration (BQSR)

Sequencers have systematic, reproducible error biases (e.g., quality drifting near cycle end, or specific sequence-context motifs being harder to call correctly) that raw Phred scores don’t fully capture. BQSR builds an empirical error model by comparing observed mismatches against a “known-sites” VCF (positions already known to be common human variation, e.g., dbSNP, so those mismatches are excluded from the “error” estimate) and re-calibrates every quality score in the BAM accordingly.

gatk BaseRecalibrator -I dedup.bam -R GRCh38.fa \
   --known-sites dbsnp.vcf.gz --known-sites Mills_indels.vcf.gz -O recal.table
gatk ApplyBQSR -I dedup.bam -R GRCh38.fa --bqsr-recal-file recal.table -O analysis_ready.bam

Reference: DePristo M.A. et al. (2011). “A framework for variation discovery and genotyping using next-generation DNA sequencing data.” Nature Genetics 43:491–498.

What every bioinformatician should know about BAM QC: always check mean coverage depth (samtools depth/mosdepth), % properly paired reads, insert size distribution, and duplication rate (Picard’s dedup_metrics.txt) before trusting a callset — e.g., WGS typically targets ≥30× coverage, WES ≥50–100× on target, and duplication rates above ~20–30% usually indicate low library input.

File journey: .fastq.gz → .sam → .bam (sorted, .bai indexed) → .bam (dedup) → .bam (BQSR, “analysis-ready”)

3. Variant Calling

3a. Germline calling — GATK HaplotypeCaller

Rather than calling variants position-by-position, HaplotypeCaller performs local de novo reassembly: for each active region (where evidence of variation exists), it builds a De Bruijn graph of all overlapping reads, walks the graph to enumerate candidate haplotypes, then re-aligns each read against every candidate haplotype to compute genotype likelihoods via a pair-HMM (Hidden Markov Model). This local-assembly approach correctly resolves nearby/overlapping indels that simple pileup-based callers often get wrong.

gatk HaplotypeCaller -R GRCh38.fa -I analysis_ready.bam -O sample.g.vcf.gz -ERC GVCF

Reference: Poplin R. et al. (2018). “Scaling accurate genetic variant discovery to tens of thousands of samples.” bioRxiv 201178 (describes GATK4’s HaplotypeCaller pipeline in full).

DeepVariant takes a fundamentally different approach: it renders the read pileup at each candidate site as an RGB image (encoding base identity, quality, and strand as pixel channels) and classifies genotype using a convolutional neural network trained on millions of labeled examples (e.g., Genome in a Bottle truth sets) — it often outperforms statistical callers on indels. Reference: Poplin R. et al. (2018). “A universal SNP and small-indel variant caller using deep neural networks.” Nature Biotechnology 36:983–987.

3b. Joint genotyping (cohorts)

gatk GenomicsDBImport --genomicsdb-workspace-path db -V s1.g.vcf.gz -V s2.g.vcf.gz -L intervals.list
gatk GenotypeGVCFs -R GRCh38.fa -V gendb://db -O cohort.vcf.gz

Calling many samples jointly (rather than merging separately-called VCFs) is important because it lets the model borrow statistical strength across samples — a site that’s ambiguous in one sample alone becomes confidently callable if many other samples in the cohort clearly show the same variant.

3c. Filtering

Hard-filtering thresholds recommended by GATK Best Practices include: QD (Quality by Depth) < 2.0, FS (Fisher Strand bias) > 60.0 for SNPs, MQ (Mapping Quality) < 40.0, ReadPosRankSum < -8.0. Alternatively, VQSR trains a Gaussian mixture model on known-true sites (HapMap, Omni, 1000 Genomes) to learn a multi-dimensional quality surface, then scores every variant by how “true-like” it looks.

gatk VariantFiltration -V cohort.vcf.gz --filter-expression "QD<2.0" --filter-name "lowQD" -O filtered.vcf.gz

Reference: GATK Best Practices workflow documentation, Broad Institute, gatk.broadinstitute.org.

3d. Somatic (tumor-normal) calling

Mutect2 models the tumor sample as a mixture of a normal-genotype component plus somatic-variant-allele-fraction component, explicitly accounting for tumor purity/ploidy and sequencing artifacts (e.g., orientation bias from FFPE-preserved tissue), and subtracts anything also seen in the matched normal sample or a Panel of Normals (PoN) built from unrelated healthy samples to remove recurrent artifacts.

gatk Mutect2 -R GRCh38.fa -I tumor.bam -I normal.bam -normal normal_sample_name \
  --germline-resource gnomad.vcf.gz -O somatic.vcf.gz

Reference: Cibulskis K. et al. (2013). “Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples.” Nature Biotechnology 31:213–219; Benjamin D. et al. (2019). “Calling Somatic SNVs and Indels with Mutect2.” bioRxiv 861054.

3e. Structural & copy-number variants

Manta detects large deletions/duplications/inversions/translocations by clustering discordant read pairs (mates mapping too far apart or in the wrong orientation) and split reads (a single read that partially aligns to two different locations, indicating a breakpoint). Reference: Chen X. et al. (2016). “Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications.” Bioinformatics 32(8):1220–1222. CNVkit infers copy-number gains/losses from normalized read-depth ratios across binned genomic windows, comparing a sample against a pooled reference of normals. Reference: Talevich E., Shain A.H., Botton T., Bastian B.C. (2016). “CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing.” PLOS Computational Biology 12(4):e1004873.

3f. VCF format — full anatomy

##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
#CHROM  POS      ID   REF  ALT  QUAL  FILTER  INFO                          FORMAT   sample1
chr7    140453136  .  A    T    850   PASS    DP=42;AF=0.5;AC=1;AN=2       GT:AD:DP:GQ  0/1:20,22:42:99

CHROM/POS: genomic coordinate (1-based).
REF/ALT: reference allele and the variant allele(s) observed.
QUAL: Phred-scaled probability the site is NOT a true variant (higher = more confident).
FILTER: PASS or the name of a failed filter.
INFO: site-level annotations (DP=total depth, AF=allele frequency, AC=allele count, AN=allele number, etc.) — a semicolon-separated key=value list.
FORMAT/sample columns: per-sample fields. GT = genotype (0/1 = heterozygous for the ALT allele, 1/1 = homozygous ALT, 0/0 = homozygous reference); AD = allele depth (ref,alt read counts); DP = per-sample depth; GQ = genotype quality.

Reference: Danecek P. et al. (2011). “The variant call format and VCFtools.” Bioinformatics 27(15):2156–2158 (original VCF spec paper); current formal spec maintained at github.com/samtools/hts-specs.

What every bioinformatician should know: always sanity-check a callset’s Ti/Tv ratio (transition-to-transversion ratio) — for whole-genome human data it should be ~2.0–2.1, and for whole-exome ~3.0–3.3 (coding regions are enriched for transitions); a value far off this range signals a systematic calling problem. Reference: DePristo M.A. et al. (2011), Nature Genetics 43:491–498 (defines Ti/Tv as a standard QC metric).

File journey: .bam → .g.vcf.gz (per-sample, all sites) → .vcf.gz (joint-genotyped) → .vcf.gz (filtered, final callset)

4. Annotation & Interpretation

VEP (Variant Effect Predictor) intersects each variant against Ensembl gene models to predict consequence (missense, synonymous, stop_gained, splice_donor, etc.), using Sequence Ontology terms, and can layer on protein-domain, conservation, and pathogenicity-prediction (SIFT, PolyPhen) annotations.

vep -i filtered.vcf.gz -o annotated.vcf --cache --offline --assembly GRCh38

Reference: McLaren W. et al. (2016). “The Ensembl Variant Effect Predictor.” Genome Biology 17:122.

gnomAD aggregates and jointly calls sequencing data from >800,000 individuals (as of gnomAD v4) to give a population allele frequency for nearly every possible variant — the default reference for “is this variant rare enough to plausibly cause a rare disease.” Reference: Chen S. et al. (2024) “A genomic mutational constraint map using variation in 76,156 human genomes” (gnomAD v4), Nature 625:92–100 (supersedes the original Karczewski K. et al. 2020, Nature 581:434–443 gnomAD v2/v3 paper).

ClinVar is a public NIH-hosted archive of variant-phenotype relationships submitted by clinical laboratories and researchers, with ACMG/AMP-standardized pathogenicity classifications (Pathogenic, Likely Pathogenic, VUS, Likely Benign, Benign). Reference: Landrum M.J. et al. (2018). “ClinVar: improving access to variant interpretations and supporting evidence.” Nucleic Acids Research 46(D1):D1062–D1067; Richards S. et al. (2015), “Standards and guidelines for the interpretation of sequence variants” (the ACMG/AMP framework), Genetics in Medicine 17:405–424.

MAF format (Mutation Annotation Format) is a simplified, one-row-per-mutation TSV widely used in cancer cohorts (e.g., TCGA), generated from VCF+VEP output via vcf2maf, containing columns like Hugo_Symbol, Variant_Classification, Tumor_Sample_Barcode. Reference: NCI Genomic Data Commons, MAF Format specification, docs.gdc.cancer.gov.

File journey: .vcf.gz (filtered) → annotated .vcf / .tsv / .maf

5. Bulk RNA-seq

Library prep note: poly-A selection captures only mature, polyadenylated mRNA (missing non-coding RNAs), while ribo-depletion removes the highly abundant ribosomal RNA (>80% of total cellular RNA) and retains a broader RNA population — the choice matters for what biology is even measurable. Reference: Illumina TruSeq Stranded mRNA / Total RNA Library Prep guides, Illumina.com.

5a. Splice-aware alignment

Because mature mRNA has introns removed, an RNA read spanning an exon-exon junction will not appear contiguous in genomic DNA — the aligner must be able to “skip” the intron. STAR solves this with a two-pass, suffix-array-based algorithm (“Maximal Mappable Prefix” search) that identifies splice junctions de novo and can also detect novel/unannotated junctions.

STAR --genomeDir star_index --readFilesIn trim_R1.fastq.gz trim_R2.fastq.gz \
  --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts --sjdbGTFfile annotation.gtf

CIGAR strings for spliced reads use the N operation (e.g., 50M2000N50M = a 50-base exon, a 2000-base intron skip, another 50-base exon). Reference: Dobin A. et al. (2013). “STAR: ultrafast universal RNA-seq aligner.” Bioinformatics 29(1):15–21.

Salmon/kallisto skip genome alignment entirely and instead pseudo-align/quasi-map reads directly against a transcriptome index, using a k-mer-based approach and an Expectation-Maximization algorithm to resolve reads that map ambiguously to multiple isoforms of the same gene — much faster, and directly outputs transcript abundance estimates.

salmon quant -i salmon_index -l A -1 trim_R1.fastq.gz -2 trim_R2.fastq.gz -o quant_out --validateMappings

Reference: Patro R., Duggal G., Love M.I., Irizarry R.A., Kingsford C. (2017). “Salmon provides fast and bias-aware quantification of transcript expression.” Nature Methods 14:417–419.

5b. GTF/GFF3 annotation format

chr7  HAVANA  gene  140719327  140924929  .  -  .  gene_id "ENSG00000157764"; gene_name "BRAF";
chr7  HAVANA  exon  140924566  140924703  .  -  .  gene_id "ENSG00000157764"; transcript_id "ENST00000288602"; exon_number "1";

9 tab-delimited fields: chromosome, source, feature type, start, end, score, strand, frame, and an attributes field of key-value pairs — this is the lookup map aligners and counters use to know where every gene/exon/transcript lies. Reference: GENCODE Consortium GTF specification, gencodegenes.org; Ensembl GFF3 spec, ensembl.org.

5c. Quantification and count matrix

featureCounts -p -a annotation.gtf -o counts.txt sample.sorted.bam

Output: a matrix where rows = genes (Ensembl IDs), columns = samples, values = raw integer read counts assigned to each gene’s exons. Reference: Liao Y., Smyth G.K., Shi W. (2014). “featureCounts: an efficient general purpose program for assigning sequence reads to genomic features.” Bioinformatics 30(7):923–930.

5d. Differential expression

RNA-seq counts follow an overdispersed distribution (variance > mean), which a standard Poisson model underestimates — DESeq2 and edgeR both use a Negative Binomial generalized linear model, and shrink per-gene dispersion estimates toward a fitted trend across all genes (borrowing information across genes to stabilize estimates for genes with few reads/low replication), then test each gene with a Wald test (DESeq2) or exact/quasi-likelihood test (edgeR), followed by Benjamini-Hochberg FDR correction for the tens of thousands of simultaneous tests.

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"))

Reference: Love M.I., Huber W., Anders S. (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15:550; Robinson M.D., McCarthy D.J., Smyth G.K. (2010). “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics 26(1):139–140.

What every bioinformatician should know: always check the PCA plot of samples before trusting DE results — if replicates from the same condition don’t cluster together, there’s a likely batch effect that needs to be modeled (e.g., adding a batch term to the DESeq2 design formula) or the experiment has a technical problem.

5e. Enrichment analysis

GSEA ranks all genes by a differential statistic (not just a significance cutoff) and asks whether members of a predefined gene set (e.g., a KEGG pathway) are non-randomly clustered toward the top/bottom of that ranking, using a Kolmogorov-Smirnov-like running-sum statistic — this avoids the arbitrary “significant gene list” cutoff and retains statistical power from genes that individually don’t reach significance but move together as a pathway. Reference: Subramanian A. et al. (2005). “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.” PNAS 102(43):15545–15550.

File journey: .fastq.gz → .bam (STAR) or quant.sf (Salmon) → gene × sample count matrix (.txt/.csv) → DE results table → enrichment report/plots

6. Single-Cell RNA-seq

6a. Raw read anatomy (10x v3 chemistry)

R1 (28 bp): 16 bp cell barcode + 12 bp UMI (Unique Molecular Identifier — tags each captured RNA molecule so PCR duplicate reads of that same molecule can be collapsed to one count).
R2 (~90 bp): actual cDNA sequence, aligned to the transcriptome.
I1/I2: sample index reads, used for demultiplexing pooled lanes back into individual samples. Reference: 10x Genomics, “Chromium Single Cell 3’ Reagent Kits v3 User Guide,” 10xgenomics.com.

6b. Barcode processing, alignment, counting

Cell Ranger (10x’s proprietary pipeline, wrapping STAR internally) or STARsolo (an open-source, faster reimplementation integrated directly into STAR) correct sequencing errors in barcodes against a known whitelist of valid barcodes (allowing 1 mismatch), align cDNA reads to the transcriptome, deduplicate by UMI (collapsing reads sharing a cell barcode + UMI + gene, since these almost certainly originated from PCR copies of one original molecule), and tally a final count per gene per cell.

STAR --soloType CB_UMI_Simple --soloCBwhitelist 3M-february-2018.txt \
  --soloUMIlen 12 --readFilesIn R2.fastq.gz R1.fastq.gz --genomeDir star_index

References: 10x Genomics Cell Ranger documentation, support.10xgenomics.com; Kaminow B., Yunusov D., Dobin A. (2021). “STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data.” bioRxiv 2021.05.05.442755.

6c. MatrixMarket / HDF5 output format

matrix.mtx (MatrixMarket sparse format — since most genes are not detected in most cells, storing only non-zero entries saves enormous space):

%%MatrixMarket matrix coordinate integer general
33538 5000 12000000     ← n_genes n_cells n_nonzero_entries
34 1 3                  ← gene_row 34, cell_column 1, count 3

paired with features.tsv (gene IDs/names/types, one per row) and barcodes.tsv (cell barcode sequences, one per column). The .h5 (HDF5) format bundles all three into one binary file with faster random access. Reference: NIST MatrixMarket format specification, math.nist.gov/MatrixMarket; 10x Genomics Feature-Barcode Matrices documentation.

6d. Empty droplet removal

Because droplet capture is imperfect, most droplets contain no real cell — just ambient RNA floating in the reagent mix — which still generates a small nonzero count profile. EmptyDrops models the ambient RNA profile from very-low-count barcodes, then uses a Dirichlet-multinomial statistical test to ask whether each candidate barcode’s profile deviates significantly from pure ambient contamination, controlling FDR across all tested barcodes. Reference: Lun A.T.L. et al. (2019). “EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data.” Genome Biology 20:63.

6e. Loading, QC filtering, doublets

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[(adata.obs.n_genes_by_counts > 200) & (adata.obs.pct_counts_mt < 10)]

High mitochondrial-read percentage is a proxy for cell damage/death — dying cells lose cytoplasmic mRNA through membrane rupture but retain the mitochondria’s own RNA transcripts, skewing the ratio. Reference: Wolf F.A., Angerer P., Theis F.J. (2018). “SCANPY: large-scale single-cell gene expression data analysis.” Genome Biology 19:15; Luecken M.D. & Theis F.J. (2019), “Current best practices in single-cell RNA-seq analysis: a tutorial,” Molecular Systems Biology 15:e8746 (community QC-threshold guidance).

Doublets (two cells captured in one droplet, appearing as one hybrid transcriptome) are detected by simulating artificial doublets from the real data and checking which real cells’ expression profiles resemble those simulated doublets more than they resemble genuine single cells. Reference: Wolock S.L., Lopez R., Klein A.M. (2019). “Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data.” Cell Systems 8(4):281–291.

6f. Normalization, dimensionality reduction, clustering

SCTransform models each gene’s counts with a regularized negative binomial regression against total sequencing depth per cell, removing the technical depth effect while preserving biological variance (an improvement over simple log-normalization, which can distort variance structure for genes with different mean expression). Reference: Hafemeister C., Satija R. (2019). “Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression.” Genome Biology 20:296.

After PCA reduces ~2,000 highly variable genes down to ~30–50 principal components, UMAP builds a k-nearest-neighbor graph in that reduced space and optimizes a 2D layout that preserves local neighborhood structure (which cells are close to which), using a fuzzy topological framework grounded in Riemannian geometry and algebraic topology. Reference: McInnes L., Healy J., Melville J. (2018). “UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction.” arXiv:1802.03426.

Leiden clustering partitions the same k-nearest-neighbor graph into communities by optimizing a modularity-like quality function, guaranteeing (unlike its predecessor Louvain) that every resulting cluster is actually internally well-connected. Reference: Traag V.A., Waltman L., van Eck N.J. (2019). “From Louvain to Leiden: guaranteeing well-connected communities.” Scientific Reports 9:5233.

6g. Cell type annotation & downstream

Marker gene detection (e.g., Scanpy’s rank_genes_groups, a Wilcoxon rank-sum test per cluster vs. rest) identifies genes uniquely enriched per cluster, which are then matched against known cell-type marker panels manually or via automated reference mapping. SingleR correlates each cell’s (or cluster’s) full expression profile against a labeled reference dataset of pure, sorted cell types and assigns the label of the best-correlating reference. Reference: Aran D. et al. (2019). “Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage.” Nature Immunology 20:163–172 (introduces SingleR).

RNA velocity (scVelo) estimates each cell’s future transcriptional trajectory from the ratio of unspliced (intron-containing, newly transcribed) to spliced (mature) mRNA counts per gene, since a high unspliced:spliced ratio indicates a gene that is actively being turned on. Reference: Bergen V., Lange M., Peidli S., Stassen S.V., Theis F.J. (2020). “Generalizing RNA velocity to transient cell states through dynamical modeling.” Nature Biotechnology 38:1408–1414.

File journey: .fastq.gz → raw count matrix (.mtx/.h5) → filtered (real-cell) matrix → .h5ad/.rds object → normalized → clustered → annotated object → figures/tables

7. Nextflow / nf-core Orchestration

Nextflow is a domain-specific language (built on Groovy/JVM) that represents a pipeline as a dataflow graph: each process block declares its inputs/outputs as typed “channels,” and Nextflow automatically figures out the dependency order, parallelizes independent samples, resumes from cache after a crash (-resume), and can execute the exact same pipeline definition on a laptop, an HPC cluster (SLURM/PBS), or cloud batch systems (AWS Batch, Google Cloud Life Sciences) by swapping the executor configuration only.

process ALIGN {
    input: tuple val(sample), path(fastq1), path(fastq2)
    output: tuple val(sample), path("${sample}.bam")
    script: "bwa mem -t ${task.cpus} ref.fa ${fastq1} ${fastq2} | samtools sort -o ${sample}.bam"
}

Reference: Di Tommaso P. et al. (2017). “Nextflow enables reproducible computational workflows.” Nature Biotechnology 35:316–319.

nf-core provides community-vetted, continuously tested Nextflow pipelines with standardized structure, containerization (every process runs inside a versioned Docker/Singularity/Conda environment, so results are bit-for-bit reproducible years later), and automated documentation. - nf-core/sarek — germline & somatic variant calling (Sections 2–4). - nf-core/rnaseq — bulk RNA-seq (Section 5). - nf-core/scrnaseq — single-cell RNA-seq (Section 6), wrapping Cell Ranger/STARsolo/alevin-fry.

Reference: Ewels P.A. et al. (2020). “The nf-core framework for community-curated bioinformatics pipelines.” Nature Biotechnology 38:276–278.

Full File Format Reference Table

Format	Full name	Stage	Structure
`.fastq.gz`	FASTQ, gzip-compressed	Raw/trimmed reads	4 lines/read: ID, sequence, `+`, Phred quality
`.sam`/`.bam`	Sequence Alignment/Map (text/binary)	Aligned reads	11 mandatory tab-delimited fields per read; `.bai` companion index
`.cram`	Compressed reference-oriented alignment map	Aligned reads (space-efficient)	Like BAM but stores only differences from the reference genome
`.g.vcf.gz`	genomic VCF	Per-sample variant evidence	Every site, including confident non-variant (“reference”) blocks
`.vcf.gz`	Variant Call Format	Final variant calls	CHROM/POS/REF/ALT/QUAL/FILTER/INFO + per-sample genotype columns
`.maf`	Mutation Annotation Format	Cancer variant summary	One row per mutation, simplified gene/consequence columns
`.gtf`/`.gff3`	Gene Transfer Format / General Feature Format v3	Reference annotation	9 tab-delimited fields per gene/exon/transcript feature
`.mtx`	MatrixMarket sparse matrix	scRNA-seq raw counts	header + nonzero (row, col, value) triplets
`.h5`/`.h5ad`	Hierarchical Data Format 5 / AnnData-on-HDF5	scRNA-seq matrices/analysis objects	Binary, hierarchical, supports partial/random-access reads
`.rds`	R Data Serialization	Seurat analysis object	Binary serialized R object
`.bed`	Browser Extensible Data	Genomic intervals (e.g., exome capture regions)	chrom, chromStart (0-based), chromEnd

Complete Reference List

Andrews S. (2010). FastQC. Babraham Bioinformatics.
Cock P.J.A. et al. (2010). The Sanger FASTQ file format. Nucleic Acids Research 38(6):1767–1771.
Ewing B. & Green P. (1998). Base-calling of automated sequencer traces using phred. Genome Research 8:186–194.
Chen S. et al. (2018). fastp. Bioinformatics 34(17):i884–i890.
Ewels P. et al. (2016). MultiQC. Bioinformatics 32(19):3047–3048.
Li H. (2013). Aligning sequence reads with BWA-MEM. arXiv:1303.3997.
Li H. (2018). Minimap2. Bioinformatics 34(18):3094–3100.
Li H. et al. (2009). SAMtools. Bioinformatics 25(16):2078–2079.
DePristo M.A. et al. (2011). GATK framework. Nature Genetics 43:491–498.
Van der Auwera G.A. & O’Connor B.D. (2020). Genomics in the Cloud. O’Reilly.
Poplin R. et al. (2018). HaplotypeCaller scaling. bioRxiv 201178.
Poplin R. et al. (2018). DeepVariant. Nature Biotechnology 36:983–987.
Danecek P. et al. (2011). VCF and VCFtools. Bioinformatics 27(15):2156–2158.
Cibulskis K. et al. (2013). MuTect. Nature Biotechnology 31:213–219.
Benjamin D. et al. (2019). Mutect2. bioRxiv 861054.
Chen X. et al. (2016). Manta. Bioinformatics 32(8):1220–1222.
Talevich E. et al. (2016). CNVkit. PLOS Comp Biol 12(4):e1004873.
McLaren W. et al. (2016). VEP. Genome Biology 17:122.
Chen S. et al. (2024). gnomAD v4. Nature 625:92–100.
Landrum M.J. et al. (2018). ClinVar. Nucleic Acids Research 46(D1):D1062–D1067.
Richards S. et al. (2015). ACMG/AMP standards. Genetics in Medicine 17:405–424.
Dobin A. et al. (2013). STAR. Bioinformatics 29(1):15–21.
Patro R. et al. (2017). Salmon. Nature Methods 14:417–419.
Liao Y., Smyth G.K., Shi W. (2014). featureCounts. Bioinformatics 30(7):923–930.
Love M.I., Huber W., Anders S. (2014). DESeq2. Genome Biology 15:550.
Robinson M.D., McCarthy D.J., Smyth G.K. (2010). edgeR. Bioinformatics 26(1):139–140.
Subramanian A. et al. (2005). GSEA. PNAS 102(43):15545–15550.
Zheng G.X.Y. et al. (2017). 10x droplet method. Nature Communications 8:14049.
Kaminow B., Yunusov D., Dobin A. (2021). STARsolo. bioRxiv 2021.05.05.442755.
Lun A.T.L. et al. (2019). EmptyDrops. Genome Biology 20:63.
Wolf F.A., Angerer P., Theis F.J. (2018). Scanpy. Genome Biology 19:15.
Luecken M.D. & Theis F.J. (2019). scRNA-seq best practices tutorial. Molecular Systems Biology 15:e8746.
Wolock S.L., Lopez R., Klein A.M. (2019). Scrublet. Cell Systems 8(4):281–291.
Hafemeister C., Satija R. (2019). SCTransform. Genome Biology 20:296.
McInnes L., Healy J., Melville J. (2018). UMAP. arXiv:1802.03426.
Traag V.A., Waltman L., van Eck N.J. (2019). Leiden algorithm. Scientific Reports 9:5233.
Aran D. et al. (2019). SingleR. Nature Immunology 20:163–172.
Bergen V. et al. (2020). scVelo. Nature Biotechnology 38:1408–1414.
Hao Y. et al. (2021). Seurat v4. Cell 184(13):3573–3587.
Di Tommaso P. et al. (2017). Nextflow. Nature Biotechnology 35:316–319.
Ewels P.A. et al. (2020). nf-core. Nature Biotechnology 38:276–278.
Wang Y. et al. (2021). Nanopore sequencing. Nature Biotechnology 39:1348–1365.
Wenger A.M. et al. (2019). PacBio HiFi. Nature Biotechnology 37:1155–1162.
Nurk S. et al. (2022). T2T-CHM13 genome. Science 376:44–53.

Note: exact commands shown are illustrative of standard/documented usage — flag names and defaults change between tool versions, so always check current --help/official docs before running in production.

PART 2 — Interview-Style Q&A, Tumor Microenvironment Case Study, and a Buildable Nextflow DSL2 Pipeline

A note on sourcing for this section: I could not verify specific proprietary interview question sets from “HermoScience,” “microA1,” or Mercor-hosted assessments, and I won’t attribute fabricated questions to named companies. What follows instead is a rigorous, referenced Q&A set built from the concepts that public interview-prep guides (Glassdoor’s aggregated Bioinformatics interview data, compbiojobs.com’s guide covering Genentech/Illumina/Amgen-style processes, and ABRF/ASGCT career-center guides) confirm are actually tested in industry bioinformatics interviews — i.e., the same core competencies, asked as concrete solved problems rather than vague trivia.

A. Foundational / Statistics Questions (solved)

Q1. You ran 20,000 differential expression tests. Why can’t you just use p < 0.05? At α = 0.05, you’d expect ~1,000 false positives (20,000 × 0.05) by chance alone even if nothing were truly different — this is the multiple testing problem. The Bonferroni correction (α/n) controls the family-wise error rate but is very conservative for genomics. The standard approach is the Benjamini-Hochberg (BH) procedure, which controls the False Discovery Rate (expected proportion of false positives among your significant results, not among all tests) — it sorts p-values ascending, finds the largest k such that p(k) ≤ (k/n)×α, and calls everything up to k significant. This is why DESeq2/edgeR output a padj column (BH-adjusted) rather than raw p-values, and why bioinformaticians filter on padj < 0.05, not pvalue < 0.05. Reference: Benjamini Y., Hochberg Y. (1995). “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society Series B 57(1):289–300.

Q2. Why does RNA-seq use a Negative Binomial model instead of Poisson? Under Poisson, variance = mean. Real RNA-seq counts, even among biological replicates of the “same” condition, show variance > mean (“overdispersion”) because of biological variability between individuals/samples on top of pure sampling noise. The Negative Binomial adds a dispersion parameter to explicitly model this excess variance; DESeq2 further “shrinks” per-gene dispersion estimates toward a mean-dispersion trend fit across all genes, which stabilizes estimates for genes with few replicates. Reference: Love M.I., Huber W., Anders S. (2014). Genome Biology 15:550; Anders S., Huber W. (2010). “Differential expression analysis for sequence count data.” Genome Biology 11:R106 (introduces the overdispersion argument).

Q3. Given a VCF with GT=0/1, AD=20,22, DP=42, GQ=99 — walk through what each field tells you and whether you trust this call. - GT=0/1: heterozygous — one reference allele, one alternate allele copy. - AD=20,22: 20 reads supported the reference allele, 22 supported the alternate — nearly a clean 50/50 split, consistent with true heterozygosity (if it were skewed like 38,4, you’d suspect an artifact or somatic mosaicism rather than a clean germline het). - DP=42: total depth at this site — moderate; below typical WGS 30× consensus recommendation you’d want ≥10–15 reads minimum to call confidently, so 42 is comfortable. - GQ=99: genotype quality, Phred-scaled — 99 is the practical ceiling GATK reports, meaning extremely high confidence this specific genotype (over the alternatives 0/0 or 1/1) is correct. Conclusion: this is a high-confidence heterozygous call. Reference: VCF specification, hts-specs (github.com/samtools/hts-specs); GATK documentation on genotype fields, gatk.broadinstitute.org.

Q4. Explain Ti/Tv ratio and why it’s used as a QC metric. Transitions (purine↔︎purine: A↔︎G, or pyrimidine↔︎pyrimidine: C↔︎T) are chemically more likely to occur than transversions (purine↔︎pyrimidine) due to mutational mechanisms (e.g., deamination of methylated cytosine preferentially creates C→T transitions). True human variation therefore shows a consistent Ti/Tv ratio (~2.0–2.1 genome-wide, ~3.0–3.3 in exomes because coding-region selection further favors transitions). Random sequencing/alignment errors are not biased this way (~0.5 expected ratio), so a callset’s Ti/Tv dropping toward 1.0–1.5 signals contamination by false-positive calls. Reference: DePristo M.A. et al. (2011). Nature Genetics 43:491–498.

B. Alignment/Pipeline Design Questions (solved)

Q5. A colleague ran BWA-MEM without setting the -R read-group tag. What breaks downstream, and why? GATK’s MarkDuplicates, BQSR, and HaplotypeCaller all rely on the @RG header (specifically SM for sample name) to know which reads belong to which sample/library/lane. Without it, multi-sample or multi-lane BAMs can’t be correctly deduplicated per-library, BQSR can’t build a correct per-sample error model, and joint genotyping tools may mislabel or merge samples incorrectly. The fix is trivial (-R "@RG\tID:...\tSM:...\tPL:ILLUMINA" at alignment time) but re-running alignment on terabytes of data after the fact is expensive — a classic “cheap to prevent, expensive to fix” pipeline design lesson. Reference: GATK Best Practices, “read groups” documentation, gatk.broadinstitute.org.

Q6. Your WGS sample shows 45% duplication rate. Diagnose and recommend next steps. Typical WGS libraries show 5–20% duplication. 45% strongly suggests low library complexity — usually caused by too little input DNA going into library prep, requiring excessive PCR amplification cycles to generate enough material, which copies the same original fragments repeatedly. This is a wet-lab problem, not something a bioinformatics pipeline can fix after the fact (MarkDuplicates flags but can’t recover lost complexity/coverage). Recommendation: check input DNA quantity/quality (e.g., via Qubit/Bioanalyzer metrics if available) for that sample, and consider re-prepping the library with more input material rather than just re-sequencing deeper (deeper sequencing of a low-complexity library mostly just re-sequences the same duplicates). Reference: Picard MarkDuplicates documentation; Van der Auwera G.A. & O’Connor B.D. (2020), Genomics in the Cloud, O’Reilly (QC chapter).

Q7. Why does GATK recommend joint genotyping across a cohort rather than calling each sample separately and merging VCFs afterward? Per-sample calling with a naive union/merge afterward creates a serious bias: if sample A has strong evidence for a real variant but sample B has only weak, sub-threshold evidence at the same site, separate calling might call it confidently in A and simply not report it (rather than reporting 0/0) in B — leaving a misleading gap (“no-call” looks identical to “confidently reference”) in the merged VCF. GVCF + joint genotyping (GenomicsDBImport/GenotypeGVCFs) instead retains likelihood information at every site for every sample regardless of whether a variant was called, so the joint model can correctly borrow statistical evidence across the cohort and produce an honest 0/0 vs. missing distinction. Reference: Poplin R. et al. (2018). bioRxiv 201178.

C. Tumor Microenvironment (TME) Single-Cell Case Study — fully worked

Scenario: You’ve received 10x Genomics scRNA-seq data from a tumor biopsy (say, a non-small-cell lung cancer resection) and a matched adjacent-normal tissue sample. The clinical question: characterize the composition and state of immune infiltration in the tumor microenvironment.

Q8. What cell types would you expect to find in a solid tumor’s microenvironment, and what marker genes define them? The TME classically comprises: malignant epithelial/tumor cells (marked by tissue-specific epithelial markers plus copy-number-inferred aneuploidy — see Q10), T cells (CD3D/CD3E/CD3G pan-T; CD8A/CD8B for cytotoxic; CD4/IL7R for helper; FOXP3/IL2RA for regulatory T cells, Tregs), B cells (MS4A1/CD79A), plasma cells (MZB1/JCHAIN/SDC1), NK cells (NKG7/GNLY/KLRD1, CD3-negative), myeloid cells — tumor-associated macrophages (TAMs: CD68/CD163/CSF1R, often further split into M1-like pro-inflammatory vs. M2-like immunosuppressive states, though this binary is now considered an oversimplification), dendritic cells (LAMP3/CLEC9A), and stromal cells — cancer-associated fibroblasts (CAFs: COL1A1/PDGFRB/FAP) and endothelial cells (PECAM1/VWF). Reference: Tirosh I. et al. (2016). “Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.” Science 352:189–196 (foundational TME scRNA-seq marker panel paper); Chen Z. et al. (2021), “Single-cell RNA sequencing highlights the role of tumor-associated macrophages…,” Journal for ImmunoTherapy of Cancer 9:e001233 (TAM marker/state discussion).

Q9. How do you distinguish a genuine biological cell state from a technical batch effect when comparing tumor vs. adjacent-normal samples processed on different 10x lanes/days? Batch effects manifest as samples separating primarily along technical covariates (processing date, lane, operator) rather than biological ones in an unsupervised embedding (PCA/UMAP colored by batch, not just by expected cell type). The fix is integration: methods like Harmony iteratively adjust each cell’s PCA embedding to remove batch-associated variance while preserving cell-type structure (using a soft, iterative clustering + correction algorithm), or Seurat’s canonical correlation analysis (CCA)/reciprocal PCA (RPCA) anchor-based integration, which finds pairs of mutually-matched cells (“anchors”) across batches to align the datasets before joint clustering. A useful sanity check: after integration, a UMAP should show each expected cell type forming one cluster that contains cells from both batches mixed together, not two batch-specific sub-clusters of “the same” cell type. Reference: Korsunsky I. et al. (2019). “Fast, sensitive and accurate integration of single-cell data with Harmony.” Nature Methods 16:1289–1296; Stuart T. et al. (2019). “Comprehensive Integration of Single-Cell Data.” Cell 177(7):1888–1902 (Seurat v3 CCA integration).

Q10. Malignant epithelial cells often don’t have a clean unique marker — how do you confidently identify the tumor cell cluster vs. normal epithelium? Beyond marker genes (which malignant cells often still express, since they’re derived from the normal tissue), the standard approach is inferring copy number variation from expression data: tools like inferCNV or CopyKAT average expression across sliding windows of genomically-adjacent genes, using non-malignant cells (immune cells, which are reliably diploid) as a reference baseline — a chromosomal region that’s been amplified in the tumor genome shows systematically elevated average expression across that whole window (because more gene copies → more transcript), and deleted regions show systematically reduced expression. This produces a “CNV profile” per cell, and cells with large-scale aneuploidy patterns (versus the flat, quiet profile of the diploid reference immune cells) are called malignant. Reference: Tirosh I. et al. (2016), Science 352:189–196 (introduces expression-based CNV inference for tumor cell identification); Gao R. et al. (2021), “Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes,” Nature Biotechnology 39:599–608 (CopyKAT).

Q11. How would you quantify and statistically test whether Tregs are significantly enriched in the tumor vs. adjacent-normal sample? First, compute the proportion of Tregs out of total cells (or out of total T cells) per sample. Because compositional data (proportions across cell types within a sample sum to 1) can create spurious correlations if analyzed with naive tests, dedicated compositional-analysis tools are preferred over a plain t-test/chi-square on raw proportions. scCODA models cell-type counts per sample with a Bayesian Dirichlet-multinomial framework specifically designed to avoid this bias; a simpler and very common approach is a generalized linear mixed model (or Fisher’s exact test on counts) per cell type comparing condition groups, with FDR correction across the number of cell types tested. If you have only one tumor and one normal sample (no biological replicates), be honest about the fact that you cannot make a statistically powered population-level claim — you can describe the observed difference, but true replication (multiple patients) is required before it becomes a defensible enrichment claim. Reference: Büttner M. et al. (2021). “scCODA is a Bayesian model for compositional single-cell data analysis.” Nature Communications 12:6876.

Q12. What is “cell-cell communication” inference and how does it work mechanistically? Tools like CellChat and CellPhoneDB use a curated database of known ligand-receptor pairs (e.g., PD-L1/CD274 on tumor or myeloid cells binding PD-1/PDCD1 on exhausted T cells), then, for every pair of cell-type clusters, test whether the ligand is significantly expressed in the “sender” cluster and the receptor significantly expressed in the “receiver” cluster simultaneously (often via a permutation test that shuffles cell-type labels to build a null distribution). A statistically significant, highly-expressed PD-L1(tumor)→PD-1(T cell) interaction inferred this way is a computational hypothesis consistent with immune checkpoint-mediated T cell exhaustion in that tumor — directly relevant to predicting response to checkpoint-inhibitor immunotherapy, though it remains a hypothesis requiring functional/clinical validation, not proof of an active signaling event. Reference: Jin S. et al. (2021). “Inference and analysis of cell-cell communication using CellChat.” Nature Communications 12:1088; Efremova M. et al. (2020), “CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes,” Nature Protocols 15:1484–1506.

Q13. T cell “exhaustion” is a key TME concept — how is it identified computationally? Exhausted CD8 T cells (a dysfunctional state arising from chronic antigen stimulation in the tumor, as opposed to healthy effector/memory T cells) are identified by co-expression of inhibitory checkpoint receptor genes — PDCD1 (PD-1), HAVCR2 (TIM-3), LAG3, CTLA4, TIGIT — often combined into an “exhaustion score” (e.g., via Scanpy’s score_genes, which averages expression of a gene set against a randomly sampled control gene set to account for baseline expression level), alongside reduced expression of effector cytokines (IFNG, GZMB) at the most terminally exhausted end of the spectrum. This is typically visualized along a pseudotime trajectory (Monocle3) from naive/effector through progenitor-exhausted to terminally-exhausted states. Reference: Wherry E.J., Kurachi M. (2015). “Molecular and cellular insights into T cell exhaustion.” Nature Reviews Immunology 15:486–499 (defines the exhaustion marker panel); Miller B.C. et al. (2019), “Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade,” Nature Immunology 20:326–336 (progenitor vs. terminal exhaustion states).

D. Coding / Algorithmic Questions (as seen in general bioinformatics technical screens)

Q14. Given a sorted BAM, write pseudocode to compute per-base coverage depth across a region without loading the whole BAM into memory.

# Streaming approach using a "sweep line" / active-interval-count algorithm — O(n) in reads, O(1) extra memory per position
import pysam
bam = pysam.AlignmentFile("analysis_ready.bam", "rb")
for pileupcolumn in bam.pileup("chr7", 140453000, 140453500, truncate=True):
    depth = pileupcolumn.n   # number of reads overlapping this position
    print(pileupcolumn.pos, depth)

This uses pysam’s pileup engine (a Python wrapper around htslib, the same C library underlying samtools), which internally streams reads in coordinate order and maintains a sliding window of “currently overlapping” reads — conceptually identical to the classic “meeting rooms” / interval-overlap-counting sweep-line algorithm taught in general coding interviews, applied to genomic intervals instead of time intervals. Reference: htslib/pysam documentation, pysam.readthedocs.io; the underlying algorithmic pattern is standard interval-scheduling (“sweep line”) as described in Cormen T.H. et al., Introduction to Algorithms, 3rd ed., MIT Press, 2009.

Q15. What’s the time/space complexity advantage of BWA’s FM-index (Burrows-Wheeler Transform) over a naive substring search when aligning millions of reads to a 3-billion-base genome? A naive search for each read against the genome is O(genome_length × read_length) per read — completely infeasible at scale. The Burrows-Wheeler Transform, combined with an FM-index (a compressed suffix array supporting fast rank/select queries), allows exact substring search in roughly O(read_length) time regardless of genome size, using backward search that narrows a range in the suffix array one query character at a time. This is the same core data structure used in bzip2 compression, repurposed for genomics; BWA-MEM extends this to allow approximate (mismatch/indel-tolerant) matching via seed-and-extend. Reference: Li H., Durbin R. (2009). “Fast and accurate short read alignment with Burrows-Wheeler transform.” Bioinformatics 25(14):1754–1760 (the original BWA paper, predating BWA-MEM); Burrows M., Wheeler D.J. (1994), “A block-sorting lossless data compression algorithm,” Digital Equipment Corporation Technical Report 124 (the original BWT).

E. Building a Nextflow DSL2 Pipeline — Complete, Replicable, Step-by-Step

This builds a minimal but fully functional germline variant-calling pipeline (FASTQ → analysis-ready BAM → VCF) using Nextflow DSL2, modeled on the same logical structure as nf-core/sarek. DSL2 (the current Nextflow syntax, replacing the deprecated DSL1) organizes code into reusable, independently-testable process and workflow blocks connected by channels. Reference: Nextflow DSL2 documentation, nextflow.io/docs/latest/dsl2.html; Di Tommaso P. et al. (2017), Nature Biotechnology 35:316–319.

Step 1 — Install Nextflow and a container engine

curl -s https://get.nextflow.io | bash
chmod +x nextflow && sudo mv nextflow /usr/local/bin/
nextflow -version                    # confirm ≥23.x, which defaults to DSL2
# also install Docker or Singularity/Apptainer — containers are what make the pipeline reproducible

Reference: Nextflow installation docs, nextflow.io/docs/latest/install.html.

Step 2 — Project structure

my-variant-pipeline/
├── main.nf                # entrypoint workflow
├── nextflow.config         # resource/container/profile configuration
├── modules/
│   ├── fastqc.nf
│   ├── fastp.nf
│   ├── bwa_align.nf
│   ├── mark_duplicates.nf
│   ├── bqsr.nf
│   └── haplotype_caller.nf
└── data/
    └── samplesheet.csv      # sample_id,fastq_1,fastq_2

This modular layout (one process per file under modules/) mirrors nf-core’s convention and lets each module be unit-tested independently. Reference: nf-core pipeline template/module structure, nf-co.re/docs/contributing/modules.

Step 3 — Define a module (example: `modules/fastqc.nf`)

process FASTQC {
    tag "$sample_id"
    container 'biocontainers/fastqc:v0.11.9_cv8'
    publishDir "results/fastqc", mode: 'copy'

    input:
    tuple val(sample_id), path(reads)

    output:
    tuple val(sample_id), path("*.zip"), path("*.html")

    script:
    """
    fastqc -t ${task.cpus} ${reads}
    """
}

tag labels this task in logs/reports with the sample ID for readability.
container pins the exact software version — this is what guarantees the same result on any machine, years later.
input/output declare typed channels; tuple val(sample_id), path(reads) means “a sample name paired with one or more files.”
publishDir copies designated outputs to a results folder (the working directory itself is a temporary hashed cache, not meant for direct browsing). Reference: Nextflow process reference documentation, nextflow.io/docs/latest/process.html.

Step 4 — Define the alignment module (`modules/bwa_align.nf`)

process BWA_ALIGN {
    tag "$sample_id"
    container 'staphb/bwa:0.7.17'
    cpus 8
    memory '16 GB'

    input:
    tuple val(sample_id), path(reads)
    path reference
    path reference_index   // .amb .ann .bwt .pac .sa, .fai, .dict — bundled

    output:
    tuple val(sample_id), path("${sample_id}.sorted.bam"), path("${sample_id}.sorted.bam.bai")

    script:
    """
    bwa mem -t ${task.cpus} -R "@RG\\tID:${sample_id}\\tSM:${sample_id}\\tPL:ILLUMINA" \\
        ${reference} ${reads[0]} ${reads[1]} | \\
        samtools sort -@ ${task.cpus} -o ${sample_id}.sorted.bam -
    samtools index ${sample_id}.sorted.bam
    """
}

cpus/memory directives declare per-process resource needs, which Nextflow passes to the executor (local, SLURM, AWS Batch) — this is how the same pipeline scales from a laptop to a cluster without code changes. Reference: Nextflow process resource directives, nextflow.io/docs/latest/process.html#resource-directives.

Step 5 — Wire modules together in `main.nf`

#!/usr/bin/env nextflow
nextflow.enable.dsl = 2

include { FASTQC }           from './modules/fastqc.nf'
include { FASTP }            from './modules/fastp.nf'
include { BWA_ALIGN }        from './modules/bwa_align.nf'
include { MARK_DUPLICATES }  from './modules/mark_duplicates.nf'
include { BQSR }             from './modules/bqsr.nf'
include { HAPLOTYPE_CALLER } from './modules/haplotype_caller.nf'

workflow {
    // Build a channel from the samplesheet: each row -> (sample_id, [fastq_1, fastq_2])
    Channel
        .fromPath(params.samplesheet)
        .splitCsv(header: true)
        .map { row -> tuple(row.sample_id, [file(row.fastq_1), file(row.fastq_2)]) }
        .set { reads_ch }

    reference = file(params.reference)
    ref_index = file(params.reference_index_dir)

    FASTQC(reads_ch)
    FASTP(reads_ch)
    BWA_ALIGN(FASTP.out, reference, ref_index)
    MARK_DUPLICATES(BWA_ALIGN.out)
    BQSR(MARK_DUPLICATES.out, reference, params.known_sites)
    HAPLOTYPE_CALLER(BQSR.out, reference)
}

include { PROCESS } from './modules/x.nf' is DSL2’s module system — this is the key feature DSL2 added over DSL1, enabling reuse of the same module across multiple pipelines. Channels (reads_ch, FASTP.out, etc.) are the “wires” connecting each process’s output directly to the next process’s input; Nextflow automatically figures out the correct execution order and parallelizes independent branches. Reference: Nextflow DSL2 include and module documentation, nextflow.io/docs/latest/module.html.

Step 6 — Configure execution in `nextflow.config`

params {
    samplesheet          = "data/samplesheet.csv"
    reference             = "data/GRCh38.fa"
    reference_index_dir   = "data/ref_index/"
    known_sites            = "data/dbsnp.vcf.gz"
    outdir                 = "results"
}

profiles {
    docker {
        docker.enabled = true
    }
    singularity {
        singularity.enabled = true
    }
    slurm {
        process.executor = 'slurm'
        process.queue    = 'normal'
    }
    aws {
        process.executor = 'awsbatch'
        process.queue    = 'my-batch-queue'
        aws.region        = 'us-east-1'
    }
}

process {
    withName: BWA_ALIGN {
        cpus   = 8
        memory = '16 GB'
    }
    withName: HAPLOTYPE_CALLER {
        cpus   = 4
        memory = '8 GB'
    }
}

profiles let you switch the entire execution environment (local Docker → university SLURM cluster → AWS cloud) with a single command-line flag, without touching any pipeline logic — this environment/logic separation is the central reproducibility feature of Nextflow. Reference: Nextflow configuration documentation, nextflow.io/docs/latest/config.html.

Step 7 — Run it

nextflow run main.nf -profile docker --samplesheet data/samplesheet.csv -resume

-profile docker selects the Docker execution profile defined above.
-resume re-uses Nextflow’s content-addressed cache (each task’s inputs are hashed; if unchanged, the cached result is reused instead of re-running) — critical for iterating on a pipeline without re-running every upstream step after every small fix. Reference: Nextflow CLI reference, nextflow.io/docs/latest/cli.html.

Step 8 — Inspect execution reports

nextflow run main.nf -profile docker -with-report report.html -with-trace trace.txt -with-timeline timeline.html -with-dag flowchart.png

These built-in flags generate: an HTML resource-usage report (CPU/memory/time per task — critical for right-sizing cluster resource requests), a machine-readable trace log, a Gantt-chart-style timeline of parallel execution, and an auto-generated diagram of the pipeline’s actual dependency graph. Reference: Nextflow tracing/reporting documentation, nextflow.io/docs/latest/tracing.html.

Step 9 — Testing with nf-test (production-grade practice)

// tests/fastqc.test.nf
nextflow_process {
    name "Test FASTQC"
    script "modules/fastqc.nf"
    process "FASTQC"

    test("Should run FastQC on paired-end reads") {
        when {
            process {
                """
                input[0] = [ 'sample1', [file('test_data/R1.fastq.gz'), file('test_data/R2.fastq.gz')] ]
                """
            }
        }
        then {
            assert process.success
            assert process.out.size() > 0
        }
    }
}

nf-test (adopted as the standard testing framework across nf-core) lets each module be validated in isolation with small test fixtures — this is what enables nf-core’s continuous integration to catch a broken module before it reaches production users. Reference: nf-test documentation, code.askimed.com/nf-test; nf-core module testing guidelines, nf-co.re/docs/contributing/tutorials/nf-test_assertions.

Step 10 — Replicate an existing, validated pipeline instead of writing from scratch

For real production use, the recommended practice is not to hand-roll every module but to run the already-validated, community-tested nf-core pipeline directly:

nextflow run nf-core/sarek -profile docker \
  --input samplesheet.csv --genome GATK.GRCh38 --tools haplotypecaller,vep --outdir results/

This single command reproduces the entire Sections 2–4 pipeline above (alignment through annotation), using exactly the tool versions and parameters validated by the nf-core community, with full provenance tracking (a pipeline_info folder records every software version and parameter used, satisfying reproducibility requirements for publication or clinical use). Reference: nf-core/sarek documentation, nf-co.re/sarek; Ewels P.A. et al. (2020), Nature Biotechnology 38:276–278.

Additional References for Part 2

Benjamini Y., Hochberg Y. (1995). JRSS B 57(1):289–300.
Anders S., Huber W. (2010). Genome Biology 11:R106.
Tirosh I. et al. (2016). Science 352:189–196.
Chen Z. et al. (2021). Journal for ImmunoTherapy of Cancer 9:e001233.
Korsunsky I. et al. (2019). Nature Methods 16:1289–1296.
Stuart T. et al. (2019). Cell 177(7):1888–1902.
Gao R. et al. (2021). Nature Biotechnology 39:599–608.
Büttner M. et al. (2021). Nature Communications 12:6876.
Jin S. et al. (2021). Nature Communications 12:1088.
Efremova M. et al. (2020). Nature Protocols 15:1484–1506.
Wherry E.J., Kurachi M. (2015). Nature Reviews Immunology 15:486–499.
Miller B.C. et al. (2019). Nature Immunology 20:326–336.
Li H., Durbin R. (2009). Bioinformatics 25(14):1754–1760.
Burrows M., Wheeler D.J. (1994). DEC Technical Report 124.
Cormen T.H., Leiserson C.E., Rivest R.L., Stein C. (2009). Introduction to Algorithms, 3rd ed. MIT Press.
Nextflow documentation, nextflow.io/docs/latest/.
nf-core documentation, nf-co.re.
nf-test documentation, code.askimed.com/nf-test.

PART 3 — Landmark Tumor Microenvironment Research, Deeper RNA-seq/scRNA-seq Algorithms, and Field History

Scope note: “everything from every lab/site” isn’t something I can responsibly claim to cover exhaustively — the TME field spans thousands of labs and papers. What I’ve done instead is cover the specific consortia, landmark studies, and algorithmic building blocks that a masters graduate would actually be expected to know by name in an interview, each with a real, checkable citation. Where I’m not confident a claim is accurate, I’ve left it out rather than filling the gap.

1. Landmark Consortia and Achievements in TME Research

The Human Tumor Atlas Network (HTAN) — an NCI-funded, multi-institutional consortium (launched 2018) building 3D, multi-omic, spatially-resolved atlases of human tumors as they evolve from precancerous lesions through advanced/metastatic disease and treatment resistance, explicitly integrating scRNA-seq, spatial transcriptomics, and imaging across ~10 major cancer types. Reference: Rozenblatt-Rosen O. et al. (2020). “The Human Tumor Atlas Network: Charting Tumor Transitions across Space and Time at Single-Cell Resolution.” Cell 181(2):236–249.

The Human Cell Atlas (HCA) — a global initiative (launched 2016) to map every cell type in the healthy human body as a reference, which TME studies rely on as the baseline for identifying what’s “abnormal” in a tumor’s cellular composition. Reference: Regev A. et al. (2017). “The Human Cell Atlas.” eLife 6:e27041.

Discovery of exhausted T cells and the mechanistic basis of checkpoint blockade — foundational work (originally in chronic viral infection models, later directly translated to the TME) established that chronic antigen exposure drives T cells into a distinct, epigenetically fixed dysfunctional state (rather than simple “fatigue”), which underlies why PD-1/PD-L1 checkpoint-blockade immunotherapy (e.g., pembrolizumab, nivolumab) works in some tumors and not others — this discovery, built substantially on immunology + computational profiling of T cell states, contributed to the 2018 Nobel Prize in Physiology or Medicine awarded to James P. Allison and Tasuku Honjo for the discovery of cancer therapy by inhibition of negative immune regulation. Reference: Wherry E.J. (2011). “T cell exhaustion.” Nature Immunology 12:492–499; The Nobel Prize in Physiology or Medicine 2018, NobelPrize.org.

Single-cell dissection of melanoma TME (Tirosh et al. 2016) — the study that essentially established the modern methodology described in Part 2C (Q10): applying scRNA-seq to freshly dissociated tumor biopsies across multiple patients, using expression-inferred CNV to separate malignant from non-malignant cells, and characterizing the immune/stromal composition — this became the template nearly all subsequent tumor scRNA-seq atlas papers follow. Reference: Tirosh I. et al. (2016). Science 352:189–196.

Spatial transcriptomics entering the TME field — technologies like 10x Genomics Visium (barcoded spots on a slide, capturing gene expression while preserving tissue location) and imaging-based methods (MERFISH, Xenium, CosMx) added back the spatial context that dissociated scRNA-seq loses (you know a cell’s expression profile but not literally where it sat in the tissue relative to other cells) — critical for TME work because spatial proximity of, e.g., exhausted T cells to tumor cells vs. being excluded to the tumor margin is itself informative about immune evasion mechanisms. Reference: Ståhl P.L. et al. (2016). “Visualization and analysis of gene expression in tissue sections by spatial transcriptomics.” Science 353:78–82 (foundational spatial transcriptomics method); Moses L., Pachter L. (2022), “Museum of spatial transcriptomics,” Nature Methods 19:534–546 (survey of the field’s major methods).

CAR-T cell therapy — while not itself a “TME bioinformatics” achievement, the clinical success of chimeric antigen receptor T cell therapy for blood cancers (FDA approval of tisagenlecleucel in 2017) created major downstream demand for TME-informed bioinformatics, because CAR-T efficacy in solid tumors (unlike blood cancers) is limited substantially by the immunosuppressive TME — this is one of the primary translational motivations driving current TME single-cell research funding. Reference: June C.H., Sadelain M. (2018). “Chimeric Antigen Receptor Therapy.” New England Journal of Medicine 379:64–73.

2. Deeper RNA-seq: the statistics and math a graduate should be able to derive, not just name

PCA, precisely. Principal Component Analysis finds the orthogonal directions of maximum variance in the (genes × samples) expression matrix. Formally: center the data, compute the covariance matrix Σ = (1/(n-1)) XᵀX, then eigendecompose Σ = VΛVᵀ. The eigenvectors V are the principal component directions; the eigenvalues Λ give the variance explained by each. In practice this is computed via Singular Value Decomposition (SVD) of X directly (X = UDVᵀ) for numerical stability, since it avoids explicitly forming the covariance matrix. A “scree plot” of eigenvalues (variance explained per PC) is the standard way to decide how many PCs (typically 10–50 for scRNA-seq) carry real biological signal before the remainder is just noise. Reference: Jolliffe I.T., Cadima J. (2016). “Principal component analysis: a review and recent developments.” Philosophical Transactions of the Royal Society A 374:20150202.

DESeq2’s shrinkage estimator, precisely. Rather than trusting a single gene’s raw maximum-likelihood log2 fold-change estimate (which is noisy for lowly-expressed genes or genes with few replicates), DESeq2 applies empirical Bayes shrinkage: it fits a prior distribution on log fold changes across all genes, then computes each gene’s posterior estimate as a weighted compromise between its own noisy MLE estimate and the cohort-wide prior — genes with low counts/high uncertainty get pulled harder toward zero (shrunk), while high-confidence, high-count genes are barely shrunk at all. This is what the lfcShrink() function does, and is why DESeq2’s shrunken LFC values are recommended for ranking/visualization over the raw log2FoldChange column. Reference: Love M.I., Huber W., Anders S. (2014). Genome Biology 15:550; Zhu A., Ibrahim J.G., Love M.I. (2019), “Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences,” Bioinformatics 35(12):2084–2092 (the apeglm shrinkage estimator now recommended by DESeq2).

Leiden clustering’s objective function, precisely. Leiden (like its predecessor Louvain) optimizes the Constant Potts Model (CPM) or modularity — modularity Q is defined as Q = (1/2m) Σᵢⱼ [Aᵢⱼ − (kᵢkⱼ)/2m] δ(cᵢ,cⱼ), where Aᵢⱼ is the edge weight between cells i and j in the k-nearest-neighbor graph, kᵢ is the degree of node i, m is total edge weight, and δ(cᵢ,cⱼ) = 1 if cells i,j are in the same cluster. In plain terms: it rewards putting cells in the same cluster if they’re more connected to each other than you’d expect by random chance given their overall connectivity. Leiden’s specific improvement over Louvain is a refinement step that guarantees every cluster returned is actually a single connected subgraph (Louvain could produce disconnected “clusters” as an artifact of its local-move heuristic). Reference: Traag V.A., Waltman L., van Eck N.J. (2019). Scientific Reports 9:5233; Newman M.E.J. (2006), “Modularity and community structure in networks,” PNAS 103(23):8577–8582 (defines modularity).

UMAP’s mathematical basis, briefly. UMAP constructs a weighted k-nearest-neighbor graph in high-dimensional (PCA) space, treats edge weights as fuzzy set membership strengths (grounded in a Riemannian-manifold assumption that data lies on a locally-uniform manifold), then optimizes a low-dimensional (2D) layout via stochastic gradient descent to minimize cross-entropy between the high-dimensional and low-dimensional fuzzy topological representations — the attractive/repulsive force balance is analogous to but mathematically distinct from t-SNE’s KL-divergence-based approach, and UMAP is generally faster and better preserves more global structure. Reference: McInnes L., Healy J., Melville J. (2018). arXiv:1802.03426.

3. Deeper scRNA-seq: additional algorithms/concepts graduate-level interviews probe

Ambient RNA correction (beyond just EmptyDrops). Even “real” cell barcodes contain some contamination from ambient RNA (lysed cells’ RNA floating freely in the droplet-generation reagent before encapsulation), which can create spurious low-level “expression” of genes that aren’t truly expressed in that cell (e.g., hemoglobin genes appearing at low levels in non-erythroid cells from a blood-contaminated sample). SoupX and CellBender model and subtract this ambient contamination profile per cell before downstream analysis. Reference: Young M.D., Behjati S. (2020). “SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data.” GigaScience 9(12):giaa151; Fleming S.J. et al. (2023), “Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender,” Nature Methods 20:1323–1335.

Pseudobulk analysis — a critical, often-tested concept. A common mistake is running standard bulk differential-expression tools (DESeq2) treating each cell as a replicate — this dramatically inflates statistical significance because cells from the same patient are not independent biological replicates (they’re pseudoreplicates); the correct approach is to sum/average counts per cell type per patient into a “pseudobulk” sample, then run DESeq2 with patients (not cells) as the unit of replication. Reference: Squair J.W. et al. (2021). “Confronting false discoveries in single-cell differential expression.” Nature Communications 12:5692 (the definitive paper demonstrating this pseudoreplication problem and its fix).

Multimodal / CITE-seq. Beyond RNA alone, CITE-seq simultaneously captures surface protein abundance (via DNA-barcoded antibodies, read out alongside the RNA library) and RNA in the same cell, giving a combined transcriptomic + proteomic readout — commonly used in immune/TME profiling because many key immune markers (e.g., PD-1 protein vs. PDCD1 mRNA) correlate imperfectly at the RNA level, and protein-level data is often more directly interpretable clinically. Reference: Stoeckius M. et al. (2017). “Simultaneous epitope and transcriptome measurement in single cells.” Nature Methods 14:865–868.

Single-cell TCR/BCR sequencing (VDJ). For T/B cell studies specifically (highly relevant to TME immuno-oncology), 10x’s paired scRNA + V(D)J sequencing links each cell’s transcriptome to its unique T cell receptor (or B cell receptor) sequence, enabling clonal expansion analysis — i.e., identifying which specific T cell clones have expanded in the tumor (a signature of an antigen-specific anti-tumor response) versus the polyclonal background repertoire. Reference: 10x Genomics, “Chromium Single Cell V(D)J Reagent Kits” technical documentation, 10xgenomics.com; Wu T.D. et al. (2020), “Peripheral T cell expansion predicts tumour infiltration and clinical response,” Nature 579:274–278 (a landmark application in melanoma checkpoint-blockade patients).

4. Pipeline flow diagram

PART 4 — Additional Depth: More Worked Questions, Figures, and RNA-seq/scRNA-seq Detail

Continuing from Parts 1-3 above. This section adds further worked interview problems, notes on accompanying figures, and extra mechanistic detail on RNA-seq/scRNA-seq that a masters-level graduate is expected to be able to explain from first principles, not just name-drop.

1. More Worked Interview Questions

Q16. You’re given raw counts for 3 tumor replicates and 3 normal replicates. edgeR gives you 500 DE genes at FDR 5%; DESeq2 gives you 350, with 300 overlapping. How do you explain the discrepancy to a PI who wants “the right answer”? Both tools model the same underlying biology (negative binomial counts) but differ in dispersion estimation strategy and default filtering: edgeR’s tagwise empirical-Bayes-moderated dispersion and DESeq2’s median-of-ratios normalization plus independent filtering (removing low-count genes before FDR correction, increasing power by reducing the number of tests) produce numerically different but broadly concordant results. There is no single “right answer” — the overlapping 300 genes are the highest-confidence calls; genes unique to one tool sit closer to that tool’s significance boundary and warrant closer inspection (effect size, raw counts, biological plausibility) rather than being dismissed or trusted outright. Reporting the intersection is standard practice when robustness matters more than maximizing the gene count. Reference: Soneson C., Delorenzi M. (2013). “A comparison of methods for differential expression analysis of RNA-seq data.” BMC Bioinformatics 14:91.

Q17. Explain why single-cell data needs a different normalization strategy than bulk RNA-seq, concretely. Bulk RNA-seq is normalized by library-size scaling factors (DESeq2’s median-of-ratios, or CPM) because sequencing depth differences between samples are the dominant technical confound, and each bulk sample already averages over thousands of cells. Single-cell data adds two problems bulk doesn’t have: (1) extreme sparsity — a cell often has zero detected counts for 80-90% of genes, largely due to “dropout” (the molecule wasn’t captured during very low-input reverse transcription, not necessarily biological absence), and (2) technical variability per cell, driven by capture efficiency differences between individual droplets. This is why SCTransform models variance as a function of mean expression per gene per cell via regularized regression, rather than one global size factor. Reference: Hafemeister C., Satija R. (2019). Genome Biology 20:296; Vallejos C.A. et al. (2017), “Normalizing single-cell RNA sequencing data: challenges and opportunities,” Nature Methods 14:565-571.

Q18. A candidate says “I removed batch effects by regressing out the batch variable in a linear model before clustering.” What’s the concern versus Harmony/CCA integration? Linear regression-based batch correction assumes the batch effect is a uniform additive shift affecting every gene the same way, regardless of cell type — but batch effects are often cell-type-specific and non-linear. If batch and biological condition are even partially confounded (e.g., all tumor samples in batch 1, all normal in batch 2), naive regression can accidentally regress out real biological signal along with the technical effect. Harmony/CCA-style integration instead works in the reduced embedding space using local cell-neighborhood structure to align batches, without assuming one global linear correction — more robust, though not a substitute for balanced experimental design. Reference: Tran H.T.N. et al. (2020). “A benchmark of batch-effect correction methods for single-cell RNA sequencing data.” Genome Biology 21:12.

Q19. Walk through, mathematically, what “UMI collapsing” is actually correcting for. Each captured RNA molecule gets a UMI (a random ~10-12bp barcode) attached before PCR amplification, which then copies that one original molecule into many PCR-duplicate reads sharing the identical UMI + cell barcode + gene mapping. Without correction, sequencing depth would directly inflate apparent gene expression. UMI collapsing counts unique (cell barcode, UMI, gene) combinations rather than raw reads, so 500 PCR-duplicate reads sharing one UMI count as 1 molecule — converting “read counting” (confounded by amplification bias) into genuine “molecule counting.” Reference: Islam S. et al. (2014). “Quantitative single-cell RNA-seq with unique molecular identifiers.” Nature Methods 11:163-166; Kivioja T. et al. (2012), Nature Methods 9:72-74.

2. Accompanying Figures

Two supporting diagrams accompany this document, rendered inline in the conversation (interactive SVG, not static image files, since a knitting environment can’t execute them): - Figure 1 — Pipeline branch overview, shown earlier: FASTQ QC/trimming forking into WGS/WES, bulk RNA-seq, and scRNA-seq branches, converging into Nextflow/nf-core orchestration. - Figure 2 — TME single-cell workflow, shown below this document: raw count matrix through to an annotated tumor microenvironment atlas.

If you want static image files embedded directly in a knitted PDF/HTML output (via knitr::include_graphics()), export any inline diagram as PNG/SVG and I can wire in the corresponding chunk referencing that file path — R Markdown itself cannot generate novel scientific figures without either your data (for ggplot2) or a pre-existing image file.

3. RNA-seq/scRNA-seq: filling in remaining mechanistic gaps

Exact alignment scoring in STAR. STAR performs a Maximal Mappable Prefix (MMP) search: starting from a read’s 5’ end, it extends the longest exact match against the genome’s suffix array; at a mismatch or splice junction, it seeds a new search from that point, then stitches seeds into one alignment, scored by mismatch/splice-junction penalties tuned to canonical GT-AG splice motifs. This seed-and-stitch design lets STAR detect novel splice junctions not present in the input GTF. Reference: Dobin A. et al. (2013). Bioinformatics 29(1):15-21.

Why Salmon’s EM step matters, concretely. A read entirely within a shared exon of two isoforms can’t be assigned to either with certainty alone. Salmon builds an “equivalence class” for every distinct set of compatible transcripts a read could have come from, then runs Expectation-Maximization: (E-step) given current abundance estimates, compute the probability each ambiguous read belongs to each compatible transcript; (M-step) re-estimate abundances from those probabilistic assignments; iterate to convergence. Reference: Patro R. et al. (2017). Nature Methods 14:417-419; Li B., Dewey C.N. (2011), BMC Bioinformatics 12:323.

Doublet simulation, precisely. Scrublet generates synthetic doublets by randomly averaging expression profiles of two randomly-chosen real cells, builds a combined PCA/kNN space containing both real cells and synthetic doublets, and computes each real cell’s local density of synthetic-doublet neighbors as a doublet score — real cells in a neighborhood dense with synthetic doublets are themselves likely true doublets. Reference: Wolock S.L., Lopez R., Klein A.M. (2019). Cell Systems 8(4):281-291.

PART 5 — Embedded Figures, Wet-Lab/Spatial Detail, and Further Worked Questions

1. Embedded reference figures

These are real, generated image files (not chat-only widgets), embedded so they render when this document is knitted to HTML/PDF.

Reference for the underlying relationship: Ewing B., Green P. (1998). “Base-calling of automated sequencer traces using phred. II. Error probabilities.” Genome Research 8:186-194.

Reference: DePristo M.A. et al. (2011). Nature Genetics 43:491-498 (defines Ti/Tv as a standard variant-calling QC metric).

Note on additional figures: the interactive pipeline-branch diagram and the TME single-cell workflow diagram referenced in Part 4 render as inline interactive SVG in the chat conversation itself, not as files — they cannot be embedded in this knitted document directly. If you want them as static image files inside this .Rmd, export them (e.g., screenshot or re-render as matplotlib/ggplot equivalents) and I can wire in the corresponding knitr::include_graphics() chunk pointing at that file.

2. Wet-lab context every computational person should know (so pipeline choices make sense)

Library prep determines what your FASTQ even represents. A computational error at Step 1 is often actually a wet-lab design choice showing up downstream: poly-A selection biases toward the 3’ end of transcripts (fine for standard gene-level quantification, useless for full-length isoform work), while ribo-depletion retains more RNA classes but at higher cost; total DNA input into library prep (measured typically via Qubit fluorometric quantification, with fragment-size distribution via a Bioanalyzer/TapeStation trace) directly determines PCR cycle number, which directly determines duplication rate seen in Step 2 QC. Reference: Illumina, “Ribo-Zero rRNA Removal Kit” and “TruSeq Stranded mRNA” technical notes, Illumina.com; Agilent TapeStation/Bioanalyzer application notes, agilent.com.

Exome capture vs. whole genome — a cost/coverage tradeoff every bioinformatician explains to wet-lab collaborators. WES uses hybridization capture probes (e.g., Agilent SureSelect, Illumina Exome) to pull down only the ~1-2% of the genome that is protein-coding exons, allowing much deeper sequencing (50-100x+) at lower cost per sample than WGS (typically 30x) — the tradeoff is that WES misses regulatory/intronic/intergenic variants entirely and has uneven capture efficiency across GC-rich/poor regions, which shows up as coverage dropout in Step 2 QC that WGS doesn’t have. Reference: Clark M.J. et al. (2011). “Performance comparison of exome DNA sequencing technologies.” Nature Biotechnology 29:908-914.

Single-cell dissociation bias — a known TME-specific artifact. Enzymatic tissue dissociation (needed to get single-cell suspensions from a solid tumor biopsy before loading onto the Chromium controller) is not equally gentle on every cell type — neurons and some epithelial subtypes are more fragile and systematically under-represented, and the dissociation process itself induces a stress-response transcriptional signature (immediate-early genes like FOS, JUN) in surviving cells that is a technical artifact, not real tumor biology, and needs to be recognized/regressed out or the sample processed via a cold-active protease protocol designed to minimize it. Reference: van den Brink S.C. et al. (2017). “Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations.” Nature Methods 14:935-936.

3. Spatial transcriptomics analysis pipeline (extending the scRNA-seq branch)

Because Part 3 introduced spatial transcriptomics as a TME-relevant technology, here is its actual computational pipeline, parallel in structure to Section 6 above:

Image + expression co-registration. Spatial platforms (Visium, Xenium) output both a histology image (H&E or immunofluorescence) and a spot/cell x gene expression matrix with x,y coordinates — the first computational step aligns these two coordinate systems so gene expression can be overlaid precisely on tissue morphology.
Spot deconvolution (Visium specifically). Visium spots (~55 micron diameter) typically capture multiple cells at once, so the raw spot-level expression is a mixture; deconvolution tools estimate the cell-type composition of each spot by referencing a matched scRNA-seq atlas of the same tissue as a signature basis. Reference: Cable D.M. et al. (2022). “Robust decomposition of cell type mixtures in spatial transcriptomics.” Nature Biotechnology 40:517-526 (RCTD).
Spatial domain/niche identification. Rather than clustering purely on expression (as in Section 6), spatial clustering also incorporates the x,y coordinates so that a “niche” cluster represents both similar expression and physical proximity — revealing structures like a tumor-immune boundary zone that pure expression clustering would miss. Reference: Dries R. et al. (2021). “Giotto: a toolbox for integrative analysis and visualization of spatial expression data.” Genome Biology 22:78.
Neighborhood enrichment analysis. Statistically tests whether specific cell-type pairs (e.g., exhausted T cells and TAMs) are found spatially closer together than expected by chance across the tissue — directly testing hypotheses about immune exclusion or immunosuppressive niches that non-spatial CellChat-style inference (Part 2C, Q12) can only suggest, not confirm spatially. Reference: Palla G. et al. (2022). “Squidpy: a scalable framework for spatial omics analysis.” Nature Methods 19:171-178.

4. Further worked interview questions

Q20. Why can’t you directly compare a UMAP plot’s axis distances as a meaningful quantity (e.g., “cluster A is twice as far from cluster B as from cluster C”)? UMAP explicitly optimizes to preserve local neighborhood structure (which points are near which), not global inter-cluster distances — the algorithm’s cross-entropy objective has no constraint forcing consistent global-distance semantics, and different random initializations can produce visually different relative cluster placements while preserving the same local neighbor relationships. This is a commonly-tested “gotcha” precisely because UMAP plots are so visually persuasive that people over-interpret them; the correct practice is to treat cluster identity/membership as the reliable output, and use quantitative tools (differential expression, marker scores) rather than plot geometry to make claims about how “different” two clusters are. Reference: McInnes L., Healy J., Melville J. (2018). arXiv:1802.03426 (discusses interpretation caveats); Chari T., Pachter L. (2023), “The specious art of single-cell genomics,” PLOS Computational Biology 19(8):e1011288 (a widely-cited critical discussion of dimensionality-reduction over-interpretation in the field).

Q21. Given a somatic VCF from Mutect2 with a variant at 8% variant allele frequency (VAF) in the tumor sample, what could explain this beyond “8% of tumor cells carry this mutation”? Several non-exclusive explanations a bioinformatician should consider: (1) tumor purity — if the biopsy is only 40% tumor cells (rest normal stroma/immune infiltrate), a mutation present in 20% of tumor cells alone would appear as ~8% VAF in bulk sequencing of the mixed sample; (2) subclonality — the mutation may be a genuinely late/subclonal event present in only a minority of tumor cells, informative about tumor evolution; (3) copy number state — if the locus has undergone a copy-number gain, VAF calculations must account for local ploidy, since simple VAF assumes a diploid background; (4) sequencing/mapping artifact at low VAF, particularly in low-complexity or repetitive regions, which Mutect2’s built-in filters (and orientation-bias filtering for FFPE artifacts) are specifically designed to catch. Distinguishing these requires tumor purity/ploidy estimation tools (e.g., ABSOLUTE, FACETS) applied jointly with the VAF. Reference: Carter S.L. et al. (2012). “Absolute quantification of somatic DNA alterations in human cancer.” Nature Biotechnology 30:413-421 (ABSOLUTE); Shen R., Seshan V.E. (2016), “FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing,” Nucleic Acids Research 44(16):e131.

Q22. You’re asked to estimate tumor mutational burden (TMB) from an exome — walk through the calculation and common pitfalls. TMB is typically reported as (number of nonsynonymous somatic mutations passing filters) / (total exonic territory sequenced, in megabases), giving mutations/Mb — used clinically as a biomarker for likely response to checkpoint-blockade immunotherapy, since a higher mutation burden generally correlates with more neoantigens for T cells to recognize. Pitfalls: the denominator (captured exonic Mb) differs between capture kits, making raw TMB values non-comparable across different sequencing panels without harmonization; germline variants must be correctly subtracted (requiring a matched normal sample or population-database filtering if unmatched); and different labs’ filtering thresholds (VAF cutoffs, depth minimums) materially shift the count, which is why standardization efforts (Friends of Cancer Research TMB Harmonization Project) exist. Reference: Chalmers Z.R. et al. (2017). “Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden.” Genome Medicine 9:34; Merino D.M. et al. (2020), “Establishing guidelines to harmonize tumor mutational burden (TMB): in silico assessment of variation in TMB quantification across diagnostic platforms,” Journal for ImmunoTherapy of Cancer 8:e000147.

Q23. What’s the difference between “sensitivity” and “precision” for a variant caller, and why does a caller’s default configuration usually favor one over the other? Sensitivity (recall) = true positives / (true positives + false negatives) — the fraction of real variants the caller successfully found. Precision = true positives / (true positives + false positives) — the fraction of called variants that are actually real. Germline diagnostic pipelines typically tune toward higher sensitivity (accepting some false positives that get filtered downstream via annotation/frequency/inheritance pattern in Part 1 Section 4) because missing a real pathogenic variant in a patient is a more serious clinical failure than flagging an extra false lead for a clinician to review and dismiss; conversely, somatic calling in a high-throughput cancer screening context may tune toward precision to avoid overwhelming downstream review with sequencing artifacts. Benchmarking against a truth set (e.g., Genome in a Bottle’s curated reference genomes) with tools like hap.py is the standard way to actually measure both quantities for a given caller/parameter configuration rather than assuming. Reference: Zook J.M. et al. (2019). “An open resource for accurately benchmarking small variant and reference calls.” Nature Biotechnology 37:561-566 (Genome in a Bottle); Krusche P. et al. (2019), “Best practices for benchmarking germline small-variant calls in human genomes,” Nature Biotechnology 37:555-560 (hap.py).

NGS & Computational Genomics Pipelines: FASTQ to VCF and Beyond

A detailed, referenced field guide — WGS/WES, Bulk RNA-seq, scRNA-seq, Tumor Microenvironment, and a buildable Nextflow DSL2 pipeline

Compiled reference document

2026-08-13