bioinformatics-analysis

Bioinformatics Analysis

Overview

This skill covers standard bioinformatics analysis workflows for sequencing, single-cell, and structural biology data. It complements the exploratory-data-analysis skill (which handles file format detection) with domain-specific analysis guidance.

Use this skill after: Running the EDA skill to understand your data's format and structure.

When to Use This Skill

• Analyzing RNA-seq or scRNA-seq count matrices
• Running differential expression analysis
• Performing gene set enrichment (GO/KEGG)
• Working with VCF variant files
• Analyzing protein structures (PDB files)
• Querying NCBI, Ensembl, UniProt, or STRING databases

---

1. RNA-seq Analysis Pipeline

Step 1: Quality Control

# FastQC for individual files
fastqc sample.fastq.gz -o qc_reports/
# MultiQC to aggregate
multiqc qc_reports/ -o multiqc_report/

What to check in QC reports:

• Per-base quality scores: should be > Q30 across most positions
• Adapter contamination: trim with Trimmomatic or fastp if > 5% reads affected
• GC content: should match expected organism GC content; bimodal suggests contamination
• Duplication rate: > 60% for polyA-selected RNA-seq may indicate issues

Step 2: Alignment

# STAR alignment (recommended for splice-aware alignment)
STAR --runThreadN 8 \
     --genomeDir /path/to/genome_index \
     --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
     --readFilesCommand zcat \
     --outSAMtype BAM SortedByCoordinate \
     --outFileNamePrefix results/sample_

# Alternative: Salmon (quasi-mapping, much faster)
salmon quant -i /path/to/salmon_index -l A \
    -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
    -p 8 -o results/sample_quant

Step 3: Quantification

# featureCounts (for STAR BAM files)
featureCounts -T 8 -p -a genome.gtf \
    -o counts.txt results/*.bam

Step 4: Differential Expression with pyDESeq2

import pandas as pd
import numpy as np
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# Load count matrix (genes × samples)
counts = pd.read_csv("counts.csv", index_col=0)
metadata = pd.read_csv("metadata.csv", index_col=0)
# metadata must have a column matching your design variable

# Create DESeq2 dataset
dds = DeseqDataSet(
    counts=counts.T,  # samples × genes
    metadata=metadata,
    design_factors="condition",  # column in metadata
)
dds.deseq2()

# Run statistical test
stat_res = DeseqStats(dds, contrast=["condition", "treatment", "control"])
stat_res.summary()
results = stat_res.results_df

# Filter significant genes
sig = results[(results["padj"] < 0.05) & (abs(results["log2FoldChange"]) > 1)]
print(f"Significant DEGs: {len(sig)} (padj<0.05, |log2FC|>1)")
sig.to_csv("DEGs.csv")