# BWA Alignment Tutorial: Step-by-Step Guide By the end of this bwa alignment tutorial you will have aligned paired-end Illumina reads to a reference genome, produced a sorted and indexed BAM file, measured alignment quality with samtools, and called variants from the result. Every command is real and runnable. The practice dataset is generated on your own machine, so the workflow is fully reproducible. You need: - A Linux or macOS shell (Windows works through WSL2) - About 4 GB of free disk space for the practice dataset - Either a conda installation or a C compiler with zlib headers - A laptop is enough for the bacterial genome used here The reason for a careful, stepwise approach is practical. Modern sequencers produce enormous amounts of data at low cost, and the bottleneck for most new groups is knowing how to process what comes off the instrument [1]. Alignment is the first computational step in nearly every analysis that follows, from mapping-by-sequencing experiments that locate causal mutations to routine clinical variant calling. Do it wrong and everything downstream inherits the error. ## How BWA Works BWA (Burrows-Wheeler Aligner) maps short DNA reads against a large reference genome. Heng Li and Richard Durbin published it in 2009, and it remains the default short-read aligner in most variant-calling pipelines. A naive aligner would compare every read against every position of the genome, which is infeasible at gigabase scale. BWA avoids this with the Burrows-Wheeler transform (BWT), a reversible rearrangement of the reference sequence that supports fast string search with little memory. BWA builds an index from the BWT, uses it to find exact-match seeds for each read, then extends and scores candidate positions with a local alignment algorithm (Smith-Waterman). This design lets a laptop align reads against a 3-gigabase mammalian genome within a few gigabytes of RAM. BWA ships three algorithms: - **bwa aln**, the original backtrack algorithm for reads up to about 100 bp. It runs in two stages (aln, then samse or sampe) and is rarely used today. - **bwasw**, a local aligner for long reads, superseded by minimap2. - **bwa mem**, the current default for reads from roughly 70 bp into the megabase range. It is fast, multithreaded, and detects split alignments. BWA writes output in SAM format (Sequence Alignment/Map), a tab-separated text standard. SAM compresses to BAM (the binary form), and BAM indexes for random access. The samtools suite handles sorting, indexing, and inspection. The rest of this guide uses bwa mem with samtools and bcftools. ## BWA Alignment Tutorial: The Workflow at a Glance 1. Install BWA and its companion tools. 2. Download a reference genome in FASTA format. 3. Index the reference with bwa index. 4. Obtain reads in [FASTQ format](/knowledge/bioinformatics/fastq-file-format-phred-scores-qc) (this tutorial simulates them with wgsim). 5. Align with bwa mem and sort to BAM in one pipe. 6. Inspect the alignment with samtools flagstat, idxstats, and stats. 7. Judge the results with mapping rates, pairing rates, and coverage. 8. Mark duplicates and call variants with samtools and bcftools. ## Step 1: Install BWA and Its Companion Tools The fastest route is conda with the bioconda channel, the standard distribution point for bioinformatics software: ``` conda create -n align -c conda-forge -c bioconda bwa samtools bcftools conda activate align bwa 2>&1 | head -n 3 ``` Expected output: ``` Program: bwa (alignment via Burrows-Wheeler transformation) Version: 0.7.17-r1188 Contact: Heng Li ``` Any 0.7.17 or newer release is fine. If you prefer to build from source: ``` git clone https://github.com/lh3/bwa.git cd bwa make ``` This needs a C compiler and the zlib development headers (the zlib1g-dev package on Debian or Ubuntu). The binary appears as ./bwa in that directory, and you can copy it somewhere on your PATH. Install wgsim the same way, because it generates the practice reads: ``` git clone https://github.com/lh3/wgsim.git cd wgsim gcc -O2 -o wgsim wgsim.c -lz -lm ``` Verify everything is reachable: ``` which bwa samtools bcftools ``` Each command should print an install path. If bwa prints nothing, the conda environment is not active or the binary is not on your PATH. ## Step 2: Get a Reference Genome The aligner needs the reference in a single FASTA file, the plain-text standard for [nucleotide sequence](/knowledge/molecular-biology/nucleotide-sequence). Reference sequences are called contigs (named sequences, whether whole chromosomes or fragments). Three rules matter: - Sequence names must be unique. Duplicates break indexing. - Pick one genome build per project and never mix builds. The human builds hg19 and GRCh38 use different coordinate systems, so a position in one does not describe the same base in the other. - Keep the file stable. Renaming it after indexing breaks the link to the index (Step 3). For this tutorial, download the reference genome of *Escherichia coli* K-12 MG1655 from NCBI RefSeq: ``` wget https://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Escherichia_coli/reference/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz mv GCF_000005845.2_ASM584v2_genomic.fna ecoli.fna head -n 1 ecoli.fna ``` Expected header line: ``` >NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome ``` This genome is 4,641,652 bp, small enough that every step runs in seconds, which makes it ideal for learning. For human, mouse, or plant work, download the primary assembly FASTA from a major archive. Ensembl and the [UCSC genome browser](/knowledge/bioinformatics/ucsc-genome-browser-beginner-tutorial-and-tips) have served as the standard retrieval points for genomic sequence in [comparative genomics](/blog/guides/comparative-genomics) for two decades [2], and NCBI RefSeq distributes the same underlying assemblies. Record the URL, the assembly version, and the download date in a project README. Your future self and your reviewers will need them. ## Step 3: Index the Reference BWA cannot search raw FASTA efficiently. It needs an index built once per reference: ``` bwa index ecoli.fna samtools faidx ecoli.fna ls -lh ecoli.fna.* ``` Expected output from bwa index (condensed): ``` [bwa_index] Pack FASTA... 0.08 sec [bwa_index] Construct BWT for the packed sequence... [BWTIncCreate] text length=4641652, bucket size=4096 [bwt_gen] Sorted 4641652 characters in 0.4 sec(s) [bwa_index] Update BWT... 0.02 sec [bwa_index] Construct SA from BWT and Occ... 0.30 sec [main] Version: 0.7.17-r1188 [main] CMD: bwa index ecoli.fna ``` The command produces five files, and bwa mem locates them by the reference prefix: | File | Contents | |--|--| | ecoli.fna.amb | Positions of ambiguous bases | | ecoli.fna.ann | Sequence names and lengths | | ecoli.fna.bwt | The Burrows-Wheeler transform | | ecoli.fna.pac | The packed reference sequence | | ecoli.fna.sa | The suffix array | samtools faidx creates ecoli.fna.fai, a per-contig offset table that lets samtools and variant callers jump to regions of the FASTA. Scale expectations: the E. coli index builds in seconds. A [human genome](/blog/guides/human-genome) takes on the order of an hour, several gigabytes of RAM, and a few gigabytes of disk. Build it once on the machine with the most memory, then share the reference directory with your group. Do not rename the FASTA afterward, because the index files carry the original prefix in their names. ## Step 4: Prepare Your Reads ### FASTQ, the language of reads Sequencers deliver reads as FASTQ. Each read occupies four lines: ``` @A00345:123:HCY7W:1:1101:242:1000 1:N:0:ATCACG GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG + FFFFFFFF:FFFFFF,FFFFFF:F,FFF:FFFFFFF,FFFFFFF,FFFFFF,F:FFFFFFFFF, ``` Line 1 is the read name, which marks the first or second read of a pair. Line 2 is the sequence. Line 3 is a separator. Line 4 is the quality string, one character per base, encoded as Phred+33: the ASCII value of the character minus 33 gives the Phred score Q, where Q20 means a 1 percent chance the base is wrong and Q30 means 0.1 percent. Modern Illumina data uses Phred+33. Legacy data encoded as Phred+64 will score nonsense, so convert it first (fastp and seqtk both handle this). ### Generate practice reads with wgsim Rather than depend on an external download, simulate reads with known truth: ``` ./wgsim/wgsim -N 500000 -1 150 -2 150 ecoli.fna r1.fq r2.fq gzip r1.fq r2.fq ``` This writes 500,000 paired-end reads of 150 bp each. wgsim introduces sequencing errors and de novo mutations, tunable with the -e and -r options, which gives variant calling something real to find. Because the simulator knows where every read came from, you can check your pipeline against ground truth. Simulation is a standard technique for benchmarking genomics tools when labeled real data is scarce, used everywhere from phylogenomic pipeline testing [3] to the evaluation of single-cell copy-number callers [4]. ### Work out the coverage you are creating ![Worked coverage sum: 500,000 read pairs of 150 bp over a 4,641,652 bp genome gives 32.3x depth.](/article-images/illustrations/bwa-alignment-tutorial-coverage-calc.webp) Coverage (depth) is the average number of reads covering each base: ``` depth = (number of reads x read length) / genome size depth = (500000 pairs x 2 x 150 bp) / 4641652 bp = 32.3x ``` Thirty-fold coverage is a common target for small-genome variant calling. To hit 30x on this genome you would need about 464,000 pairs of 150 bp reads. Real libraries rarely distribute coverage this evenly, which is why you measure it again after alignment. ## Step 5: Align with bwa mem The alignment command, with a read group and sorting in the same pipe: ``` bwa mem -t 4 -R '@RG\tID:ecoli_test\tSM:ecoli_test\tPL:ILLUMINA\tLB:wgsim_sim' ecoli.fna r1.fq.gz r2.fq.gz | samtools sort -@ 4 -o ecoli.sorted.bam samtools index ecoli.sorted.bam ``` Breaking it down: - bwa mem reads the two FASTQ files as a pair. They must be in the same order, mate for mate. BWA reads gzipped FASTQ directly. - -t 4 uses four threads. Set this to the cores you can spare. - -R adds a read group, a metadata line attached to every alignment record. ID identifies the run, SM is the sample name (variant callers group files by SM), PL is the platform, and LB is the library. Set these on the first run. GATK and several bcftools workflows refuse BAM files without a read group, and adding one after the fact is tedious. - The pipe into samtools sort converts the SAM text stream into a coordinate-sorted BAM without writing a huge intermediate SAM file. If you want the raw SAM, redirect to a file instead. - samtools index builds ecoli.sorted.bam.bai, the random-access index that region queries and genome browsers require. BWA writes progress to stderr and alignments to stdout, which is what makes the pipe work. Expect output like this: ``` [M::bwa_load_idx] The BWT index was loaded in 0.25 sec. [M::process] read 80000 sequences (24000000 bp)... [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (0, 499999, 0, 0) [M::mem_pestat] skip orientation FF as there are not enough pairs [M::mem_pestat] skip orientation RF as there are not enough pairs [M::mem_pestat] skip orientation RR as there are not enough pairs [M::mem_process_seqs] Processed 500000 reads in 14.0 sec. [main] Version: 0.7.17-r1188 [main] CMD: bwa mem -t 4 -R ... ecoli.fna r1.fq.gz r2.fq.gz ``` The mem_pestat lines show BWA estimating insert size and orientation from the data itself. Here nearly all pairs vote for the FR (forward-reverse) orientation, typical for standard Illumina libraries. The final CMD line matters more than it looks: BWA records the exact command into the BAM header as a @PG line, giving you free provenance. If older Picard or GATK versions mishandle split reads in your pipeline, add -M, which marks the shorter hit of a split alignment as secondary. Modern tools do not need it. ## Step 6: Sort, Index, and Inspect the BAM ### Verify the file ``` samtools quickcheck -v ecoli.sorted.bam ``` No output means the BAM is structurally sound. Then read the header: ``` samtools view -H ecoli.sorted.bam | head -n 4 ``` Expected (tabs shown as spaces for readability): ``` @HD VN:1.6 SO:coordinate @SQ SN:NC_000913.3 LN:4641652 @RG ID:ecoli_test SM:ecoli_test PL:ILLUMINA LB:wgsim_sim @PG ID:bwa PN:bwa VN:0.7.17-r1188 CL:bwa mem -t 4 ... ``` @SQ lines list the reference contigs and lengths. If they do not match your FASTA, something is wrong. ### The one command to run every time ``` samtools flagstat ecoli.sorted.bam ``` Expected output for the simulated data (your numbers will differ slightly): ``` 1000000 + 0 in total (QC-passed reads + QC-failed reads) 1000000 + 0 primary 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 0 + 0 primary duplicates 999200 + 0 mapped (99.92% : N/A) 1000000 + 0 paired in sequencing 500000 + 0 read1 500000 + 0 read2 998400 + 0 properly paired (99.84% : N/A) 999000 + 0 with itself and mate mapped 200 + 0 singletons (0.02% : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5) ``` How to read it: - **mapped**, the fraction of reads placed anywhere on the reference. This is the headline number. - **properly paired**, reads whose mate mapped nearby in the expected orientation. Low values point to wrong insert sizes, chimeric pairs, or contaminated libraries. - **singletons**, reads that mapped while their mate did not. - **secondary** and **supplementary**, alternative and split alignments. They are not independent observations, so never count them as extra evidence. ### Per-contig counts and depth ``` samtools idxstats ecoli.sorted.bam ``` ``` NC_000913.3 4641652 999200 200 * 0 0 600 ``` Columns are sequence name, sequence length, mapped reads, and unmapped reads placed by their mate, with fully unplaced reads aggregated in the asterisk row. Mean depth: ``` samtools depth -a ecoli.sorted.bam | awk '{s+=$3} END {print s/NR}' ``` This prints a number close to 32, matching the coverage you designed in Step 4. The -a flag includes zero-coverage positions so the denominator is the true genome size. For a fuller picture, samtools stats reports error rate, mapping quality distribution, and the insert size histogram (the IS section): ``` samtools stats ecoli.sorted.bam > ecoli.stats grep '^IS' ecoli.stats | head -n 10 ``` The insert size peak should sit near your library preparation target. A peak far from expectation suggests sizing problems or contamination. ## Step 7: Judge Whether the Alignment Worked Absolute numbers depend on the organism, the library, and the reference. The patterns matter more than any single threshold. - **Mapping rate.** Clean whole-genome Illumina data against a matching reference commonly maps above 95 percent. Substantially lower rates usually mean the wrong build, the wrong organism, heavy adapter contamination, or a divergent sample. A perfect 100 percent on real data is itself suspicious and often indicates reference contamination with the sample. - **Properly paired rate.** Standard short-insert libraries pair properly at a high rate. A collapse here with a normal mapping rate suggests insert size problems or read-through into adapters. - **Mapping quality.** MAPQ is a Phred-scaled estimate of how likely an alignment is wrong. Q30 corresponds to a 1 in 1,000 chance. Filter downstream steps by MAPQ, for example with samtools view -q 20. - **Coverage uniformity.** Real libraries dip and spike. Large zero-coverage blocks mean the sample is not what you think, or the reference contains regions absent from the sample. Two habits keep you honest. First, inspect a region you know: ``` samtools view ecoli.sorted.bam NC_000913.3:1000000-1010000 | head -n 5 ``` Each line is one alignment record with the read name, FLAG, reference name, position, MAPQ, and the CIGAR string (the M, I, and D operations describing how the read aligns). Second, remember that alignment choices propagate. Decisions made during mapping, from reference choice to parameter settings, demonstrably influence downstream results such as variant calls and phasing [3], so a mapping rate you did not investigate is a result you do not understand. ## Step 8: Mark Duplicates and Call Variants PCR and optical duplicates stack reads at the same coordinates and inflate the pileup (the column of reads a variant caller examines at each position). BWA does not mark them, so do it after alignment. The samtools-native route: ``` samtools collate -O ecoli.sorted.bam | samtools fixmate -m - - | samtools sort -o ecoli.fixmate.bam samtools markdup ecoli.fixmate.bam ecoli.markdup.bam samtools flagstat ecoli.markdup.bam | head -n 6 ``` fixmate -m adds the tags markdup needs, and markdup flags duplicates so downstream tools can ignore them. Picard MarkDuplicates does the same job and remains common in published pipelines. With duplicates flagged, call variants: ``` bcftools mpileup -f ecoli.fna -q 20 -Q 20 ecoli.markdup.bam | bcftools call -mv -Oz -o ecoli.vcf.gz bcftools index ecoli.vcf.gz ``` The output is a VCF (Variant Call Format) file listing every position where the sample differs from the reference. With wgsim's default mutation rate of 0.001, you should recover on the order of 4,600 SNPs across the 4.64 Mb genome, a built-in sanity check for the whole pipeline. For long-term storage, CRAM saves roughly half the space of BAM by referencing the FASTA: ``` samtools view -C -T ecoli.fna -o ecoli.sorted.cram ecoli.sorted.bam ``` CRAM requires the reference to be available whenever the file is read, so keep both together. ## bwa mem Options That Matter | Option | Default | What it does | When to change it | |--|--|--|--| | -t | 1 | Number of threads | Set to the cores you can spare | | -k | 19 | Minimum seed length | Lower to 11 to 15 for short or degraded reads | | -w | 100 | Band width for gap filling | Raise for reads spanning large indels | | -T | 30 | Minimum score to output an alignment | Lower for divergent species or damaged DNA | | -M | off | Mark shorter split hits as secondary | Add only for older Picard/GATK compatibility | | -R | none | Read group line | Always set ID, SM, PL, and LB | | -x | none | Read type preset | Use pacbio or ont2d for those platforms | | -I | estimated | Insert size override | Rarely needed, BWA measures it from the data | Two parameters deserve thought. Lowering -k lets BWA seed on shorter exact matches, which helps divergent or fragmented reads but increases spurious hits and runtime. Lowering -T keeps weaker alignments, useful when the sample is distant from the reference, at the cost of more noise. Change one parameter at a time and compare flagstat output between runs. ## Practical Tips for Reliable Alignments - Keep the reference and its index files together in one directory, and never rename the FASTA after indexing. - Use one conda environment per project with pinned versions, and record the versions in the README. - Trust the @PG header for provenance instead of shell history. It stores the exact command that produced each BAM. - For large batches, bwa-mem2 produces essentially identical alignments about twice as fast at the cost of more RAM. The commands are bwa-mem2 index followed by bwa-mem2 mem. - Test parameter changes on a slice first. For gzipped FASTQ, `zcat r1.fq.gz | head -n 400000 | gzip > r1.small.fq.gz` yields 100,000 reads. - Watch stderr during long runs. The Processed lines report progress in reads, so a stalled job is obvious. - Keep read group conventions consistent. Use the same SM for technical replicates of one sample and a different ID for each run. - For human work against GRCh38 with its ALT contigs, follow the alt-aware pipeline described in the BWA documentation, because naive alignment to ALT contigs can distort variant calls. ## Common Mistakes and Limitations 1. **bwa: command not found.** The conda environment is not active, or the compiled binary is not on your PATH. Run conda activate align, or add the build directory to PATH. 2. **[E::bwa_idx_load_from_disk] fail to locate the index files** (older versions phrase it slightly differently). BWA cannot find the index. Either you never ran bwa index, or you renamed or moved the FASTA afterward, because BWA locates the index by the reference prefix. Re-run bwa index on the exact path you pass to bwa mem. 3. **samtools index fails with an error containing "Unsorted positions".** The BAM is not coordinate-sorted. Run samtools sort first, then index the sorted file. 4. **GATK or bcftools complains about a missing read group.** The BAM lacks an @RG line. Re-align with -R, or add one with samtools addreplacerg. Going forward, set -R at alignment time. 5. **Zero or very low mapping rate.** Check, in order: the reference build matches the sample, contig naming matches (chr1 versus 1), the FASTQ files are the organism you think they are, and quality encoding is Phred+33. Mixing hg19 and GRCh38 within a project is the classic version of this failure, because coordinates differ between builds by up to hundreds of kilobases in places. 6. **Desynchronized paired FASTQ files.** If you filtered or trimmed one file independently of its mate, the pairs no longer line up and BWA pairs the wrong reads together, often silently. Re-sync the files before aligning, or filter both files together. 7. **Out of memory during indexing.** The Linux OOM killer prints Killed and the process dies. A human genome index needs several gigabytes of RAM. Build it on the largest machine available and share the directory. 8. **Duplicate sequence names in the reference.** samtools faidx reports [E::fai_insert_index] Ignoring duplicate sequence, and downstream behavior becomes unpredictable. Make every FASTA name unique before indexing. 9. **Counting secondary and supplementary alignments as evidence.** flagstat lists them separately for a reason. Coverage and variant tools should see primary alignments, filtered by MAPQ. 10. **Expecting BWA to do everything.** BWA aligns. It does not trim adapters, mark duplicates, remove contaminants, call variants, or assess quality. Those are separate steps with separate tools, and skipping them is how bad variant calls happen. The broader limitation is that every result depends on the reference. Aligning to a poor or distant reference produces low rates and biased calls no matter how well you run the commands, and each project needs its own reference choice and validation. ## When BWA Is Not the Right Tool ![Comparison of bwa mem, STAR or HISAT2, and minimap2 by read type and when each is the right tool.](/article-images/illustrations/bwa-alignment-tutorial-aligner-choice.webp) Three situations call for different software. **RNA-seq.** Reads cross exon-exon junctions, so a contiguous aligner like BWA will drop or misplace them. Use a spliced aligner such as STAR or HISAT2. **Very long reads.** For Oxford Nanopore or PacBio data, minimap2 is the standard aligner. bwa mem can process such reads with its presets, but minimap2 is faster and better tested at those lengths. Mixed short- and long-read strategies are now routine in genomics, for example recent [mitochondrial genome](/blog/guides/mitochondrial-genome) assemblies built from both second- and third-generation data [5], and the long-read half of such workflows belongs to long-read tools. **Metagenomics.** When reads come from many species, aligning to a single reference answers the wrong question. Use tools built for taxonomic classification. **Non-model organisms.** If no good reference exists, alignment rates will suffer against a distant relative. Consider assembling a reference from your own data first, then aligning to it. ## Where to Go Next You now have the core loop of read alignment. If your goal is to learn bioinformatics, extend the same project rather than starting a new one. Add FastQC before alignment and MultiQC after, pull real data for your organism from SRA or ENA and rerun the pipeline unchanged, visualize the BAM in IGV, and build on the bcftools variant calling you already ran. A structured [bioinformatics course](/blog/careers/bioinformatics-courses-and-certifications-a-practical-guide) will cover the same ground with more theory, but nothing substitutes for running a pipeline end to end on data you generated yourself, then breaking it on purpose to see what the failure looks like. ## Frequently Asked Questions ### What is BWA used for? BWA aligns short [DNA sequencing](/blog/guides/dna-sequencing) reads against a reference genome. It is the standard first step for resequencing projects, variant calling, and mapping experiments on Illumina-style data. ### Which BWA command should I use, aln, sw, or mem? Use bwa mem for essentially all modern work. bwa aln is a legacy algorithm for reads of about 100 bp or shorter, and bwasw has been superseded by minimap2 for long reads. ### How much memory and time does bwa index need? Bacterial genomes index in seconds on a laptop. A human genome takes on the order of an hour, several gigabytes of RAM, and a few gigabytes of disk, once per reference. ### What mapping rate should I expect? Clean whole-genome Illumina data against a matching reference commonly maps above 95 percent. Lower rates point to the wrong build, contamination, adapter problems, or a divergent sample. ### What is a read group and why does it matter? A read group is a metadata line with ID, SM, PL, and LB fields attached to each alignment. Variant callers use the sample name to group files, and many tools refuse BAM files without one. ### Should I use BWA-MEM2 instead of BWA? Yes for large batches. BWA-MEM2 produces essentially identical alignments roughly twice as fast while using more memory, with the commands bwa-mem2 index and bwa-mem2 mem. ### Can BWA align Nanopore or PacBio reads? bwa mem can handle long reads with its read type presets, but minimap2 is the standard and better-tested aligner for Nanopore and PacBio data. ### Do I need to trim adapters before BWA-MEM? BWA soft-clips overhanging sequence, so light adapter content is tolerated. Heavy contamination lowers mapping rates, so run FastQC and trim with fastp when adapters are flagged. ## Sources 1. [Getting started in mapping-by-sequencing.](https://pubmed.ncbi.nlm.nih.gov/25359627/) 2. [Comparative bioinformatics for mouse and human genes: getting started.](https://pubmed.ncbi.nlm.nih.gov/18428352/) 3. [NGSphy: phylogenomic simulation of next-generation sequencing data.](https://pubmed.ncbi.nlm.nih.gov/29534152/) 4. [SCSilicon: a tool for synthetic single-cell DNA sequencing data generation.](https://pubmed.ncbi.nlm.nih.gov/35546390/) 5. [Assembly and analysis of stephania japonica mitochondrial genome provides new insights into its identification and energy metabolism.](https://pubmed.ncbi.nlm.nih.gov/39994544/)