DNA/RNA序列分析

# 使用 BioPython 读取序列
from Bio import SeqIO

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

# 统计序列数量
print(f"序列数量: {len(sequences)}")

# 统计长度分布
lengths = [len(seq) for seq in sequences]
print(f"平均长度: {sum(lengths) / len(lengths):.2f}")
print(f"最短: {min(lengths)}, 最长: {max(lengths)}")

# 计算GC含量
from Bio.SeqUtils import gc_fraction
for seq in sequences[:5]:  # 展示前5个
    gc = gc_fraction(seq.seq) * 100
    print(f"{seq.id}: GC含量 = {gc:.2f}%")

from Bio.Blast import NCBIWWW

# 在线BLAST比对（需要联网）
sequence = "ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"
result = NCBIWWW.qblast("blastn", "nt", sequence)

# 解析比对结果
from Bio.Blast import NCBIXML
blast_records = NCBIXML.parse(result)

for record in blast_records:
    for alignment in record.alignments:
        for hsp in alignment.hsps:
            print(f"匹配: {alignment.title}")
            print(f"相似度: {hsp.identities / hsp.align_length * 100:.2f}%")
            print(f"E值: {hsp.expect}")

from Bio import Seq
from Bio.SeqUtils import seq1

# 创建DNA序列对象
dna_seq = Seq("ATGGCCATTGTAATGGGCCGCCTGAAAGGGTGCCCGATAG")

# 翻译成蛋白质
protein = dna_seq.translate()
print(f"蛋白质序列: {protein}")

# 查找所有ORF（开放阅读框）
def find_orfs(sequence, min_len=100):
    orfs = []
    for frame, strand in [(1, 1), (2, 1), (3, 1),
                              (1, -1), (2, -1), (3, -1)]:  # 6个读码框
        if strand == 1:
            seq = sequence[frame-1:]
        else:
            seq = sequence.reverse_complement()[frame-1:]

        for i in range(0, len(seq) - 2, 3):
            codon = seq[i:i+3]
            if len(codon) == 3 and codon == "ATG":  # 起始密码子
                for j in range(i+3, len(seq) - 2, 3):
                    if seq[j:j+3] in ["TAA", "TAG", "TGA"]:  # 终止密码子
                        orf_seq = seq[i:j+3]
                        if len(orf_seq) >= min_len:
                            orfs.append(orf_seq.translate())
                        break
    return orfs

orfs = find_orfs(dna_seq)
print(f"找到 {len(orfs)} 个ORF")

import pysam

# 打开BAM文件
bamfile = pysam.AlignmentFile("aligned.bam", "rb")

# 统计比对率
total = 0
mapped = 0
for read in bamfile.fetch():
    total += 1
    if not read.is_unmapped:
        mapped += 1

mapping_rate = mapped / total * 100
print(f"比对率: {mapping_rate:.2f}%")

# 提取特定区域的reads
chr_name = "chr1"
start = 1000000
end = 1005000
reads = bamfile.count_coverage(chr_name, start, end)
print(f"区域 {chr_name}:{start}-{end} 的覆盖度: {reads}")

# 计算平均覆盖度
avg_coverage = sum(reads[0]) / len(reads[0])
print(f"平均覆盖度: {avg_coverage:.2f}x")

bamfile.close()

import matplotlib.pyplot as plt
import pysam

# 1. 覆盖度图
bamfile = pysam.AlignmentFile("aligned.bam", "rb")
chr_name = "chr1"
start, end = 1000000, 1005000

coverage = []
for pos in range(start, end, 100):
    cov = bamfile.count(chr_name, pos, pos+100)
    coverage.append(cov)

plt.figure(figsize=(12, 4))
plt.plot(range(start, end, 100), coverage)
plt.xlabel("Position")
plt.ylabel("Coverage")
plt.title(f"Coverage of {chr_name}:{start}-{end}")
plt.savefig("coverage.png", dpi=300)

# 2. GC含量分布图
from Bio.SeqUtils import gc_fraction
gc_contents = [gc_fraction(seq.seq) for seq in sequences]

plt.figure(figsize=(8, 6))
plt.hist(gc_contents, bins=20, edgecolor='black')
plt.xlabel('GC Content')
plt.ylabel('Frequency')
plt.title('GC Content Distribution')
plt.savefig("gc_distribution.png", dpi=300)