超长片段测序 UltraLong¶

# 1. basecalling（碱基识别）—— 用 Dorado 模型
dorado basecaller sup input_pod5/ > calls.bam  # sup = 超高精度模型

# 2. 质控统计
NanoPlot --bam calls.bam -o nanoplot_out  # 查看读长分布和质量

# 3. 过滤超长读段（只保留 >100 kb 的）
samtools view -h calls.bam | \
  awk 'length($10) > 100000 || $1 ~ /^@/' | \  # 筛选序列长度>100kb
  samtools view -bS - > ultralong.bam           # 输出为BAM

# 4. 比对到参考基因组
minimap2 -ax map-ont ref.fa ultralong.bam | \  # map-ont 模式适配ONT数据
  samtools sort -o aligned.bam                  # 排序输出

# 5. 结构变异检测
sniffles --input aligned.bam --vcf sv.vcf \    # Sniffles检测SV
  --minsvlen 50                                 # 最小SV长度50bp

# 统计 UltraLong 读段长度分布（Python + pysam）
import pysam                          # BAM文件处理库
import numpy as np                    # 数值计算

bam = pysam.AlignmentFile("calls.bam", "rb")  # 打开BAM文件
lengths = []                          # 存储读段长度

for read in bam:                      # 遍历每条读段
    if not read.is_secondary:         # 跳过次要比对
        lengths.append(read.query_length)  # 记录读段长度

lengths = np.array(lengths)           # 转为numpy数组
print(f"总读段数: {len(lengths)}")     # 打印总数
print(f"N50: {calculate_n50(lengths)}")  # 打印N50
print(f">100kb读段数: {np.sum(lengths > 100000)}")  # 超长读段计数
print(f"最长读段: {np.max(lengths)/1e6:.2f} Mb")     # 最长的那条

def calculate_n50(lengths):           # 计算N50的函数
    sorted_lens = np.sort(lengths)[::-1]  # 从大到小排序
    cumsum = np.cumsum(sorted_lens)   # 累积和
    total = cumsum[-1]                # 总碱基数
    return sorted_lens[cumsum >= total/2][0]  # 累积到一半时的长度