Genome Repository
基因组数据分析和可视化的Jupyter Notebook集合,包含多种生物信息学分析方法
技术栈
Genome Repository
项目简介
这是一个专门用于基因组数据分析和可视化的Jupyter Notebook集合。项目包含了多种生物信息学分析方法的实现,从基础的序列分析到复杂的基因组比较研究,为研究人员提供了完整的分析工具包。
技术框架
🧬 核心技术栈
- Python: 主要编程语言
- Jupyter Notebook: 交互式分析环境
- Pandas: 数据处理和分析
- BioPython: 生物序列处理
- Matplotlib/Seaborn: 数据可视化
- NumPy: 数值计算
- Scikit-learn: 机器学习算法
📊 分析流程
原始数据 → 数据清洗 → 序列分析 → 统计分析 → 可视化 → 结果解释
↓ ↓ ↓ ↓ ↓ ↓
FASTA DataFrame Alignment Statistics Plots Report
主要模块
1. 序列分析模块
import pandas as pd
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqUtils import GC, molecular_weight
class SequenceAnalyzer:
def __init__(self, fasta_file):
self.sequences = list(SeqIO.parse(fasta_file, "fasta"))
def basic_stats(self):
"""计算基本序列统计信息"""
stats = []
for record in self.sequences:
seq_stats = {
'id': record.id,
'length': len(record.seq),
'gc_content': GC(record.seq),
'molecular_weight': molecular_weight(record.seq, 'DNA'),
'a_count': record.seq.count('A'),
'c_count': record.seq.count('C'),
'g_count': record.seq.count('G'),
't_count': record.seq.count('T')
}
stats.append(seq_stats)
return pd.DataFrame(stats)
def find_orfs(self, min_length=100):
"""寻找开放阅读框"""
orfs = []
for record in self.sequences:
seq = record.seq
for frame in range(3):
for start in range(frame, len(seq), 3):
codon = seq[start:start+3]
if len(codon) == 3 and codon == 'ATG': # 起始密码子
for end in range(start+3, len(seq), 3):
stop_codon = seq[end:end+3]
if stop_codon in ['TAA', 'TAG', 'TGA']:
orf_length = end - start
if orf_length >= min_length:
orfs.append({
'sequence_id': record.id,
'frame': frame + 1,
'start': start + 1,
'end': end + 3,
'length': orf_length,
'sequence': seq[start:end+3]
})
break
return pd.DataFrame(orfs)
2. 基因组比较分析
import matplotlib.pyplot as plt
import seaborn as sns
from Bio.Blast import NCBIXML
from Bio.Blast.Applications import NcbiblastnCommandline
class GenomeComparator:
def __init__(self, genome1_file, genome2_file):
self.genome1 = list(SeqIO.parse(genome1_file, "fasta"))
self.genome2 = list(SeqIO.parse(genome2_file, "fasta"))
def calculate_similarity_matrix(self):
"""计算序列相似性矩阵"""
from Bio.Align import PairwiseAligner
aligner = PairwiseAligner()
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -2
aligner.extend_gap_score = -0.5
similarity_matrix = []
for seq1 in self.genome1:
row = []
for seq2 in self.genome2:
alignment = aligner.align(seq1.seq, seq2.seq)
score = alignment.score / max(len(seq1.seq), len(seq2.seq))
row.append(score)
similarity_matrix.append(row)
return np.array(similarity_matrix)
def plot_similarity_heatmap(self, similarity_matrix):
"""绘制相似性热图"""
plt.figure(figsize=(12, 8))
sns.heatmap(similarity_matrix,
annot=True,
cmap='viridis',
xticklabels=[seq.id for seq in self.genome2],
yticklabels=[seq.id for seq in self.genome1])
plt.title('Genome Similarity Matrix')
plt.xlabel('Genome 2 Sequences')
plt.ylabel('Genome 1 Sequences')
plt.tight_layout()
plt.show()
3. 变异分析模块
class VariantAnalyzer:
def __init__(self, vcf_file):
self.variants = self.parse_vcf(vcf_file)
def parse_vcf(self, vcf_file):
"""解析VCF文件"""
variants = []
with open(vcf_file, 'r') as f:
for line in f:
if not line.startswith('#'):
fields = line.strip().split('\t')
variant = {
'chrom': fields[0],
'pos': int(fields[1]),
'ref': fields[3],
'alt': fields[4],
'qual': float(fields[5]) if fields[5] != '.' else None,
'filter': fields[6],
'info': fields[7]
}
variants.append(variant)
return pd.DataFrame(variants)
def variant_distribution(self):
"""分析变异分布"""
# 按染色体统计变异数量
chrom_counts = self.variants['chrom'].value_counts()
# 按变异类型统计
variant_types = []
for _, row in self.variants.iterrows():
if len(row['ref']) == len(row['alt']) == 1:
variant_types.append('SNP')
elif len(row['ref']) > len(row['alt']):
variant_types.append('Deletion')
elif len(row['ref']) < len(row['alt']):
variant_types.append('Insertion')
else:
variant_types.append('Complex')
self.variants['type'] = variant_types
type_counts = self.variants['type'].value_counts()
return chrom_counts, type_counts
def plot_variant_distribution(self):
"""绘制变异分布图"""
chrom_counts, type_counts = self.variant_distribution()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# 染色体分布
chrom_counts.plot(kind='bar', ax=ax1)
ax1.set_title('Variants per Chromosome')
ax1.set_xlabel('Chromosome')
ax1.set_ylabel('Number of Variants')
ax1.tick_params(axis='x', rotation=45)
# 变异类型分布
type_counts.plot(kind='pie', ax=ax2, autopct='%1.1f%%')
ax2.set_title('Variant Type Distribution')
plt.tight_layout()
plt.show()
分析案例
案例1:人类基因组变异分析
# 加载数据
analyzer = VariantAnalyzer('human_variants.vcf')
# 基本统计
print(f"总变异数: {len(analyzer.variants)}")
print(f"涉及染色体: {analyzer.variants['chrom'].nunique()}")
# 质量分布
high_quality_variants = analyzer.variants[analyzer.variants['qual'] > 30]
print(f"高质量变异: {len(high_quality_variants)}")
# 可视化
analyzer.plot_variant_distribution()
案例2:物种间基因组比较
# 比较两个物种的基因组
comparator = GenomeComparator('species_a.fasta', 'species_b.fasta')
# 计算相似性
similarity_matrix = comparator.calculate_similarity_matrix()
# 可视化结果
comparator.plot_similarity_heatmap(similarity_matrix)
# 统计分析
mean_similarity = np.mean(similarity_matrix)
print(f"平均相似性: {mean_similarity:.3f}")
案例3:基因功能预测
class GeneFunctionPredictor:
def __init__(self, protein_sequences):
self.sequences = protein_sequences
def predict_domains(self):
"""预测蛋白质结构域"""
# 使用HMM模型预测结构域
domains = []
for seq in self.sequences:
# 简化的结构域预测逻辑
if 'GXXGXGK' in str(seq.seq):
domains.append('ATP_binding')
elif 'CXXC' in str(seq.seq):
domains.append('Zinc_finger')
else:
domains.append('Unknown')
return domains
def functional_classification(self):
"""功能分类"""
domains = self.predict_domains()
classification = pd.DataFrame({
'sequence_id': [seq.id for seq in self.sequences],
'predicted_domain': domains
})
return classification
数据可视化
基因组浏览器
def create_genome_browser(variants_df, genes_df):
"""创建简单的基因组浏览器"""
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 10), sharex=True)
# 绘制基因
for _, gene in genes_df.iterrows():
ax1.barh(0, gene['end'] - gene['start'],
left=gene['start'], height=0.5,
color='blue', alpha=0.7)
ax1.text(gene['start'] + (gene['end'] - gene['start'])/2, 0,
gene['name'], ha='center', va='center')
ax1.set_ylim(-0.5, 0.5)
ax1.set_title('Genes')
ax1.set_ylabel('Genes')
# 绘制变异
snps = variants_df[variants_df['type'] == 'SNP']
indels = variants_df[variants_df['type'].isin(['Insertion', 'Deletion'])]
ax2.scatter(snps['pos'], [1]*len(snps),
color='red', alpha=0.6, s=20, label='SNPs')
ax2.scatter(indels['pos'], [2]*len(indels),
color='green', alpha=0.6, s=20, label='InDels')
ax2.set_ylim(0.5, 2.5)
ax2.set_title('Variants')
ax2.set_ylabel('Variant Type')
ax2.set_xlabel('Genomic Position')
ax2.legend()
plt.tight_layout()
plt.show()
项目结构
genome-repo/
├── notebooks/
│ ├── 01_sequence_analysis.ipynb
│ ├── 02_variant_calling.ipynb
│ ├── 03_genome_comparison.ipynb
│ ├── 04_phylogenetic_analysis.ipynb
│ └── 05_functional_annotation.ipynb
├── data/
│ ├── sample_genomes/
│ ├── reference_data/
│ └── test_datasets/
├── scripts/
│ ├── data_preprocessing.py
│ ├── analysis_utils.py
│ └── visualization.py
├── results/
│ ├── figures/
│ └── reports/
└── requirements.txt
使用方法
环境设置
# 克隆仓库
git clone https://github.com/quanxquan/genome-repo.git
cd genome-repo
# 安装依赖
pip install -r requirements.txt
# 启动Jupyter
jupyter notebook
数据准备
# 数据格式要求
# FASTA格式的基因组序列
# VCF格式的变异数据
# GFF格式的基因注释
# 示例数据加载
from scripts.analysis_utils import load_genome_data
genome_data = load_genome_data('data/sample_genomes/sample.fasta')
variant_data = load_variant_data('data/variants/sample.vcf')
应用领域
1. 比较基因组学
- 物种间基因组比较
- 进化关系分析
- 基因家族研究
2. 医学基因组学
- 疾病相关变异分析
- 药物基因组学研究
- 个性化医疗
3. 农业基因组学
- 作物改良
- 抗性基因挖掘
- 品质性状分析
项目成果
📊 分析统计
- 处理基因组数量:50+
- 分析变异位点:100万+
- 生成报告:20+
🔬 科研贡献
- 支持研究项目:5个
- 发现新变异:30+
- 方法改进:3项
未来发展
技术改进
- 集成更多分析算法
- 优化大数据处理性能
- 增加交互式可视化
- 支持云端计算
功能扩展
- 单细胞基因组分析
- 表观基因组学分析
- 多组学数据整合
- 自动化报告生成
学习价值
这个项目帮助我掌握了:
-
生物信息学核心技能
- 序列分析算法
- 统计方法应用
- 数据可视化技巧
-
编程实践能力
- Python科学计算
- Jupyter Notebook使用
- 代码模块化设计
-
科研思维培养
- 问题分析方法
- 实验设计思路
- 结果解释能力
-
工具开发经验
- 用户需求分析
- 文档编写规范
- 开源项目管理
通过这个项目,我不仅提升了技术能力,更重要的是培养了将复杂生物学问题转化为计算问题的能力,为后续的科研工作奠定了坚实基础。