合成生物学与生物计算安全取证深度分析

#!/bin/bash
STORAGE_FASTA=$1
REFERENCE_FASTA=$2
OUTPUT_DIR=$3

mkdir -p "$OUTPUT_DIR"

blastn -query "$STORAGE_FASTA" -subject "$REFERENCE_FASTA" \
    -outfmt 6 -out "$OUTPUT_DIR/blast_results.tsv" \
    -evalue 1e-10 -word_size 11 -num_threads 8

python3 << 'PYEOF'
import sys
from Bio import SeqIO
from Bio.Blast import NCBIXML
import subprocess

def calculate_gc_content(seq):
    seq_str = str(seq).upper()
    gc = sum(1 for b in seq_str if b in ('G', 'C'))
    return gc / len(seq_str) if len(seq_str) > 0 else 0

def detect_homopolymer_runs(seq, max_run=4):
    seq_str = str(seq).upper()
    runs = []
    current_base = seq_str[0]
    run_len = 1
    for i in range(1, len(seq_str)):
        if seq_str[i] == current_base:
            run_len += 1
        else:
            if run_len >= max_run:
                runs.append((current_base, run_len, i - run_len))
            current_base = seq_str[i]
            run_len = 1
    return runs

def analyze_encoding_pattern(fasta_file):
    records = list(SeqIO.parse(fasta_file, "fasta"))
    gc_contents = [calculate_gc_content(r.seq) for r in records]
    mean_gc = sum(gc_contents) / len(gc_contents)
    std_gc = (sum((gc - mean_gc)**2 for gc in gc_contents) / len(gc_contents))**0.5
    return {
        "total_sequences": len(records),
        "mean_gc": mean_gc,
        "std_gc": std_gc,
        "gc_skewed": abs(mean_gc - 0.5) > 0.1,
        "records": records
    }

def scan_for_anomalies(fasta_file, output_file):
    analysis = analyze_encoding_pattern(fasta_file)
    anomalies = []
    for record in analysis["records"]:
        seq_str = str(record.seq).upper()
        if 'N' in seq_str:
            anomalies.append({"id": record.id, "type": "ambiguous_base", "count": seq_str.count('N')})
        runs = detect_homopolymer_runs(record.seq)
        if runs:
            anomalies.append({"id": record.id, "type": "homopolymer_run", "runs": runs})
        if seq_str not in set("ATGC"):
            bad_chars = set(seq_str) - set("ATGC")
            anomalies.append({"id": record.id, "type": "invalid_base", "chars": list(bad_chars)})
    with open(output_file, 'w') as f:
        f.write(f"Total sequences: {analysis['total_sequences']}\n")
        f.write(f"GC content: mean={analysis['mean_gc']:.4f}, std={analysis['std_gc']:.4f}\n")
        f.write(f"GC skewed: {analysis['gc_skewed']}\n")
        f.write(f"Anomalies found: {len(anomalies)}\n\n")
        for a in anomalies:
            f.write(f"[{a['type']}] {a['id']}: {a}\n")
    return anomalies

scan_for_anomalies(sys.argv[1], sys.argv[2] if len(sys.argv) > 2 else "anomalies_report.txt")
PYEOF

echo "[DONE] DNA存储完整性验证完成，报告已输出至 $OUTPUT_DIR"

#!/usr/bin/env python3
import sys
import gzip
import hashlib
import statistics
from collections import Counter, defaultdict

def is_gzipped(filepath):
    with open(filepath, 'rb') as f:
        return f.read(2) == b'\x1f\x8b'

def read_fastq(filepath):
    opener = gzip.open if is_gzipped(filepath) else open
    with opener(filepath, 'rt') as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            f.readline()
            qual = f.readline().strip()
            yield header, seq, qual

def analyze_fastq(filepath):
    stats = {
        "total_reads": 0,
        "read_lengths": [],
        "quality_scores": [],
        "base_freqs": Counter(),
        "gc_contents": [],
        "anomalies": [],
        "duplicate_headers": 0,
        "header_set": set(),
        "length_variants": defaultdict(int),
        "quality_distribution": defaultdict(int),
    }
    phred_offset = 33
    for header, seq, qual in read_fastq(filepath):
        stats["total_reads"] += 1
        if stats["total_reads"] % 100000 == 0:
            print(f"  Processed {stats['total_reads']} reads...")
        if header in stats["header_set"]:
            stats["duplicate_headers"] += 1
        stats["header_set"].add(header)
        stats["read_lengths"].append(len(seq))
        stats["length_variants"][len(seq)] += 1
        gc = sum(1 for b in seq.upper() if b in ('G', 'C')) / len(seq) if seq else 0
        stats["gc_contents"].append(gc)
        for b in seq.upper():
            stats["base_freqs"][b] += 1
        for q in qual:
            qval = ord(q) - phred_offset
            stats["quality_scores"].append(qval)
            stats["quality_distribution"][qval] += 1
    return stats

