保姆级教程：用seqtk、bwa和bedtools从零绘制GC-depth图，诊断测序污染

从零构建GC-depth分析全流程：手把手教你诊断测序数据污染

刚拿到测序数据的生物信息学新手，常常会面临一个灵魂拷问：我的数据干净吗？GC-depth分析就像给测序数据做"体检"，通过一张图就能快速发现细菌污染、样品混杂等潜在问题。本文将用最通俗的语言，带你从fastq文件开始，一步步生成专业的GC-depth图谱。

1. 环境准备与数据检查

工欲善其事，必先利其器。在开始分析前，我们需要准备好以下工具和环境：

• 软件清单：
  • seqtk：用于fastq到fasta格式转换
  • BWA：序列比对工具
  • samtools：处理BAM文件
  • bedtools：基因组区间操作
  • R：绘图分析

推荐使用conda一键安装所有依赖：

conda create -n gcdepth_env -c bioconda seqtk bwa samtools bedtools r-ggplot2 r-argparse conda activate gcdepth_env

常见问题排查：如果遇到权限问题，可以尝试添加--prefix ~/miniconda3/envs/gcdepth_env参数。安装完成后，建议用which seqtk等命令检查各软件是否在PATH中。

2. 数据格式转换与比对

2.1 Fastq到Fasta转换

原始测序数据通常是fastq格式，我们需要先转换为fasta格式：

seqtk seq -A Sample_R1.fastq.gz > Sample_R1.fa seqtk seq -A Sample_R2.fastq.gz > Sample_R2.fa

提示：如果数据量很大，可以添加-C参数压缩输出，节省磁盘空间

2.2 建立参考基因组索引

比对前需要为参考基因组建立索引：

bwa index Sample_R1.fa

2.3 序列比对与排序

使用BWA进行比对并生成排序后的BAM文件：

bwa mem -t 8 Sample_R1.fa Sample_R1.fastq.gz Sample_R2.fastq.gz | \ samtools sort -@ 8 -o aln_sort.bam

参数说明：

• -t 8：使用8个CPU线程
• -@ 8：samtools使用8个线程排序

3. GC-depth分析核心流程

3.1 参数设置与窗口划分

创建参数配置文件config.sh：

#!/bin/bash fa=Sample_R1.fa bam=aln_sort.bam prefix=output window=40 step=10

生成基因组长度文件：

seqtk comp $fa | awk '{print $1"\t"$2}' > $prefix.len

使用bedtools划分分析窗口：

bedtools makewindows -w $window -s $step -g $prefix.len > $prefix.window.bed

3.2 GC含量计算

提取窗口序列并计算GC含量：

seqtk subseq $fa $prefix.window.bed > $prefix.window.fasta seqtk comp $prefix.window.fasta | awk '{print $1 "\t" ($4+$5)/($3+$4+$5+$6) }' > $prefix.window.gc

3.3 测序深度计算

计算每个窗口的平均深度：

bedtools coverage -b $bam -a $prefix.window.bed -mean | \ awk '{print $1":"$2+1"-"$3"\t"$4}' > $prefix.window.depth

4. 可视化分析与结果解读

4.1 R绘图脚本准备

创建gc_depth_plot.r脚本：

library(ggplot2) library(aplot) args <- commandArgs(trailingOnly=TRUE) gc_file <- args[1] depth_file <- args[2] prefix <- args[3] # 数据读取与处理 gc_data <- read.table(gc_file, header=F) depth_data <- read.table(depth_file, header=F) combined <- merge(gc_data, depth_data, by="V1") # 主绘图函数 plot_gcdepth <- function(data, title=""){ ggplot(data, aes(V2.x*100, V2.y)) + geom_point(color="gray60", alpha=0.5, size=0.8) + stat_density_2d(aes(fill=..density..), geom="raster", contour=FALSE) + scale_fill_gradientn(colors=c("transparent","yellow","red")) + labs(x="GC Content (%)", y="Sequencing Depth (X)", title=title) + theme_minimal() } # 生成最终图片 png(paste0(prefix,"_gc_depth.png"), width=8, height=6, units="in", res=300) print(plot_gcdepth(combined)) dev.off()

4.2 运行绘图脚本

执行R脚本生成可视化结果：

Rscript gc_depth_plot.r output.window.gc output.window.depth my_sample

4.3 结果解读指南

正常样本的GC-depth图应该呈现单峰分布，异常情况可能表现为：

5. 进阶技巧与优化建议

5.1 参数优化策略

不同数据类型推荐参数设置：

5.2 批量处理脚本

对于多个样本，可以使用以下批量处理脚本：

#!/bin/bash for sample in $(ls *.fastq.gz | cut -d'_' -f1 | sort -u); do # 处理每个样本 seqtk seq -A ${sample}_R1.fastq.gz > ${sample}.fa bwa mem -t 8 ${sample}.fa ${sample}_R1.fastq.gz ${sample}_R2.fastq.gz | \ samtools sort -@ 8 -o ${sample}.bam # 后续分析步骤... done

5.3 常见报错解决

• BWA索引失败：检查参考基因组是否为fasta格式
• bedtools报错：确保BED文件格式正确，不含header
• R绘图空白：检查输入文件路径是否正确

在实际项目中，我发现窗口大小的设置对结果影响很大。对于人类基因组，窗口40bp配合步长10bp通常能得到理想结果，而植物基因组可能需要更大的窗口。建议先用小样本测试不同参数组合，找到最适合自己数据的配置。