Featurecounts Rna Counting

# Install with conda (recommended)
conda install -c bioconda subread

# Verify
featureCounts -v
# featureCounts v2.0.6

# Alternative: install via apt (Ubuntu/Debian)
sudo apt-get install subread

decisions:
  - id: D1
    param: annotation
    kind: required
    source: upstream
    ask: "Which GTF/GFF annotation file should be used, and does its assembly exactly match the BAM alignment reference?"
    default: null

  - id: D2
    param: strandSpecific
    kind: required
    source: data
    ask: "Which strandedness did the representative-alignment check infer, and does it agree with the library preparation?"
    default: null

  - id: D3
    param: -t, -g
    kind: optional
    source: user
    ask: "Should expression be summarized from exons to genes, or from another feature type for a specific downstream question?"
    default: "-t exon -g gene_id"

  - id: D4
    param: -M, -O
    kind: optional
    source: user
    ask: "Should reads with multiple genomic mappings or overlapping gene assignments be included in the count matrix?"
    default: "exclude both multi-mappers and multi-overlap reads"

  - id: D5
    param: --minOverlap, --fracOverlap
    kind: optional_conditional
    source: user
    ask: "Does the assay or feature definition require a minimum overlap stricter than one aligned base?"
    default: "--minOverlap 1 --fracOverlap 0"

# Count reads for multiple samples (unstranded paired-end RNA-seq)
featureCounts \
    -a gencode.v47.annotation.gtf \
    -o counts/gene_counts.txt \
    -T 8 \
    -p --countReadPairs \
    results/sample1/Aligned.sortedByCoord.out.bam \
    results/sample2/Aligned.sortedByCoord.out.bam

echo "Count matrix: counts/gene_counts.txt"
head -3 counts/gene_counts.txt

# Verify BAM files are sorted
samtools view -H results/sample1/Aligned.sortedByCoord.out.bam | grep "SO:"
# Expected: SO:coordinate

# List BAMs to count
ls results/*/Aligned.sortedByCoord.out.bam | head -5

# Download GENCODE GTF (same version used for STAR indexing)
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.primary_assembly.annotation.gtf.gz
gunzip gencode.v47.primary_assembly.annotation.gtf.gz

echo "GTF lines: $(wc -l < gencode.v47.primary_assembly.annotation.gtf)"

# Quick strandedness check: count 1 sample with all 3 modes
# Compare assigned rates: highest = correct mode
for strand in 0 1 2; do
    echo "=== Strandedness -s $strand ==="
    featureCounts \
        -a gencode.v47.primary_assembly.annotation.gtf \
        -o /tmp/test_s${strand}.txt \
        -T 4 \
        -p --countReadPairs \
        -s $strand \
        results/sample1/Aligned.sortedByCoord.out.bam 2>&1 \
        | grep "Successfully assigned"
done
# Rule: 0=unstranded if similar rates; 1 or 2 if one is much higher

mkdir -p counts

# Multi-sample batch counting: pass all BAMs as positional arguments
featureCounts \
    -a gencode.v47.primary_assembly.annotation.gtf \
    -o counts/gene_counts.txt \
    -T 8 \
    -p \
    --countReadPairs \
    -s 0 \
    -t exon \
    -g gene_id \
    results/ctrl_1/Aligned.sortedByCoord.out.bam \
    results/ctrl_2/Aligned.sortedByCoord.out.bam \
    results/treat_1/Aligned.sortedByCoord.out.bam \
    results/treat_2/Aligned.sortedByCoord.out.bam

echo "Count matrix: $(wc -l < counts/gene_counts.txt) genes"
# Also generates: counts/gene_counts.txt.summary (assignment stats)
cat counts/gene_counts.txt.summary

# Reverse-stranded library (most TruSeq Stranded protocols): -s 2
featureCounts \
    -a gencode.v47.primary_assembly.annotation.gtf \
    -o counts/gene_counts_stranded.txt \
    -T 8 \
    -p --countReadPairs \
    -s 2 \
    results/*/Aligned.sortedByCoord.out.bam

# Forward-stranded (e.g., Lexogen QuantSeq, Takara SMARTer): -s 1
# featureCounts ... -s 1 ...

echo "Stranded count complete."
head -2 counts/gene_counts_stranded.txt

import pandas as pd

# featureCounts output has 6 metadata columns before count columns
counts_raw = pd.read_csv("counts/gene_counts.txt", sep="\t", comment="#")
print(f"Columns: {list(counts_raw.columns)}")

# Metadata columns: Geneid, Chr, Start, End, Strand, Length
# Count columns start at index 6
count_cols = counts_raw.columns[6:]  # BAM file paths as column names
counts = counts_raw.set_index("Geneid")[count_cols].copy()

# Rename columns to sample names (strip path and file extension)
import re
counts.columns = [re.sub(r".*/|Aligned\.sortedByCoord\.out\.bam", "", col)
                  for col in counts.columns]

print(f"Count matrix shape: {counts.shape}")  # (genes × samples)
print(f"Samples: {list(counts.columns)}")
print(f"Genes with counts > 0: {(counts.sum(axis=1) > 0).sum()}")
counts.to_csv("gene_count_matrix.tsv", sep="\t")
print("Saved: gene_count_matrix.tsv")

import pandas as pd
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
from pydeseq2.ds import DeseqStats

# Load count matrix (genes × samples)
counts = pd.read_csv("gene_count_matrix.tsv", sep="\t", index_col=0).T
print(f"Count matrix: {counts.shape} (samples × genes)")

# Sample metadata
metadata = pd.DataFrame({
    "condition": ["control", "control", "treated", "treated"]
}, index=counts.index)

# Filter low-count genes (recommended before DESeq2)
counts_filtered = counts.loc[:, counts.sum() > 10]
print(f"Genes after low-count filter: {counts_filtered.shape[1]}")

# Run DESeq2
dds = DeseqDataSet(counts=counts_filtered, metadata=metadata,
                   design_factors="condition",
                   inference=DefaultInference(n_cpus=8))
dds.deseq2()

stat_res = DeseqStats(dds, contrast=["condition", "treated", "control"],
                      inference=DefaultInference())
stat_res.summary()
results = stat_res.results_df
sig = results[results["padj"] < 0.05]
print(f"DE genes (padj < 0.05): {len(sig)}")
print(sig.sort_values("log2FoldChange", ascending=False).head())