00-Python在生信中的应用

from Bio import SeqIO
import csv

# 读取FASTA文件
input_file = "sequences.fasta"
sequences = list(SeqIO.parse(input_file, "fasta"))

# 将FASTA转换为CSV
output_file = "sequences.csv"
with open(output_file, "w", newline='') as csvfile:
    csv_writer = csv.writer(csvfile)
    # 写入CSV文件的标题行
    csv_writer.writerow(["ID", "Sequence", "Length"])

    # 遍历序列并将ID、序列和长度写入CSV文件
    for seq_record in sequences:
        csv_writer.writerow([seq_record.id, str(seq_record.seq), len(seq_record)])

from Bio import SeqIO
from Bio import pairwise2
from Bio.SubsMat import MatrixInfo as matlist

# 两个待比对的DNA序列
seq1 = "ATGCTAGCTAGCTAGCTGAC"
seq2 = "ATGCTTGCTAGCTTGCTGAC"

# 使用BLOSUM62得分矩阵
matrix = matlist.blosum62

# 设置gap penalties
gap_open = -10
gap_extend = -0.5

# 计算全局pairwise序列比对
alignments = pairwise2.align.globalds(seq1, seq2, matrix, gap_open, gap_extend)

# 获取最佳比对结果
best_alignment = alignments[0]
aligned_seq1, aligned_seq2, score, begin, end = best_alignment

# 打印比对结果
print("序列1: {}\n序列2: {}\n得分: {}".format(aligned_seq1, aligned_seq2, score))

import gffutils

# 读取GFF3文件并创建数据库
input_file = "annotations.gff3"
db = gffutils.create_db(input_file, dbfn="annotations.db", force=True, keep_order=True)

# 查询所有基因
genes = db.features_of_type("gene")

# 打印基因ID及其坐标
for gene in genes:
    print("基因ID: {}, 染色体: {}, 起始坐标: {}, 结束坐标: {}".format(
        gene.id, gene.chrom, gene.start, gene.end))

# 查询某个基因的所有外显子
gene_id = "example_gene_id"
exons = db.children(gene_id, featuretype="exon")

# 打印外显子ID及其坐标
for exon in exons:
    print("外显子ID: {}, 染色体: {}, 起始坐标: {}, 结束坐标: {}".format(
        exon.id, exon.chrom, exon.start, exon.end))

import numpy as np
import matplotlib.pyplot as plt

# 示例数据
genes = ["Gene1", "Gene2", "Gene3", "Gene4", "Gene5"]
expression_values = [10, 15, 7, 12, 4]

# 创建柱状图
fig, ax = plt.subplots()
index = np.arange(len(genes))
bar_width = 0.5

rects = ax.bar(index, expression_values, bar_width, label="Expression")

# 设置图表标题和坐标轴标签
ax.set_xlabel("Genes")
ax.set_ylabel("Expression Values")
ax.set_title("Gene Expression Levels")
ax.set_xticks(index)
ax.set_xticklabels(genes)

# 添加数据标签
for rect in rects:
    height = rect.get_height()
    ax.annotate("{}".format(height),
                xy=(rect.get_x() + rect.get_width() / 2, height),
                xytext=(0, 3),  # 3 points vertical offset
                textcoords="offset points",
                ha="center", va="bottom")

# 显示图形
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans

# 示例数据：10个样本，每个样本5个基因的表达值
gene_expression_data = np.array([
    [10, 15, 7, 12, 4],
    [9, 14, 8, 11, 5],
    [11, 16, 6, 13, 3],
    [30, 35, 27, 32, 24],
    [29, 34, 28, 31, 25],
    [31, 36, 26, 33, 23],
    [50, 55, 47, 52, 44],
    [49, 54, 48, 51, 45],
    [51, 56, 46, 53, 43],
    [12, 17, 5, 14, 2]
])

# PCA降维
pca = PCA(n_components=2)
pca_result = pca.fit_transform(gene_expression_data)

# K-means聚类
kmeans = KMeans(n_clusters=3, random_state=0).fit(pca_result)

# 绘制PCA散点图，颜色表示聚类结果
plt.scatter(pca_result[:, 0], pca_result[:, 1], c=kmeans.labels_, cmap='viridis')

# 设置图表标题和坐标轴标签
plt.xlabel("Principal Component 1")
plt.ylabel("Principal Component 2")
plt.title("PCA and Clustering of Gene Expression Data")

# 显示图形
plt.show()

import networkx as nx
import matplotlib.pyplot as plt

# 创建一个空的有向图
G = nx.DiGraph()

# 添加节点（这里我们使用基因名称作为节点）
genes = ["Gene1", "Gene2", "Gene3", "Gene4", "Gene5"]
G.add_nodes_from(genes)

# 添加边（表示基因间的相互作用关系）
edges = [("Gene1", "Gene2"), ("Gene1", "Gene3"), ("Gene2", "Gene3"),
         ("Gene3", "Gene4"), ("Gene3", "Gene5"), ("Gene4", "Gene5")]
G.add_edges_from(edges)

# 计算每个节点的入度和出度
in_degrees = G.in_degree()
out_degrees = G.out_degree()

# 打印节点的入度和出度
for gene in genes:
    print("基因: {}, 入度: {}, 出度: {}".format(gene, in_degrees[gene], out_degrees[gene]))

# 绘制网络图
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True, node_color="skyblue", font_size=10, font_weight="bold")
plt.title("Gene Interaction Network")
plt.show()

from Bio.PDB import PDBParser

# 读取PDB文件
pdb_file = "example.pdb"
parser = PDBParser()
structure = parser.get_structure("example_protein", pdb_file)

# 遍历蛋白质结构中的链、残基和原子
for model in structure:
    for chain in model:
        print("链: ", chain.id)
        for residue in chain:
            print("  残基: ", residue.resname, residue.id[1])
            for atom in residue:
                print("    原子: ", atom.name, atom.coord)

pip install snakemake

# Snakefile
import os

# 定义输入文件夹，输出文件夹和样本
input_dir = "raw_data"
output_dir = "results"
samples = ["sample1", "sample2", "sample3"]

# 定义规则：FastQC质量控制
rule fastqc:
    input:
        os.path.join(input_dir, "{sample}.fastq")
    output:
        os.path.join(output_dir, "fastqc", "{sample}_fastqc.html")
    shell:
        "fastqc {input} --outdir={output_dir}/fastqc"

# 定义规则：比对到参考基因组
rule align:
    input:
        os.path.join(input_dir, "{sample}.fastq")
    output:
        os.path.join(output_dir, "aligned", "{sample}.sam")
    shell:
        "hisat2 -x path/to/reference_genome -U {input} -S {output}"

# 定义规则：计算表达量
rule quantification:
    input:
        os.path.join(output_dir, "aligned", "{sample}.sam")
    output:
        os.path.join(output_dir, "counts", "{sample}.txt")
    shell:
        "htseq-count -f sam {input} path/to/annotation.gtf > {output}"

# 定义工作流程的目标
rule all:
    input:
        expand(os.path.join(output_dir, "counts", "{sample}.txt"), sample=samples)

snakemake --cores 4