跳转至

Python

按照图片高度拆分为多个图 DeepSeek(元宝) 按需求写代码,零修改

  • 需求:写个python 脚本,读取图片,按照图片高度拆分为多个图,比如 输入5 时,按高度截取为5个子图
Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import argparse
import os
from PIL import Image

def validate_arguments(input_path, num_split):
    if not os.path.isfile(input_path):
        raise FileNotFoundError(f"输入文件 {input_path} 不存在")
    if num_split <= 0:
        raise ValueError("拆分数量必须大于0")

def split_image(input_path, output_dir, num_split):
    # 验证参数有效性
    validate_arguments(input_path, num_split)

    # 打开原始图片
    img = Image.open(input_path)
    width, height = img.size

    # 检查拆分数量是否超过图片高度
    if num_split > height:
        raise ValueError(f"拆分数量({num_split})超过图片高度({height})")

    # 计算每个子图的高度分配
    split_height = height // num_split
    remainder = height % num_split

    # 生成高度列表(前remainder个子图高度+1)
    heights = [split_height + 1 if i < remainder else split_height 
               for i in range(num_split)]

    # 准备输出目录和文件名
    os.makedirs(output_dir, exist_ok=True)
    filename = os.path.basename(input_path)
    name, ext = os.path.splitext(filename)

    # 拆分并保存图片
    current_height = 0
    for i, h in enumerate(heights):
        # 计算裁剪区域
        box = (0, current_height, width, current_height + h)
        cropped_img = img.crop(box)

        # 生成输出路径
        output_path = os.path.join(output_dir, f"{name}_{i}{ext}")
        cropped_img.save(output_path)

        current_height += h

if __name__ == "__main__":
    # 配置命令行参数解析
    parser = argparse.ArgumentParser(description="将图片按高度拆分为多个子图")
    parser.add_argument("-i", "--input",  required=True, help="输入图片路径")
    parser.add_argument("-o", "--output", required=True, help="输出目录")
    parser.add_argument("-n", "--num_split", type=int, required=True, help="拆分数量")

    # 解析参数并执行拆分
    args = parser.parse_args()
    split_image(args.input, args.output, args.num_split)

运行

Text Only
1
python aa.py -i 1.jpg -o o -n 5

1. 原文参考

● 利用贝叶斯定理判断SNV 一个位点所有可能的genotype如下图: ["AA", "CC", "TT", "GG", "AC", "AT", "AG", "CT", "CG", "TG"]

利用贝叶斯定理判断SNV的方法就是根据测序数据(D)计算每一种genotype(G)的概率,即Pr(G|D)。

也就是计算所有十种可能性的概率, 那么,概率最高的就是该位点的genotype。

根据贝叶斯定理计算Pr(G|D)公式如下:

\[Pr(G | D)=Pr(G)*Pr(D | G)/Pr(D)\]

我们的目的是比较不同genotype的概率,因此计算Pr(G)*Pr(D|G)/Pr(D)即可。

\[Pr(G|D)\propto Pr(G)*Pr(D|G)\]

现在有两个问题:

  1. 如何计算Pr(D|G)?
  2. 如何计算P(G)?

1.如何计算Pr(D|G)? 下面通过一个例子来说明计算Pr(D|G)的方法。data如下:

假设G为A1A2,如下:

Pr(D|G)的计算公式如下:

\[Pr(\mathcal{D}|G)=\prod_{i=1..n}Pr(b_{i}|\{A_{1},A_{2}\})=\prod_{i=1..n}\frac{Pr(b_{i}|A_{1})+Pr(b_{i}|A_{2})}{2}\]
\[b_{j}\in\{\mathrm{A},\mathrm{G},\mathrm{A}\}\]

P(b|A)的概率如下, 其中e为对应碱基的质量值:

\[\text{当 }b_i=A_j\quad Pr(b_i|A_j)=1-e_i\]
\[\text{当 }b_i\neq A_j\quad Pr(b_i|A_j)=e_i/3\]

当Genotype=AG时,套用上面的公式,首先分别计算不同碱基对应的结果:

