tooluniverse-sequence-analysis

# 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(term="TP53[Symbol] AND Homo sapiens[Organism]")` -> "7157" 2. `NCBIDatasets_get_orthologs(gene_id="7157", page_size=10)` -> mouse Trp53, rat Tp53, etc. ## Phase 5: Domain Architecture and Homology **InterPro_get_entries_for_protein**: `accession` (UniProt ID). Returns InterPro domain/family/superfamily entries with positions. **Pfam_get_protein_annotations**: `accession` (UniProt ID). Returns Pfam domain hits with exact residue coordinates and E-values. For a raw sequence with no UniProt accession yet (e.g. a novel or designed sequence), `EBI_scan_pfam_domains` runs the same Pfam scan directly on the sequence instead. **BLAST_protein_search**: `sequence` (amino acid string), `database` (default "swissprot"), `limit`. Returns homologs with alignment scores, identity, E-values. If BLAST finds nothing (remote homology, low sequence identity), `EBI_profile_search` (`method='phmmer'` default, or `'psiblast'`) uses profile/HMM methods that are more sensitive for distant relationships BLAST's pairwise scoring misses. **EBI_pairwise_align**: exactly two sequences (`algorithm='needle'` for global, `'water'` for local), returns percent identity/similarity/gaps and the alignment itself -- no existing tool in this skill does direct pairwise alignment. **EBI_translate_sequence**: nucleotide-to-protein, protein-to-nucleotide (back-translation), or 6-frame translation with ORF detection. `DNA_translate_reading_frames` above only does forward 3-frame (no reverse strand, no back-translation) -- use EBI's tool specifically for reverse-strand frames or when going protein-to-DNA. **EBI_predict_membrane_topology**: predicts transmembrane helices and signal peptides directly from a protein sequence (Phobius). Distinct from `proteins_api_get_features`'s TRANSMEM annotations (see "Reasoning for Protein Feature Questions" below) -- that's looking up *already-annotated* topology for a known UniProt entry; this predicts topology for a sequence that may have none. **EnsemblCompara_get_orthologues**: `gene` (gene symbol, e.g., "CFTR"), `species` (e.g., "human"). User-friendly alternative to NCBIDatasets_get_orthologs — accepts gene symbols directly. ## Phase 6: Variant and Clinical Context **EnsemblVEP_annotate_hgvs**: `hgvs_notation` (e.g., "NM_000492.4:c.1521_1523del"). Returns consequence, protein impact, genomic coordinates. **ClinVar_search_variants**: `gene` (gene symbol). Returns variant count and IDs for clinical significance lookup. **PubMed_search_articles**: `query`, `limit`. Literature context for gene/variant findings. --- ## Tool Parameter Quick Reference | Tool | Correct Param | Common Mistake | |------|--------------|----------------| | NCBIGene_search | `term` (with [Symbol] syntax) | `query` or `gene` | | NCBIGene_get_summary | `id` (string) | Integer type | | NCBI_fetch_accessions | `uids` (array) | `accessions` | | NCBI_get_sequence | `accession` (string) | Passing UID | | NCBIDatasets_get_orthologs | `gene_id` (string) | Gene symbol | | EnsemblSeq_get_id_sequence | `ensembl_id` | `id` | | ensembl_get_sequence | `id` + `multiple_sequences` | Omitting multiple_sequences for gene+CDS | | UniProt_get_sequence_by_accession | `accession` | Response is `result` not `data` | ## Fallbacks - Gene not found -> try NCBIDatasets_get_gene_by_symbol with explicit taxon - No accessions from search -> broaden query (remove strain/seq_type filters) - Ensembl error for gene+CDS -> use transcript ID (ENST) or set multiple_sequences=true - UniProt accession unknown -> NCBIDatasets_get_gene or UniProt_search for cross-refs - Ortholog search empty -> verify gene_id is numeric NCBI Gene ID ## Sequence Analysis Reasoning (CRITICAL) **LOOK UP DON'T GUESS** -- always fetch sequences, coordinates, and domain boundaries from databases. Do not reconstruct them from memory. ### When to Use Which Tool | Question Type | Tool Choice | Why | |--------------|------------|-----| | "Find similar sequences" | BLAST_protein_search | Homology search against databases; returns E-values and identity | | "What domains does this protein have?" | InterPro_get_entries_for_protein or Pfam_get_protein_annotations | Domain architecture with exact residue coordinates | | "Get the sequence of gene X" | NCBI_search_nucleotide -> NCBI_get_sequence | Nucleotide retrieval by gene name | | "Compare orthologs" | NCBIDatasets_get_orthologs or EnsemblCompara_get_orthologues | Cross-species gene comparison | | "What is the protein impact of variant X?" | EnsemblVEP_annotate_hgvs | Consequence prediction with protein coordinates | | "Align two sequences" | BLAST (pairwise) | Quick pairwise comparison with scoring |