🖥️
Command line & samtools
The shell basics plus the genomics toolkit you use daily.
Navigating & files
pwd | where am I (print working dir) |
ls -lh / ls -a | list (sizes) / include hidden |
cd dir / cd .. / cd ~ | enter / up one / home |
mkdir -p a/b | make folders + parents |
cp -r src dst | copy (-r for folders) |
mv src dst | move / rename |
rm file / rm -r dir | delete (no undo!) |
chmod +x script.sh | make a script executable |
find . -name "*.bam" | find files by pattern |
Text & data wrangling
cat / less file | print / page (q quits) |
head -n 20 / tail -n 5 | first / last lines |
wc -l file | count lines |
grep "pat" file | find matching lines (-v invert, -c count) |
cut -f1,3 | pick columns (tab-delim) |
sort | uniq -c | sort then count duplicates |
awk '{print $1}' | print column 1 |
sed 's/a/b/g' | find & replace |
Pipes, redirects & survival
cmd | cmd2 | pipe output into next command |
cmd > out / >> out | write / append to file |
cmd && cmd2 | run cmd2 only if cmd succeeds |
*.fastq wildcard | match all matching files |
| Tab / Ctrl-C / ↑ | autocomplete / stop / last cmd |
gzip / tar -xzvf | compress / extract archive |
samtools (BAM/SAM)
samtools view -b in.sam > out.bam | SAM → BAM (-b) |
samtools sort in.bam -o s.bam | sort by position (do before index) |
samtools index s.bam | build the .bai index |
samtools flagstat in.bam | mapping summary (% mapped…) |
samtools stats in.bam | detailed QC statistics |
samtools idxstats in.bam | reads per chromosome |
samtools view in.bam chr1:1-1000 | reads in a region (needs index) |
samtools view -f 4 / -F 4 | keep / exclude by FLAG (4 = unmapped) |
samtools view -q 30 | filter by mapping quality |
samtools markdup | mark PCR duplicates |
samtools depth / coverage | per-base depth / per-region coverage |
samtools merge out.bam *.bam | merge BAMs |
samtools faidx ref.fa | index a FASTA (.fai) |
bcftools (variants)
bcftools mpileup -f ref.fa in.bam | pile up reads vs reference |
... | bcftools call -mv -Oz -o out.vcf.gz | call variants (-m multiallelic, -v variant sites, -Oz bgzip) |
bcftools view -f PASS in.vcf.gz | view / filter records |
bcftools norm -f ref.fa | normalize & left-align indels |
bcftools filter -e 'QUAL<20' | exclude by expression |
bcftools index in.vcf.gz | index a bgzipped VCF |
bcftools stats in.vcf.gz | summary (ts/tv, counts…) |
⚠️ Order matters: sort → index before any region query. You cannot index an unsorted BAM, and most tools need the matching
.bai/.fai.Full lessons: Bash basics · The samtools survival kit · Variant calling
🧬
File formats
What each bioinformatics file holds, and the traps.
FASTA — sequences
>seq_id optional description ACGTACGTACGTACGT...
Header line starts with >; sequence may wrap over many lines. DNA, RNA, or protein. Index with samtools faidx.
FASTQ — reads + quality
@read_id # 1 name (@) ACGTACGTACGT # 2 sequence + # 3 separator (+) IIIIFFFF;;;; # 4 quality (1 char/base)
| Encoding | Phred+33 (Sanger/Illumina 1.8+); ASCII − 33 = Q |
| Q20 | 1 error in 100 (99%) |
| Q30 | 1 in 1,000 (99.9%) — common target |
| Q40 | 1 in 10,000 (99.99%) |
GTF / GFF — annotation
9 tab-separated columns (seqname, source, feature, start, end, score, strand, frame, attributes). Maps coordinates to genes/exons/CDS. 1-based.
SAM / BAM — aligned reads
SAM = text, BAM = compressed binary. 11 mandatory columns:
| QNAME | read name |
| FLAG | bitwise (4 = unmapped, 16 = reverse, 2 = proper pair, 1024 = duplicate) |
| RNAME / POS | reference & 1-based position |
| MAPQ | mapping quality (Phred) |
| CIGAR | M match, I insert, D delete, S soft-clip, N skip |
| SEQ / QUAL | read bases & qualities |
VCF — variants
| CHROM / POS | location (1-based) |
| REF / ALT | reference vs alternate allele |
| QUAL / FILTER | confidence / PASS or fail |
| INFO | site annotations (DP, AF…) |
| FORMAT + sample | genotype: 0/0 hom-ref, 0/1 het, 1/1 hom-alt |
BED — regions
chr1 200 1000 name # start is 0-based
| .fai / .bai | FASTA / BAM index files |
| .vcf.gz + .tbi | bgzipped VCF + tabix index |
⚠️ Coordinate trap: BED is 0-based (half-open); GTF, SAM, VCF are 1-based. Mixing them shifts everything by one.
Full lessons: File-format alphabet soup · Sequences, FASTA & FASTQ
📦
conda & Git
Environments that never break, and version control you can trust.
conda — one-time setup
conda config --add channels conda-forge conda config --add channels bioconda conda config --set channel_priority strict
conda — daily use
conda create -n env python=3.11 | new environment |
conda activate env | turn it on |
conda install -c bioconda samtools | install a tool |
conda list / env list | packages / all envs |
conda search bwa | find a package |
conda deactivate | turn off |
conda env remove -n env | delete an environment |
conda env export --from-history > environment.yml | save a portable spec |
conda env create -f environment.yml | recreate anywhere |
mamba install ... | faster drop-in for conda |
One env per project. Never install into
base. Pin versions that affect results.Git — the daily loop
git clone <url> | copy a repo locally |
git status | what changed |
git add file / git add . | stage changes |
git commit -m "msg" | save a snapshot |
git push / git pull | send to / get from GitHub |
git log --oneline | history |
git diff | see unstaged changes |
Git — branches & fixes
git checkout -b feature | create + switch branch |
git merge feature | merge into current branch |
git stash / stash pop | shelve / restore changes |
git restore file | discard local edits to a file |
git push -u origin main | first push of a branch |
# .gitignore — never commit these *.bam *.fastq.gz data/ .env # secrets / API keys __pycache__/
The loop: edit → add → commit → push. Commit small and often, with clear messages.
Full lessons: Conda environments · Git & GitHub for scientists
📊
Stats & DESeq2
The numbers you must interpret, and the workflow.
Statistics quick reference
| p-value | P(data this extreme | null true). Not the chance the result is real. |
| padj / FDR | p corrected for many tests (Benjamini-Hochberg). Filter on this, not raw p. |
| Bonferroni | strict: significant if p < 0.05 / (number of tests) |
| effect size | how big the change is (log2FC, Cohen's d). Report with the p-value. |
| log2FC | 1 = 2×, 2 = 4×, 3 = 8×; negative = down |
| multiple testing | 20,000 genes at p<0.05 ≈ 1,000 false hits by chance |
| neg. binomial | the count model RNA-seq uses (over-dispersed); don't t-test raw counts |
A microscopic p-value with a near-zero effect size is significant but biologically boring. Use both axes (that's what a volcano plot shows).
Volcano thresholds (typical)
| padj < 0.05 | significance cutoff (y-axis) |
| |log2FC| > 1 | at least 2-fold change (x-axis) |
DESeq2 workflow (R)
# 1. load library(DESeq2) # 2. build: counts + sample table + design dds <- DESeqDataSetFromMatrix( countData = counts, colData = coldata, design = ~ condition) # 3. (optional) pre-filter low counts dds <- dds[rowSums(counts(dds)) >= 10, ] # 4. run the analysis dds <- DESeq(dds) # 5. results: treated vs control res <- results(dds, contrast = c("condition","treated","control")) # 6. shrink LFCs for ranking & plots res <- lfcShrink(dds, coef=2, type="apeglm") # 7. significant genes sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1) # 8. handy extras plotMA(res) # MA plot vsd <- vst(dds) # for PCA / heatmaps counts(dds, normalized=TRUE) # normalized counts
Same idea in edgeR:
DGEList → calcNormFactors → estimateDisp → glmQLFit → glmQLFTest.Full lessons: Stats you can't avoid · The multiple-testing problem · Differential expression with DESeq2