\[Pr(b_1=\mathtt A|\mathtt A\mathtt G)=\quad\frac{1}{2}\left((1-10^{-2})+\frac{10^{-2}}{3}\right)=\quad0.49667\]
\[Pr(b_2=\mathtt G|\mathtt A\mathtt G)=\quad\frac{1}{2}\left(\frac{10^{-1}}{3}+(1-10^{-1})\right)=\quad0.466667\]
\[Pr(b_{3}=\mathtt A|\mathtt A\mathtt G)=\quad\frac{1}{2}\left((1-10^{-5})+\frac{10^{-5}}{3}\right)=\quad0.499997\]

最终的结果如下:

\[Pr({D}|{AG})=Pr(b_1={A}|{AG}) * Pr(b_2={G}|{AG}) * Pr(b_3={A}|{AG}) = 0.49667 * 0.466667 * 0.499997=0.115888\]
  1. 如何计算P(G)?

计算P(G), 就是计算10种可能的genotype的概率, 即Pr(AA),Pr(TT),Pr(CC)等等。

假定参考碱基为G,如果杂合突变率为0.001, 纯合突变率为0.0005。genotype的概率如下: 有了Pr(D|AG)和Pr(G)的概率, Pr(AG|D)就可以计算出来了, 如下:

\[Pr({AG}|{D})= Pr({D}|{AG})Pr({AG})=0.115888*0.001=0.000116\]

同理可以计算出其他genotype的概率, 如下:

因此Pr(AG)是10种genotype的概率最高的, 该位点的genotype为AG。

实际SNV calling的方法,可能还有很多细节, 或者其他方法,我这里分享的东西也可能存在不少错误,如果你发现了, 欢迎留言。

2. 代码实现

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import math
import itertools as it

# 计算基因型的先验概率
def calc_prior_prob(ref_base, g1, g2):
    # 假设已知的参数
    # REF_BASE = 'G'  # 参考碱基
    HET_RATE = 0.001  # 杂合突变率
    HOM_RATE = 0.0005  # 纯合突变率
    if g1 == g2 == ref_base:
        return 1.0 - HET_RATE - HOM_RATE
    elif g1 != g2:
        return HET_RATE
    else:
        return HOM_RATE

# 根据PHRED质量值计算错误率
def get_error_rate(qual):
    # 将PHRED分数转换为概率
    error_rate = 10 ** (-qual / 10.0)
    return error_rate

# 计算读取数据的条件概率
def calc_likelihood2(data, g1, g2, base_qualities):
    likelihood = 1.0
    for base, qual in zip(data, base_qualities):
        p_err = get_error_rate(qual)
        p_A1 = get_genotype_rate(base, p_err, g1)
        p_A2 = get_genotype_rate(base, p_err, g2)
        likelihood *= (p_A1 + p_A2) / 2
    return likelihood


def get_genotype_rate(base, p_err, genotype):
    if base == genotype:
        return 1.0 - p_err
    else:
        return p_err / 3

# 计算后验概率
def calc_posterior_prob(ref_base, data, base_qualities):
    probs = {}
    genotype_list = get_genotype(ref_base, data)
    for g1,g2 in genotype_list:
        prior = calc_prior_prob(ref_base, g1, g2)
        likelihood = calc_likelihood2(data, g1, g2, base_qualities)
        probs[(g1, g2)] = prior * likelihood
    return probs


def get_genotype(ref_base, data):
    res = []
    tmp = list(set(data))
    if ref_base not in tmp:
        tmp.append(ref_base)
    for e in it.combinations_with_replacement(tmp, 2):
        res.append(e)
    return res


# 使用示例
REF_BASE = 'G'  # 参考碱基
data = ['A', 'G', 'A']
base_qualities = [20, 10, 50]
posterior_probs = calc_posterior_prob(REF_BASE, data, base_qualities)

# 输出每种基因型的后验概率
for genotype, prob in sorted(posterior_probs.items(), key=lambda x: x[1], reverse=True):
    print(f"Genotype {genotype}: {prob:.6f}")