tooluniverse-sequence-analysis

What this skill does

# Biological Sequence Analysis

## ⚠️ TOP-OF-MIND RULE: Trimmomatic "reads completely discarded" = `F + R + 2*D`, summed across samples

When a question asks about Trimmomatic "reads completely discarded" / "reads thrown out" /
"reads not in any output", do NOT report the `Dropped` field alone. `Dropped` counts
PAIRS where both mates failed; each pair = 2 individual reads. Plus the Forward-only and
Reverse-only buckets also discard one read per pair.

```
reads_discarded = sum over samples of (Forward_Only + Reverse_Only + 2 * Dropped)
```

❌ WRONG: `sum(Dropped per sample)` — typically reports ~thousands, GT is 100×+ higher

✅ RIGHT: `sum(F + R + 2*D per sample)`

Full counter table is in the FASTQ section below.

---

## RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:
- `*_executed.ipynb` → read with `tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}'` and cite its cell outputs as the authoritative answer
- Pre-computed result files (CSV/TSV with names like `*results*`, `*deseq*`, `*enrich*`, `*stats*`, `*_simplified.csv`) → read directly and report the requested value
- Canonical analysis scripts (`analysis.R`, `run_*.py`, `find_*.R`, `*.Rmd`) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if **none** of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

---

Retrieve, annotate, and compare biological sequences from NCBI, Ensembl, and UniProt. Covers nucleotide search, sequence fetching, gene summaries, ortholog discovery, and protein sequence extraction.

## FASTQ QC, Trimmomatic, and read alignment (when raw reads are present)

When the data folder has `*.fastq` files and the question involves Trimmomatic, BWA, samtools, FastQC, or coverage depth, this skill is the entry point — but the actual work is shell-level (no specific ToolUniverse data tool).

### Trimmomatic PE — counting "completely discarded" reads

Trimmomatic PE classifies each input pair as:
- **Both Surviving** (B): R1 and R2 both pass → kept as paired
- **Forward Only Surviving** (F): R1 passes, R2 dropped → R1 kept as singleton, R2 fully discarded
- **Reverse Only Surviving** (R): R2 passes, R1 dropped → R2 kept as singleton, R1 fully discarded
- **Dropped** (D): both fail → BOTH R1 and R2 fully discarded

Counter selection — read the question carefully. **CRITICAL: "READS completely discarded" ≠ "PAIRS dropped".** The Trimmomatic Dropped count counts PAIRS (each = 2 individual reads). When the question asks about reads (not pairs), translate every counter to per-read terms:

| Question phrasing | Formula (per-sample, then SUM across all samples) |
|---|---|
| "**reads completely discarded**", "reads thrown out", "reads not in any output" | `F + R + 2*D` (every individual R1 or R2 not in any output FASTQ) |
| "**read pairs dropped**", "pairs where both mates failed" | `D` |
| "individual R2 reads dropped" (R1 kept as singleton) | `F` |
| "individual R1 reads dropped" (R2 kept as singleton) | `R` |
| "reads passing QC" / "surviving reads" | `2*B + F + R` |

Trimmomatic's stderr summary gives `Input Read Pairs: N Both Surviving: B (b%) Forward Only Surviving: F (f%) Reverse Only Surviving: R (r%) Dropped: D (d%)`. Always sum across ALL input sample pairs (e.g., SRR1 + SRR2 + ...).

DO NOT report just `D` as "reads completely discarded" — that's pair count, not read count, and is off by ~100×. The "Forward Only" R2 mate IS discarded; the "Reverse Only" R1 mate IS discarded; "Dropped" pairs lose BOTH reads.

### Coverage depth (samtools depth / mosdepth)

For "average coverage depth", run `samtools depth -a alignment.bam | awk '{sum+=$3; n++} END {print sum/n}'` — the `-a` flag includes positions with zero coverage (otherwise the average is biased upward). For per-chromosome coverage, group by `$1`.

## When to Use

- "Get the mRNA sequence for BRCA1"
- "Search NCBI for E. coli K-12 complete genome"
- "Find orthologs of TP53 across species"
- "Fetch the protein sequence for UniProt P04637"
- "Get the CDS sequence for Ensembl transcript ENST00000269305"

## Workflow

```
Input -> Phase 1: Gene ID resolution -> Phase 2: Nucleotide retrieval
      -> Phase 3: Protein sequences -> Phase 4: Orthologs -> Output
```