def detect_anomalies(stats):
    anomalies = []
    if stats["total_reads"] == 0:
        anomalies.append({"severity": "CRITICAL", "type": "empty_file", "detail": "FASTQ file contains no reads"})
        return anomalies
    mean_len = statistics.mean(stats["read_lengths"])
    std_len = statistics.stdev(stats["read_lengths"]) if len(stats["read_lengths"]) > 1 else 0
    if std_len > mean_len * 0.5:
        anomalies.append({"severity": "HIGH", "type": "length_variance", "detail": f"std/mean={std_len/mean_len:.3f}"})
    mean_gc = statistics.mean(stats["gc_contents"])
    if abs(mean_gc - 0.5) > 0.15:
        anomalies.append({"severity": "MEDIUM", "type": "gc_bias", "detail": f"mean_gc={mean_gc:.4f}"})
    total_bases = sum(stats["base_freqs"].values())
    for base in ('A', 'T', 'G', 'C'):
        freq = stats["base_freqs"].get(base, 0) / total_bases if total_bases > 0 else 0
        if freq > 0.40:
            anomalies.append({"severity": "HIGH", "type": "base_frequency_anomaly", "detail": f"{base}={freq:.4f}"})
    if stats["duplicate_headers"] > 0:
        dup_rate = stats["duplicate_headers"] / stats["total_reads"]
        if dup_rate > 0.01:
            anomalies.append({"severity": "HIGH", "type": "duplicate_read_ids", "detail": f"rate={dup_rate:.4f}"})
    mean_qual = statistics.mean(stats["quality_scores"]) if stats["quality_scores"] else 0
    if mean_qual < 10:
        anomalies.append({"severity": "HIGH", "type": "low_quality", "detail": f"mean_phred={mean_qual:.2f}"})
    if 'N' in stats["base_freqs"]:
        n_rate = stats["base_freqs"]['N'] / total_bases
        if n_rate > 0.01:
            anomalies.append({"severity": "MEDIUM", "type": "high_n_content", "detail": f"N_rate={n_rate:.6f}"})
    n_lengths = len(stats["length_variants"])
    if n_lengths > 3:
        sorted_lens = sorted(stats["length_variants"].items(), key=lambda x: x[1], reverse=True)
        anomalies.append({"severity": "MEDIUM", "type": "multi_modal_length", "detail": f"variants={n_lengths}, top={sorted_lens[:3]}"})
    return anomalies

def generate_report(stats, anomalies, output_file):
    with open(output_file, 'w') as f:
        f.write("=" * 60 + "\n")
        f.write("FASTQ Anomaly Detection Report\n")
        f.write("=" * 60 + "\n\n")
        f.write(f"Total reads: {stats['total_reads']}\n")
        if stats['read_lengths']:
            f.write(f"Read length: mean={statistics.mean(stats['read_lengths']):.1f}, "
                    f"std={statistics.stdev(stats['read_lengths']):.2f}, "
                    f"min={min(stats['read_lengths'])}, max={max(stats['read_lengths'])}\n")
        if stats['gc_contents']:
            f.write(f"GC content: mean={statistics.mean(stats['gc_contents']):.4f}\n")
        if stats['quality_scores']:
            f.write(f"Quality: mean_phred={statistics.mean(stats['quality_scores']):.2f}\n")
        f.write(f"Duplicate headers: {stats['duplicate_headers']}\n\n")
        f.write("ANOMALY DETECTION RESULTS:\n")
        f.write("-" * 40 + "\n")
        if not anomalies:
            f.write("  No anomalies detected.\n")
        for a in anomalies:
            f.write(f"  [{a['severity']}] {a['type']}: {a['detail']}\n")
        f.write("\n")
        f.write("BASE FREQUENCY DISTRIBUTION:\n")
        total = sum(stats["base_freqs"].values())
        for base in ('A', 'T', 'G', 'C', 'N'):
            count = stats["base_freqs"].get(base, 0)
            f.write(f"  {base}: {count/total*100:.2f}%\n" if total > 0 else f"  {base}: N/A\n")

