**基于Python的高通量测序数据质量控制与可视化全流程实战**在生物信息学

基于Python的高通量测序数据质量控制与可视化全流程实战

在生物信息学领域，高通量测序（Next-Generation Sequencing, NGS）已成为基因组研究的核心技术之一。然而，原始测序数据往往包含低质量碱基、接头污染、序列偏倚等问题，若不进行严格的质量控制（QC），后续分析结果可能严重失真。本文将以Python语言为核心工具，结合FastQC、MultiQC和自定义脚本，构建一套完整的NGS数据质量评估与可视化流程。

🔍 一、整体流程图示意

原始FASTQ文件 → FastQC检测 → 多样本汇总(MultiQC) → 自定义Python清洗/过滤 → 输出高质量数据 ↑ 批量处理脚本（bash/python） ``` 此流程可轻松扩展至数百个样本，适用于WGS、RNA-seq、ChIP-seq等多种实验设计。 --- ### 🧪 二、使用FastQC快速生成单样本质量报告 首先，我们用命令行调用 FastQC 对每个 FASTQ 文件执行基础质量评估： ```bash fastqc sample1_R1.fastq.gz -o ./qc_reports/ fastqc sample1_R2.fastq.gz -o ./qc_reports/

⚠️ 建议将所有样本放在同一目录下，并统一命名规则如sample_X_R1.fq.gz，便于后续自动化处理。
FastQC 输出为 HTML 报告和 JSON 格式元数据，可用于进一步解析。例如，获取每个碱基位置的平均质量得分（Phred Score）：

importjsonimportosdefparse_fastqc_json(json_path):withopen(json_path,'r')asf:data=json.load(f)# 获取碱基质量均值per_base_quality=data['Basic Statistics']['Per base sequence quality']avg_quality=[round(float(q),2)forqinper_base_quality]returnavg_quality# 示例调用json_file="./qc_reports/sample1_R1_fastqc.json"quality_scores=parse_fastqc_json(json_file)print("每碱基平均质量分数:",quality_scores[:10])# 显示前10位

输出示例：

每碱基平均质量分数: [35.7, 36.1, 36.5, 37.0, 37.4, 37.8, 38.1, 38.4, 38.6, 38.8]

这一步帮助我们判断是否存在质量下降趋势（如末端碱基质量显著降低），从而决定是否需要剪切或过滤。

📊 三、多样本整合：利用 MultiQC 构建综合质量概览

当项目涉及多个样本时，手动查看每个 FastQC 报告效率低下。此时推荐使用MultiQC自动生成合并报表：

multiqc ./qc_reports/-o./multiqc_output/

该命令会自动收集所有.fastqc.zip或.html文件，生成交互式网页报告，包含以下关键指标：

• 每个样本的GC含量分布
• • 接头污染比例（Adapter Content）
• • 序列长度分布（Sequence Length Distribution）
• • 碱基错误率（Base Quality Distribution）
    📌 在实际项目中，建议设置阈值：
• 平均Phred > 30（表示错误率 < 0.1%）
• • 接头污染 < 5%
• • GC含量偏离参考基因组 ±10%
    若某样本不达标，应进入下一步清洗环节。

🛠 四、Python定制化清洗流程（基于 Trimmomatic 思想）

虽然 Trimmomatic 是主流工具，但有时我们需要更灵活的控制策略。这里提供一个轻量级 Python 脚本，用于截断低质量尾部（Threshold=Q20）并去除短于 36bp 的读段：

fromBioimportSeqIOimportsysdeftrim_low_quality_reads(input_fq,output_fq,min_length=36,quality_threshold=20):trimmed_count=0total_count=0withopen(output_fq,'w')asout_handle:forrecordinSeqIO.parse(input_fq,"fastq"):total_count+=1qual_scores=record.letter_annotations["phred_quality"]# 找到第一个低于阈值的位置cut_point=next((ifori,qinenumerate(qual_scores)ifq<quality_threshold),len(qual_scores))ifcut_point>0andlen(record.seq)>=min_length:trimmed_seq=record.seq[:cut_point]trimmed_qual=qual_scores[:cut_point]new_record=record[:]new_record.seq=trimmed_seq new_record.letter_annotations["phred_quality"]=trimmed_qual SeqIO.write(new_record,out_handle,"fastq")trimmed_count+=1print(f"Processed{total_count}reads, kept{trimmed_count}after trimming.")``` 调用方式： ```python trim_low_quality_reads("sample1_R1.fastq.gz","sample1_R1_trimmed.fastq.gz")

✅ 此脚本优势在于：

• 可嵌入Pipeline（如 Snakemake 或 Nextflow）
• • 支持参数动态配置（如调整 threshold、min_length）
• • 无需额外依赖外部工具（仅需 Biopython）

📈 五、进阶技巧：绘制质量变化热力图（Matplotlib + Seaborn）

对于批量样本，我们可以将各样本的 per-base quality 绘制成热力图，直观比较质量差异：

importseabornassnsimportmatplotlib.pyplotaspltimportnumpyasnp# 假设有多个样本的质量数据（模拟数据）data_matrix=np.random.rand(5,100)*10+25# 5个样本，每个100个碱基位置sample_names=["Sample_A",'Sample_B", "Sample_C", "Sample_D", "Sample_E"]sns.set(style="whitegrid")plt.figure(figsize=(12,6))heatmap=sns.heatmap9data_matrix,xticklabels=range(1,101),yticklabels=sample_names,cmap="YlOrRd",annot=False)plt.title("Per-base Quality Heatmap Across Samples")plt.xlabel("Position (base)")plt.ylabel("Sample")plt.tight_layout()plt.savefig("quality_heatmap.png",dpi=300)

这张图能清晰展示哪些样本存在质量问题（如末端急剧下降），辅助决策是否重新测序或调整建库流程。

✅ 六、总结与最佳实践建议

💡 实践建议：

• 所有QC步骤应记录日志，推荐使用logging模块。
• • 清洗后务必再次运行 FastQC 验证效果。
• • 对于大规模项目，可用 Snakemake 封装上述流程形成可复现工作流。

通过这套基于 Python 的高质量控制流程，你不仅能显著提升数据可靠性，还能建立起标准化、可视化的 QC 工作模板，真正实现从“跑完数据”到“读懂数据”的跨越！