参考基因组下载与索引准备：Ensembl/UCSC/NCBI

1
#!/bin/bash
2
# download_ref.sh —— 完整下载 Ensembl GRCh38 + GTF 注释
3

4
RELEASE=113
5
SPECIES="homo_sapiens"
6
ASSEMBLY="GRCh38"
7
BASE="ftp://ftp.ensembl.org/pub/release-${RELEASE}"
8
OUTDIR="./ref/GRCh38_ensembl"
9

10
mkdir -p "$OUTDIR"
11

12
# 1. 基因组 FASTA
13
echo "Downloading genome FASTA..."
14
axel -n 10 -o "$OUTDIR" \
15
    "${BASE}/fasta/${SPECIES}/dna/${SPECIES^}.${ASSEMBLY}.dna.primary_assembly.fa.gz"
16

17
# 2. GTF 注释文件
18
echo "Downloading GTF annotation..."
19
axel -n 10 -o "$OUTDIR" \
20
    "${BASE}/gtf/${SPECIES}/${SPECIES^}.${ASSEMBLY}.${RELEASE}.gtf.gz"
21

22
# 3. 校验
23
echo "Downloading CHECKSUMS..."
24
wget -q -O "$OUTDIR/CHECKSUMS" "${BASE}/fasta/${SPECIES}/dna/CHECKSUMS"
25

26
echo "Done! Files in $OUTDIR"
27
ls -lh "$OUTDIR"

1
#!/bin/bash
2
# build_all_indices.sh —— 一键建好所有比对器索引
3
set -euo pipefail
4

5
GENOME="./ref/GRCh38_ensembl/Homo_sapiens.GRCh38.dna.primary_assembly.fa"
6
THREADS=16
7
OUTDIR="./ref/GRCh38_ensembl/indices"
8
mkdir -p "$OUTDIR"
9

10
echo "============================================"
11
echo "1/5: samtools faidx"
12
echo "============================================"
13
samtools faidx "$GENOME"
14
# 生成 .fai 索引文件
15

16
echo "============================================"
17
echo "2/5: GATK CreateSequenceDictionary"
18
echo "============================================"
19
# GATK 需要 .dict 文件（本质是 SAM header 格式的序列字典）
20
gatk CreateSequenceDictionary \
21
    -R "$GENOME" \
22
    -O "${GENOME%.fa}.dict"
23
# 或用 picard：
24
# picard CreateSequenceDictionary R="$GENOME" O="${GENOME%.fa}.dict"
25

26
echo "============================================"
27
echo "3/5: bwa index (BWA-MEM)"
28
echo "============================================"
29
# BWA 索引——建在基因组同目录
30
bwa index -p "${OUTDIR}/bwa/GRCh38" "$GENOME"
31
# 生成 .amb .ann .bwt .pac .sa 五个文件
32

33
echo "============================================"
34
echo "4/5: bowtie2-build"
35
echo "============================================"
36
bowtie2-build --threads "$THREADS" \
37
    "$GENOME" \
38
    "${OUTDIR}/bowtie2/GRCh38"
39

40
echo "============================================"
41
echo "5/5: hisat2-build (含 splice site + exon 信息)"
42
echo "============================================"
43
# hisat2 最好加上已知的 splice site 和 exon 信息提高比对精度
44
GTF="./ref/GRCh38_ensembl/Homo_sapiens.GRCh38.113.gtf.gz"
45

46
# 先从 GTF 提取剪接位点和外显子
47
hisat2_extract_splice_sites.py "$GTF" > "${OUTDIR}/hisat2/splicesites.txt"
48
hisat2_extract_exons.py "$GTF" > "${OUTDIR}/hisat2/exons.txt"
49

50
# 建索引（含已知剪接信息）
51
hisat2-build -p "$THREADS" \
52
    --ss "${OUTDIR}/hisat2/splicesites.txt" \
53
    --exon "${OUTDIR}/hisat2/exons.txt" \
54
    "$GENOME" \
55
    "${OUTDIR}/hisat2/GRCh38"
56

57
echo "============================================"
58
echo "All indices built!"
59
echo "============================================"
60

61
# 索引文件大小统计
62
echo ""
63
echo "Index sizes:"
64
du -sh "${OUTDIR}/"*

1
#!/bin/bash
2
# verify_indices.sh —— 验证所有索引是否可用
3

4
GENOME="./ref/GRCh38_ensembl/Homo_sapiens.GRCh38.dna.primary_assembly.fa"
5

6
echo "=== samtools faidx ==="
7
samtools faidx "$GENOME" chr1:1-100
8
# 应输出 chr1 前 100bp 序列
9

10
echo "=== bwa ==="
11
echo -e "@test\nACGTACGTACGT\n+\nIIIIIIIIIIII" > /tmp/test.fq
12
bwa mem -t 2 ./ref/GRCh38_ensembl/indices/bwa/GRCh38 /tmp/test.fq 2>&1 | head -5
13
# 不报错即可（未必能比对到）
14

15
echo "=== bowtie2 ==="
16
bowtie2 -x ./ref/GRCh38_ensembl/indices/bowtie2/GRCh38 -U /tmp/test.fq 2>&1 | head -5
17

18
echo "=== GATK dict ==="
19
head -5 "${GENOME%.fa}.dict"
20
# 应输出 @HD @SQ 等 SAM header
21

22
echo "=== Validation complete ==="
23
rm -f /tmp/test.fq