{"title":"BSA-seq 突变位点定位分析","description":"面向 Bulked Segregant Analysis 的变异检测和候选区间定位流程，覆盖数据下载、质控、比对、GATK SNP calling、过滤、SNP-index 与 ED5 下游分析。","omics_type":"BSA-seq","category_key":"variant","category_name":"基因组变异与定位","dag_json":{"nodes":[{"id":"download","label":"测序数据下载"},{"id":"qc","label":"FastQC/Trimmomatic 质控"},{"id":"align","label":"BWA-MEM 比对 TAIR10"},{"id":"markdup","label":"GATK MarkDuplicates"},{"id":"haplotype","label":"HaplotypeCaller GVCF"},{"id":"combine","label":"CombineGVCFs/GenotypeGVCFs"},{"id":"filter","label":"SNP/INDEL 硬过滤"},{"id":"qtlseq","label":"SNP-index/ED5 下游定位"}],"edges":[{"source":"download","target":"qc"},{"source":"qc","target":"align"},{"source":"align","target":"markdup"},{"source":"markdup","target":"haplotype"},{"source":"haplotype","target":"combine"},{"source":"combine","target":"filter"},{"source":"filter","target":"qtlseq"}]},"metadata_json":{"difficulty":"入门","tools":[],"inputs":["FASTQ","BAM","VCF"],"outputs":[],"estimated_time":"0.5-1 天","scenario":"遗传定位"},"content":"# BSA-seq 突变位点定位分析\n\n## 适用场景\n适用于 Bulked Segregant Analysis（BSA）或 QTL-seq 类型项目：将极端表型个体混池测序，与亲本或对照样本比较等位基因频率差异，从而定位候选突变区间。\n\n本流程来自拟南芥 BSA 分析记录，参考基因组使用 `TAIR10_chr_all.fas`，核心样本包括 `1-37`、`27-3-1`、`27N`、`Neg`、`top1α-10`。流程前半段使用 BWA/GATK 产出高质量 VCF，后半段在 R 中计算 SNP-index、ΔSNP-index 和 ED5，并绘制候选区间图。\n\n## 样本与分组\n| Sample | Role | Notes |\n| --- | --- | --- |\n| `1-37` | bulk / segregant | 用于 37 组合分析 |\n| `27-3-1` | bulk / segregant | 用于 27 组合分析 |\n| `27N` | normal/control | 27 组合对照 |\n| `Neg` | wild type/control | 37 组合对照 |\n| `top1α-10` | parent/reference mutant | 亲本或参考突变材料 |\n\n## 1. 下载原始数据\n```bash\n./lnd login -u X101SC250910114-Z01-J003 -p ggeah41n\n\n./lnd cp -d   oss://CP2023071800097/H101SC250910114/RSPD00204/X101SC250910114-Z01/X101SC250910114-Z01-J003/01.RawData   ../../../01_raw_data/\n```\n\n## 2. FastQC 质控与 Trimmomatic 裁剪\n```bash\nmodule load FastQC/0.11.9\n\nbsub -J fastqc1-37 -n 1 -R \"span[hosts=1]\" -o %J.out -e %J.err -q normal \"fastqc 1-37_*\"\nbsub -J fastqc27-3 -n 1 -R \"span[hosts=1]\" -o %J.out -e %J.err -q normal \"fastqc 27-3*\"\nbsub -J fastqc27N -n 1 -R \"span[hosts=1]\" -o %J.out -e %J.err -q normal \"fastqc 27N*\"\nbsub -J fastqcNeg -n 1 -R \"span[hosts=1]\" -o %J.out -e %J.err -q normal \"fastqc Neg*\"\nbsub -J fastqctop1 -n 1 -R \"span[hosts=1]\" -o %J.out -e %J.err -q normal \"fastqc top1α*\"\n```\n\n```bash\nfor sample in 1-37 27-3-1 27N Neg top1α-10\ndo\n  trimmomatic.sh PE -threads 10 -phred33     01_raw_data/${sample}_1.fq.gz     01_raw_data/${sample}_2.fq.gz     02_clean_data/${sample}_R1_paired.fq.gz     02_clean_data/${sample}_R1_unpaired.fq.gz     02_clean_data/${sample}_R2_paired.fq.gz     02_clean_data/${sample}_R2_unpaired.fq.gz     HEADCROP:2\ndone\n```\n\n## 3. BWA-MEM 比对到 TAIR10\n```bash\nREF=\"./ref_genomic/TAIR10_chr_all.fas\"\n\nfor sample in 1-37 27-3-1 27N Neg top1α-10\ndo\n  bwa mem -M -t 8     -R \"@RG\tID:${sample}\tSM:${sample}\"     ${REF}     02_clean_data/${sample}_R1_paired.fq.gz     02_clean_data/${sample}_R2_paired.fq.gz |     samtools sort -@ 8 -o 03_alignment/${sample}.sorted.bam -\n\n  samtools flagstat 03_alignment/${sample}.sorted.bam > 03_alignment/${sample}_flagstat.log\n  samtools index 03_alignment/${sample}.sorted.bam\ndone\n```\n\n## 4. GATK MarkDuplicates 去重\n```bash\nfor sample in 1-37 27-3-1 27N Neg top1α-10\ndo\n  gatk MarkDuplicates     -I 03_alignment/${sample}.sorted.bam     -M 03_alignment/${sample}.markup_metrics.txt     -O 03_alignment/${sample}.sorted.markup.bam\n\n  samtools index 03_alignment/${sample}.sorted.markup.bam\ndone\n```\n\n## 5. HaplotypeCaller 生成 GVCF\n```bash\nREF=\"./ref_genomic/TAIR10_chr_all.fas\"\n\nfor sample in 1-37 27-3-1 27N Neg top1α-10\ndo\n  gatk --java-options \"-Xmx10G\" HaplotypeCaller     -R ${REF}     -I 03_alignment/${sample}.sorted.markup.bam     -O 04_snp_calling/${sample}.g.vcf.gz     -ERC GVCF     -stand-call-conf 10\ndone\n```\n\n## 6. 合并 GVCF 与联合分型\n```bash\nREF=\"./ref_genomic/TAIR10_chr_all.fas\"\n\ngatk CombineGVCFs   -R ${REF}   -V 04_snp_calling/1-37.g.vcf.gz   -V 04_snp_calling/Neg.g.vcf.gz   -V 04_snp_calling/top1α-10.g.vcf.gz   -O 04_snp_calling/BSA-37.gvcf.gz\n\ngatk CombineGVCFs   -R ${REF}   -V 04_snp_calling/27-3-1.g.vcf.gz   -V 04_snp_calling/27N.g.vcf.gz   -V 04_snp_calling/top1α-10.g.vcf.gz   -O 04_snp_calling/BSA-27.gvcf.gz\n\ngatk GenotypeGVCFs   -R ${REF}   -V 04_snp_calling/BSA-27.gvcf.gz   -O 05_snp/BSA_no_filter-27.HC.vcf\n```\n\n## 7. SNP/INDEL 硬过滤\n```bash\ngatk SelectVariants   -R ./ref_genomic/TAIR10_chr_all.fas   -V 05_snp/BSA_no_filter-27.HC.vcf   -select-type SNP   -O 05_snp/SNPs-27.vcf\n\ngatk VariantFiltration   -V 05_snp/SNPs-27.vcf   --filter-expression \"MQ < 30.0\"   --filter-name \"MQ_filter_SNP\"   --filter-expression \"QD < 2.0\"   --filter-name \"QD_filter_SNP\"   -O 05_snp/SNPs_filter-27.vcf\n\ngrep -E '^#|PASS' 05_snp/SNPs_filter-27.vcf > 05_snp/BSA_SNPs_filterPASSED-27.vcf\n```\n\nINDEL 可继续使用 `SelectVariants -select-type INDEL` 后按 `MQ`、`SOR`、`QD`、`FS` 进行硬过滤，并与 SNP 结果合并为最终 VCF。\n\n## 8. R 下游：SNP-index 与 ED5 定位\n```r\nlibrary(QTLseqr)\nlibrary(vcfR)\nlibrary(ggplot2)\nlibrary(dplyr)\n\nvcf <- read.vcfR(\"BSA_SNPs_filterPASSED-27.vcf\")\nchrom <- getCHROM(vcf)\npos <- getPOS(vcf)\nref <- getREF(vcf)\nalt <- getALT(vcf)\nad <- extract.gt(vcf, \"AD\")\ngt <- extract.gt(vcf, \"GT\")\n\nref_split <- masplit(ad, record = 1, sort = 0)\nalt_split <- masplit(ad, record = 2, sort = 0)\n\ndf <- data.frame(\n  CHROM = chrom,\n  POS = pos,\n  REF = ref,\n  ALT = alt,\n  AD_REF.1_37 = ref_split[, 1],\n  AD_ALT.1_37 = alt_split[, 1],\n  AD_REF.Neg = ref_split[, 2],\n  AD_ALT.Neg = alt_split[, 2],\n  AD_REF.top1 = ref_split[, 3],\n  AD_ALT.top1 = alt_split[, 3],\n  P_Top1 = gt[, \"top1α-10\"],\n  son_1_37 = gt[, \"1-37\"],\n  son_Neg = gt[, \"Neg\"]\n)\n\ndf <- df %>%\n  filter(\n    nchar(as.character(REF)) == 1,\n    nchar(as.character(ALT)) == 1,\n    !is.na(P_Top1),\n    (AD_REF.1_37 + AD_ALT.1_37) >= 10,\n    (AD_REF.Neg + AD_ALT.Neg) >= 10,\n    (AD_REF.top1 + AD_ALT.top1) >= 10\n  )\n```\n\n```r\ndf_result <- df %>%\n  mutate(\n    SNPindex.L = AD_ALT.1_37 / (AD_REF.1_37 + AD_ALT.1_37),\n    SNPindex.S = AD_ALT.Neg / (AD_REF.Neg + AD_ALT.Neg),\n    deltaSNP = SNPindex.L - SNPindex.S\n  )\n\ndelta_extreme_points <- df_result %>%\n  filter(deltaSNP > 0.5 | deltaSNP < -0.5)\n\nggplot(df_result, aes(x = POS / 1e6, y = deltaSNP)) +\n  geom_point(color = \"#000000\", size = 0.5, alpha = 0.6) +\n  geom_point(\n    data = delta_extreme_points,\n    aes(color = ifelse(deltaSNP > 0, \"Positive\", \"Negative\")),\n    size = 1.2,\n    alpha = 0.8\n  ) +\n  geom_hline(yintercept = c(-0.5, 0.5), linetype = \"dashed\") +\n  facet_wrap(~ CHROM, scales = \"free_x\", nrow = 3) +\n  labs(x = \"Genomic Position (Mb)\", y = \"Delta SNP-index\")\n```\n\n```r\nvalid_rows <- complete.cases(df) &\n  df$AD_REF.1_37 + df$AD_ALT.1_37 > 0 &\n  df$AD_REF.Neg + df$AD_ALT.Neg > 0\n\ndf_filtered <- df[valid_rows, ]\n\nED_values <- apply(df_filtered, 1, function(row) {\n  ref_mutant <- as.numeric(row[\"AD_REF.1_37\"])\n  alt_mutant <- as.numeric(row[\"AD_ALT.1_37\"])\n  ref_wild <- as.numeric(row[\"AD_REF.Neg\"])\n  alt_wild <- as.numeric(row[\"AD_ALT.Neg\"])\n\n  depth_mutant <- ref_mutant + alt_mutant\n  depth_wild <- ref_wild + alt_wild\n\n  freq_ref_mutant <- ref_mutant / depth_mutant\n  freq_alt_mutant <- alt_mutant / depth_mutant\n  freq_ref_wild <- ref_wild / depth_wild\n  freq_alt_wild <- alt_wild / depth_wild\n\n  freq_diff <- sqrt(\n    (freq_ref_wild - freq_ref_mutant)^2 +\n      (freq_alt_wild - freq_alt_mutant)^2\n  )\n\n  freq_diff^5\n})\n\nresult_df <- data.frame(\n  CHROM = df_filtered$CHROM,\n  POS = df_filtered$POS,\n  ED = ED_values\n)\n```\n\n## 关键质控与解释\n- FASTQ 阶段关注碱基质量、接头残留和头部异常碱基，当前流程使用 `HEADCROP:2`。\n- 比对阶段保留 read group，方便 GATK 正确识别样本。\n- GVCF 阶段采用每个样本单独 HaplotypeCaller，再按组合合并和联合分型。\n- SNP 过滤至少包含 `MQ >= 30` 与 `QD >= 2`；INDEL 过滤可加入 `SOR` 和 `FS`。\n- 下游定位同时看 `Delta SNP-index` 和 `ED5`，比只看单一指标更稳。\n\n## 主要输出\n- FastQC 质控报告\n- Trimmomatic clean reads\n- 排序和去重后的 BAM\n- 样本级 GVCF\n- 合并分型 VCF\n- 过滤后的 SNP/INDEL VCF\n- SNP-index / Delta SNP-index 图\n- ED5 全基因组定位图\n","id":23,"created_at":"2026-06-03T17:48:59.764103Z"}