A sequence alignment pipeline for ChIP-seq data using BWA, developed for the "Bioinformatics II" module in the Genomic Sciences program at UNAM Juriquilla.
This project implements a complete alignment pipeline including:
- Reference genome indexing with BWA
- Read alignment for ChIP-seq experiments
- Quality control of raw sequences with FastQC
- Alignment quality assessment with samtools
The analysis uses two datasets:
- Reference genome:
Escherichia_coli_K12_MG1655.fasta - ChIP-seq reads:
SRX189773_FNR_ChIP.fastq(FNR transcription factor)
- Reference genome: GRCm38 (
mm10.fa) - ChIP-seq reads:
mus_musculus_CEBPA_liver_ERR005132.fastq.gz(CEBPA, liver tissue)
| Step | Command/Script | Description |
|---|---|---|
| 1 | bwa index |
Generate genome index files (.amb, .ann, .bwt, .pac, .sa) |
| 2 | fastqc |
Quality control of raw reads |
| 3 | bwa aln + bwa samse |
Align single-end reads (short reads) |
| 3 (alt) | bwa mem |
Align reads using MEM algorithm (longer reads) |
| 4 | samtools stats/flagstat |
Alignment quality assessment |
bwa index -p E_index -a bwtsw Escherichia_coli_K12_MG1655.fasta# Generate alignment candidates
bwa aln E_index SRX189773_FNR_ChIP.fastq > Escherichia.sai
# Convert to SAM format
bwa samse E_index Escherichia.sai SRX189773_FNR_ChIP.fastq > Escherichia.sambwa mem -M -t 8 /path/to/mm10.fa reads.fastq.gz > Mus_alg.samsamtools stats alignment.sam | grep 'SN'
samtools flagstat alignment.sam| Metric | Value |
|---|---|
| Total reads | 3,603,544 |
| Mapped reads | 2,351,791 (65%) |
| Multi-mapped (MQ0) | 85,370 (2.4%) |
| Error rate | 0.2% |
| Mean quality | 34.6 (Phred) |
| Metric | Value |
|---|---|
| Read length | 42 bp |
| GC content | 45% |
| Mapped reads | 6,530,881 |
| Unmapped reads | 10,640,249 |
| Error rate | 2.1% |
Note: The M. musculus dataset showed significant quality issues (mean Phred ~18.7), likely due to sequencing errors. Results should be interpreted with caution.
# Load modules (cluster environment)
module load bwa/0.7.15
module load samtools/1.9
module load fastqc/0.11.3- Set up working directory and copy required files
- Index reference genome (if not pre-indexed)
- Run FastQC on raw reads
- Perform alignment with BWA
- Assess alignment quality with samtools
For large genomes (M. musculus), submit jobs via SGE:
qsub MusculusAlignment.sge*.sai: BWA alignment index files*.sam: Sequence Alignment/Map format files*_fastqc.html: FastQC quality reports- Alignment statistics from samtools
Jessica Garcia, Manuel Rivera, Axel Rodriguez Perez
Genomic Sciences, UNAM · February 2023