if __name__ == "__main__":
    if len(sys.argv) < 2:
        print(f"Usage: {sys.argv[0]} <fastq_file> [output_report]")
        sys.exit(1)
    input_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) > 2 else "fastq_anomaly_report.txt"
    print(f"[*] Analyzing {input_file}...")
    stats = analyze_fastq(input_file)
    anomalies = detect_anomalies(stats)
    generate_report(stats, anomalies, output_file)
    print(f"[*] Report generated: {output_file}")
    print(f"[*] Anomalies found: {len(anomalies)}")
    for a in anomalies:
        print(f"    [{a['severity']}] {a['type']}")

#!/bin/bash
LIMS_LOG_DIR=$1
BACKUP_LOG_DIR=$2
REPORT_FILE=$3

echo "=== LIMS Audit Log Integrity Check ===" > "$REPORT_FILE"
echo "Timestamp: $(date -u +%Y-%m-%dT%H:%M:%SZ)" >> "$REPORT_FILE"
echo "" >> "$REPORT_FILE"

TOTAL_LOGS=0
MODIFIED_LOGS=0
MISSING_LOGS=0

for log_file in "$LIMS_LOG_DIR"/*.log; do
    [ -f "$log_file" ] || continue
    TOTAL_LOGS=$((TOTAL_LOGS + 1))
    BASENAME=$(basename "$log_file")
    BACKUP_FILE="$BACKUP_LOG_DIR/$BASENAME"
    CURRENT_MD5=$(md5sum "$log_file" | awk '{print $1}')
    if [ ! -f "$BACKUP_FILE" ]; then
        echo "[MISSING] $BASENAME - no backup found" >> "$REPORT_FILE"
        MISSING_LOGS=$((MISSING_LOGS + 1))
        continue
    fi
    BACKUP_MD5=$(md5sum "$BACKUP_FILE" | awk '{print $1}')
    if [ "$CURRENT_MD5" != "$BACKUP_MD5" ]; then
        echo "[MODIFIED] $BASENAME" >> "$REPORT_FILE"
        echo "  Current MD5:  $CURRENT_MD5" >> "$REPORT_FILE"
        echo "  Backup MD5:   $BACKUP_MD5" >> "$REPORT_FILE"
        diff <(stat -c "%Y %s" "$log_file") <(stat -c "%Y %s" "$BACKUP_FILE") >> "$REPORT_FILE" 2>/dev/null
        MODIFIED_LOGS=$((MODIFIED_LOGS + 1))
    fi
done

TIMESTAMP_GAPS=0
PREV_TIME=0
for log_file in "$LIMS_LOG_DIR"/*.log; do
    [ -f "$log_file" ] || continue
    MTIME=$(stat -c "%Y" "$log_file" 2>/dev/null || stat -f "%m" "$log_file")
    if [ "$PREV_TIME" -gt 0 ] 2>/dev/null; then
        DIFF=$((MTIME - PREV_TIME))
        if [ "$DIFF" -lt 0 ] 2>/dev/null; then
            echo "[TIMESTAMP_ANOMALY] $log_file - time went backward by ${DIFF#-}s" >> "$REPORT_FILE"
            TIMESTAMP_GAPS=$((TIMESTAMP_GAPS + 1))
        elif [ "$DIFF" -gt 86400 ] 2>/dev/null; then
            echo "[TIMESTAMP_GAP] $log_file - gap of ${DIFF}s from previous log" >> "$REPORT_FILE"
            TIMESTAMP_GAPS=$((TIMESTAMP_GAPS + 1))
        fi
    fi
    PREV_TIME=$MTIME
done

echo "" >> "$REPORT_FILE"
echo "=== SUMMARY ===" >> "$REPORT_FILE"
echo "Total logs checked:     $TOTAL_LOGS" >> "$REPORT_FILE"
echo "Modified logs:          $MODIFIED_LOGS" >> "$REPORT_FILE"
echo "Missing backup logs:    $MISSING_LOGS" >> "$REPORT_FILE"
echo "Timestamp anomalies:    $TIMESTAMP_GAPS" >> "$REPORT_FILE"