### Bulk RNA-Seq 完整分析流程 Bulk RNA-Seq 是一种用于研究整个细胞群体基因表达的技术。以下是基于 Python 的生物信息学工具实现 bulk RNA-Seq 数据处理和分析的完整流程。 #### 1. 质量控制 (Quality Control) 质量控制是数据分析的第一步，目的是评估测序数据的质量并过滤低质量序列。可以使用 `FastQC` 工具来生成报告，并通过 `multiqc` 进行汇总展示。 ```python import subprocess def run_fastqc(input_file, output_dir): command = f"fastqc {input_file} --outdir={output_dir}" subprocess.run(command.split(), check=True) run_fastqc("sample.fastq", "results/fastqc") # 示例调用[^2] ``` 对于进一步的数据清理，可利用 `cutadapt` 或其他类似的工具去除接头污染和其他噪声。 --- #### 2. 序列比对 (Read Alignment) 将高质量的 reads 映射到参考基因组上通常采用 STAR 或 HISAT2 等工具完成。这些工具可以通过命令行运行，也可以集成到 Python 中。 ```python def align_reads_with_star(genome_index, input_files, output_prefix): star_command = [ "STAR", "--genomeDir", genome_index, "--readFilesIn"] + list(input_files) + [ "--outFileNamePrefix", output_prefix ] subprocess.run(star_command, check=True) align_reads_with_star( "/path/to/genome/index", ["reads_1.fq.gz", "reads_2.fq.gz"], "./aligned/" ) # 示例调用[^3] ``` --- #### 3. 基因计数 (Gene Counting) 经过比对后的 BAM 文件需要转换成基因级别的读取计数矩阵。推荐使用的工具有 HTSeq 和 featureCounts。 ```python from pysam import Samfile def count_genes(samfile_path, gtf_file): samfile = Samfile(samfile_path, 'rb') counts_dict = {} with open(gtf_file) as gtf: for line in gtf: fields = line.strip().split("\t") gene_id = fields[-1].split(";")[0].strip() if gene_id not in counts_dict: counts_dict[gene_id] = 0 for read in samfile.fetch(): if not read.is_unmapped: gene_name = get_gene_from_read(read) # 自定义函数获取基因名 counts_dict[gene_name] += 1 return counts_dict counts = count_genes("./aligned/Aligned.out.bam", "genes.gtf") # 示例调用[^4] ``` --- #### 4. 差异表达分析 (Differential Expression Analysis) 差异表达分析常用 DESeq2（R包），但在 Python 生态中有替代方案如 `scikit-learn` 配合统计模型或者直接调用 R 接口。 ```python import pandas as pd from rpy2.robjects.packages import importr from rpy2.robjects import DataFrame, FloatVector deseq = importr('DESeq2') # 构建输入数据框 count_data = pd.read_csv("gene_counts.csv", index_col=0) metadata = pd.DataFrame({ "condition": ["control", "treatment"], }, index=["S1", "S2"]) dds = deseq.DESeqDataSetFromMatrix(countData=DataFrame(count_data), colData=DataFrame(metadata), design=~FloatVector([1])) res = deseq.nbinomWaldTest(dds) result_table = deseq.results(res)[^5] pd.DataFrame(result_table).to_csv("diff_expr_results.csv") ``` --- #### 5. 可视化与功能富集分析 (Visualization & Enrichment Analysis) 最后一步是对结果进行可视化以及生物学意义解释。常用的库有 Matplotlib、Seaborn 和 GSEAPy。 ```python import matplotlib.pyplot as plt import seaborn as sns sns.set_theme(style="whitegrid") data = pd.read_csv("diff_expr_results.csv") plt.figure(figsize=(8, 6)) sns.scatterplot(x=data["log2FoldChange"], y=-np.log10(data["padj"]), hue=np.where(data["padj"] < 0.05, "Significant", "Not Significant")) plt.title("Volcano Plot of Differential Gene Expression") plt.show()[^6] ``` --- ### 总结 以上展示了如何构建一个完整的 bulk RNA-Seq 分析管道，涵盖了从原始 FASTQ 到最终差异表达结果的过程。每一步都依赖于特定的软件或算法，并可通过 Python 实现自动化脚本开发。