求助：基于Python实现随机人类基因组序列生成及后续分析

Alright, let's tackle your genomic sequence workflow step by step. I'll walk you through generating the random sequences, searching for disease-linked patterns, and comparing the stats to real data—all with Python code you can tweak for your needs.

First, we need an efficient way to generate millions of ACGT sequences and save them to a text file. Using random.choices is fast for sequence generation, and writing in batches avoids slow逐行IO operations when dealing with large datasets.

import random

BASES = ['A', 'C', 'T', 'G']
SEQUENCE_LENGTH = 100  # Adjust this to match the length of your real genomic sequences
TOTAL_SEQUENCES = 1_000_000
OUTPUT_FILE = 'random_genomes.txt'

def generate_random_sequence(length):
    # Generate a single random sequence of specified length
    return ''.join(random.choices(BASES, k=length))

# Write sequences to file in batches to save memory and speed up IO
with open(OUTPUT_FILE, 'w') as f:
    batch_size = 10_000  # Tune this based on your system's memory
    batches = TOTAL_SEQUENCES // batch_size
    remaining = TOTAL_SEQUENCES % batch_size

    for _ in range(batches):
        # Generate a batch of sequences and write all at once
        batch = [generate_random_sequence(SEQUENCE_LENGTH) + '\n' for _ in range(batch_size)]
        f.writelines(batch)
    
    # Handle any leftover sequences
    if remaining > 0:
        final_batch = [generate_random_sequence(SEQUENCE_LENGTH) + '\n' for _ in range(remaining)]
        f.writelines(final_batch)

Next, we'll build a script to scan the generated sequences for your target disease-linked patterns. We'll track two key metrics: total occurrences of each pattern, and how many unique sequences contain each pattern.

DISEASE_PATTERNS = ['ATCGGTT', 'GCTTAAG', 'TTGACCA']  # Replace with your actual patterns
INPUT_FILE = 'random_genomes.txt'

def count_pattern_matches(file_path, patterns):
    # Initialize counters for total occurrences and sequence hits
    total_occurrences = {p: 0 for p in patterns}
    sequences_with_pattern = {p: 0 for p in patterns}

    with open(file_path, 'r') as f:
        for line in f:
            seq = line.strip()
            # Track if this sequence contains each pattern (to avoid double-counting sequences)
            seq_has_pattern = {p: False for p in patterns}
            
            for pattern in patterns:
                if pattern in seq:
                    # Count all instances of the pattern in the sequence
                    total_occurrences[pattern] += seq.count(pattern)
                    seq_has_pattern[pattern] = True
            
            # Update sequence hit counters
            for pattern, found in seq_has_pattern.items():
                if found:
                    sequences_with_pattern[pattern] += 1

    return total_occurrences, sequences_with_pattern

# Run the search on random sequences
random_total, random_seq_hits = count_pattern_matches(INPUT_FILE, DISEASE_PATTERNS)

print("=== Random Sequence Pattern Stats ===")
print("Total occurrences per pattern:")
for p, cnt in random_total.items():
    print(f"- {p}: {cnt}")

print("\nNumber of sequences containing each pattern:")
for p, cnt in random_seq_hits.items():
    print(f"- {p}: {cnt}")

Now we'll repeat the pattern search on your real genomic data and compare the results. This will show if your disease patterns are significantly more common in real data than in random noise.

REAL_DATA_FILE = 'real_genomes.txt'  # Path to your real genomic sequences

# Get stats from real data
real_total, real_seq_hits = count_pattern_matches(REAL_DATA_FILE, DISEASE_PATTERNS)

# Compare and print results
print("\n=== Random vs Real Data Comparison ===")
for pattern in DISEASE_PATTERNS:
    rand_total = random_total.get(pattern, 0)
    real_total_val = real_total.get(pattern, 0)
    rand_seq = random_seq_hits.get(pattern, 0)
    real_seq_val = real_seq_hits.get(pattern, 0)

    # Calculate ratios (avoid division by zero if a pattern never appears in random data)
    occ_ratio = real_total_val / rand_total if rand_total > 0 else "N/A"
    seq_ratio = real_seq_val / rand_seq if rand_seq > 0 else "N/A"

    print(f"\nPattern: {pattern}")
    print(f"  Total Occurrences: Random={rand_total}, Real={real_total_val}, Real/Random Ratio={occ_ratio:.2f}" if isinstance(occ_ratio, float) else f"  Total Occurrences: Random={rand_total}, Real={real_total_val}, Real/Random Ratio={occ_ratio}")
    print(f"  Sequences with Pattern: Random={rand_seq}, Real={real_seq_val}, Real/Random Ratio={seq_ratio:.2f}" if isinstance(seq_ratio, float) else f"  Sequences with Pattern: Random={rand_seq}, Real={real_seq_val}, Real/Random Ratio={seq_ratio}")

• If you're working with FASTA-formatted real data, use the biopython library to parse sequences more easily.
• For extremely large pattern sets or very long sequences, consider using the re module (regular expressions) or specialized sequence search algorithms like KMP for faster matching.
• Use numpy to generate sequences in bulk if you need even more speed for 1M+ entries.

内容的提问来源于stack exchange，提问作者Philip Estes