## Phase 1: Gene Identification and Summary

**NCBIGene_search**: `term` (string REQUIRED, format `"TP53[Symbol] AND Homo sapiens[Organism]"`), `retmax` (int, default 10). Returns `{status, data: {esearchresult: {idlist: ["7157"]}}}`.

**NCBIGene_get_summary**: `id` (string REQUIRED, e.g., "7157"). Returns `{status, data: {result: {"7157": {name, description, summary, chromosome, maplocation, genomicinfo, mim}}}}`. Result is keyed by gene ID string.

**NCBIDatasets_get_gene_by_symbol**: `symbol` (string REQUIRED, e.g., "BRCA1"), `taxon` (string, e.g., "human"). Returns gene ID, description, location, cross-references.

**NCBIDatasets_get_gene**: `gene_id` (string REQUIRED, e.g., "7157"). Returns comprehensive gene info.

## Phase 2: Nucleotide Sequence Search and Retrieval

**NCBI_search_nucleotide**: `query` (free-form), `organism` (string), `gene` (string), `strain` (string), `keywords` (string), `seq_type` ("complete_genome"/"mRNA"/"refseq"), `limit` (int, default 20). Returns `{status, data: {uids: [...], accessions: [...]}}`.

**NCBI_fetch_accessions**: `uids` (array REQUIRED, e.g., ["545778205"]). Returns `{status, data: ["U00096.3"], count: 1}`.

**NCBI_get_sequence**: `accession` (string REQUIRED, e.g., "NM_007294"), `format` ("fasta"/"gb"/"embl"). Returns `{status, data: "FASTA string...", accession, format, length}`.

**EnsemblSeq_get_region_sequence**: `region` (string REQUIRED, "chr:start-end", e.g., "17:7668421-7668520"), `species` (default "homo_sapiens"). Returns `{status, data: {sequence, sequence_length}}`.

**ensembl_get_sequence**: `id` (string REQUIRED, Ensembl ID), `type` ("genomic"/"cds"/"cdna"/"protein"), `multiple_sequences` (bool). Returns sequence data.

**Gotchas**:
- NCBI_search_nucleotide returns UIDs, not accessions. Use NCBI_fetch_accessions to convert.
- NCBI_fetch_accessions requires `uids` (NOT `accessions`).
- ensembl_get_sequence with gene ID (ENSG) + type != "genomic" requires `multiple_sequences=true`. Use transcript IDs (ENST) for specific sequences.

### Recipe: Get mRNA for a human gene
1. `NCBI_search_nucleotide(organism="Homo sapiens", gene="BRCA1", seq_type="mRNA", limit=5)`
2. `NCBI_fetch_accessions(uids=[first_uid])` -> accession
3. `NCBI_get_sequence(accession="NM_007294", format="fasta")`

## Phase 3: Protein Sequence Retrieval

**UniProt_get_sequence_by_accession**: `accession` (string REQUIRED, e.g., "P04637"). Returns `{result: "MEEPQSDP..."}`. **Note**: response key is `result`, NOT `data`.

**EnsemblSeq_get_id_sequence**: `ensembl_id` (string REQUIRED, e.g., "ENSP00000269305"), `type` ("protein"/"cdna"/"cds"). Returns `{status, data: {ensembl_id, molecule, sequence, sequence_length}}`.

**UniProt_get_entry_by_accession**: `accession` (string REQUIRED). Full protein annotation.

**Gotchas**:
- UniProt_get_sequence_by_accession returns `{result: "..."}`, not `{status, data}`.
- For Ensembl protein seqs, use ENSP IDs. For cDNA/CDS, use ENST IDs.
- To find UniProt accession from gene: use NCBIDatasets_get_gene_by_symbol (has cross-refs).

## Phase 4: Ortholog and Comparative Analysis

**NCBIDatasets_get_orthologs**: `gene_id` (string REQUIRED, NCBI Gene ID e.g., "7157"), `page_size` (int, default 20, max 100). Returns `{status, data: [{gene_id, symbol, description, taxname, common_name, chromosomes}]}`.

**NCBIProtein_get_summary**: `id` (string REQUIRED, GI number or accession). Returns protein title, organism, length.

**Gotcha**: NCBIDatasets_get_orthologs requires NCBI Gene ID (numeric string), not gene symbol or Ensembl ID. Resolve via Phase 1 first.

### Recipe: Compare orthologs
1. `NCBIGene_search