RNA-seq: Measuring Gene Expression

RNA-seq: Measuring Gene Expression

Every cell in your body carries the same genome, yet a hepatocyte, a neuron, and a pancreatic beta cell look and behave nothing alike. The difference lies not in which genes are present, but in which genes are actively transcribed. This set of actively transcribed RNAs—the transcriptome —changes with tissue type, developmental stage, disease state, and drug exposure. For a pharmacologist, measuring the transcriptome means asking directly: which genes does this compound upregulate, which does it silence, and what pathway does it perturb? RNA sequencing (RNA-seq) is the dominant technology for answering these questions quantitatively and genome-wide. The wet-lab portion of an RNA-seq experiment is conceptually simple. RNA is extracted from a biological sample, purified (typically enriching for messenger RNA via poly-A selection, or depleting ribosomal RNA), fragmented, and reverse-transcribed into cDNA. Sequencing adapters are ligated, and the library is sequenced on a high-throughput platform, most commonly Illumina. The output is a set of FASTQ files : millions of short reads (usually 50–150 bases), each carrying a sequence and a per-base quality score. A typical bulk RNA-seq library contains 20–40 million read pairs per sample, and a well-designed experiment includes at least three biological replicates per condition—replication, not read depth, is what powers statistical inference. The computational challenge is to convert millions of reads into a table of gene expression values. The standard pipeline proceeds in four steps: (1) quality control of raw reads, (2) trimming of adapters and low-quality bases, (3) alignment of reads to a reference genome using a splice-aware aligner, and (4) quantification , in which reads overlapping annotated genes are tallied. Splice-awareness matters because mature mRNA lacks introns: an aligner such as HISAT2 or STAR must allow reads to split across exon–exon junctions rather than discard them. The result of step 4 is the count matrix —the central data object of all downstream analysis—with genes as rows, samples as columns, and integer counts as entries. Below is a minimal command-line workflow for one sample. In practice you would wrap this in a loop or workflow manager to process all samples identically. # Step 1: Assess read quality per sample fastqc sample_R1.fastq.gz sample_R2.fastq.gz # Step 2: Trim adapters and low-quality bases (min length 25) cutadapt -a AGATCGGAAGAGC -q 20 -m 25 \ -o trimmed_R1.fastq.gz -p trimmed_R2.fastq.gz \ sample_R1.fastq.gz sample_R2.fastq.gz # Step 3: Splice-aware alignment to the human genome, then sort and index hisat2 -p 8 -x GRCh38_index \ -1 trimmed_R1.fastq.gz -2 trimmed_R2.fastq.gz \ | samtools sort -o sample_aligned.bam samtools index sample_aligned.bam # Step 4: Count reads overlapping exons, summarized per gene featureCounts -p -t exon -g gene_id \ -a genes.gtf -o sample_gene_counts.txt sample_aligned.bam Each tool plays a defined role. fastqc produces per-sample quality reports that you inspect before committing to downstream analysis; systematic adapter contamination or a drop in quality at read ends signals that trimming parameters need adjustment. cutadapt removes technical artifacts so they do not inflate mismatch rates during alignment. hisat2 outputs a SAM stream that samtools converts to a sorted, indexed BAM file—the standard format for aligned reads. Finally, featureCounts intersects each aligned read with exon annotations from a GTF file and assigns it to the overlapping gene. Reads that map equally well to multiple genes