UVA - 1606 Amphiphilic Carbon Molecules(极角排序+扫描法+计算几何)

本文介绍了一道计算几何题目的解法,通过极角排序和扫描法优化算法的时间复杂度,最终实现了高效的解题策略。

题目链接:

https://uva.onlinejudge.org/index.php?option=com_onlinejudge&Itemid=8&page=show_problem&problem=4481


题目大意:

平面上有 n 个点, 每个点为白点或黑点。现在需要放置一条隔板,使得隔板一侧的黑色点数加上另一侧白色点数的和最大。


解题过程:

这是人生第一道 AC 的计算几何题!
由于之前一直没接触过这种题,紫书上的分析也看不懂,一言不合就让我极角排序、叉乘,根本没学过啊QAQ。
然后翻了很多博客,重学了一下极坐标,可算是弄明白了。


题目分析:

首先这题最简单的想法是枚举所有可能的挡板,这样需要枚举两个点,然后统计两端的点。枚举是 n×n, 统计也是 n,总共时间复杂度是 n^3。显然是超时的。

然后扫描法的思想是,按照一定顺序枚举,这样可以动态的维护一个量,省去了扫描时的 n。

所以这题正解是枚举每一个点当作基点,然后计算每个点相对于基点的极角,进行极角排序。然后按极角的大小进行旋转时的枚举。每旋转一次可以根据之前统计的点的个数进行动态的更新(这就是“维护”)。

这题有个可以优化的地方,如果一个点是黑点的话,可以把他映射到关于基点对称的地方去,视为白点。这样扫描的时候只需要统计白点就够了。然后扫描的时候只从扫描 180 度即可。

判断是否超过 180 度用叉乘,a×b = |a| * |b| * sin<a,b>,两个点叉乘用坐标计算是,x1*y2 - y1*x2 因为向量的模一定是正的,所以只要判断这个叉乘正负,就可以判断 sin 的正负了,即是否超过 180 度。


AC代码:

#include<bits/stdc++.h>
using namespace std;

const int MAX = 1005;

struct Node {
    int x, y;
    int color;
    double rad;
}raw_data[MAX], data[MAX];

bool operator < (const struct Node &a, const struct Node &b) {
    return a.rad < b.rad;
}

int n, ans;

//用叉乘判断角度是否大于180度,即sin值。
bool Turn(Node &a, Node &b) {
    return a.x * b.y >= a.y * b.x;
}

void solve() {
    if (n <= 3) {
        ans = n;
        return;
    }
    ans = 0;

    //枚举每一个点当多基点
    for (int i = 0; i < n; i++) {
        int k = 0;
        for (int j = 0; j < n; j++) {
            if (i == j)
                continue;

            //计算每个点关于基点的相对坐标。
            data[k].x = raw_data[j].x - raw_data[i].x;
            data[k].y = raw_data[j].y - raw_data[i].y;

            //如果是黑点的话,那么映射到关于基点对称的地方,这样扫描的时候只记录白点就够了。
            if (raw_data[j].color) {
                data[k].x *= -1;
                data[k].y *= -1;
            }
            //求极角
            data[k].rad = atan2(data[k].y, data[k].x);
            k++;
        }

        //极角排序
        sort(data, data + k);

        //L为与基点相连的点,R为扫描线。
        int L = 0, R = 0, cnt = 1;

        //枚举每一个点与基点相连当作挡板。
        while (L < k) {
            if (L == R) {
                R = (R + 1) % k;
                cnt++;
            }

            //扫描的时候只扫描一侧,保证扫描线与挡板的夹角不超过180度。
            while (L != R && Turn(data[L], data[R])) {
                cnt++;
                R = (R + 1) % k;
            }
            ans = max(ans, cnt);
            //挡板旋转。
            L++;
            //由于挡板旋转,与基点相连的点不在挡板上了,减去1。
            cnt--;
        }
    }
}

int main() {
    while (cin >> n && n) {
        for (int i = 0; i < n; i++) {
            cin >> raw_data[i].x >> raw_data[i].y >> raw_data[i].color;
        }
        solve();
        cout << ans << endl;
    }
}
<think>我们使用PseAAC(伪氨基酸组成)方法对蛋白质序列进行特征提取。PseAAC不仅包含氨基酸组成信息,还包含序列顺序信息和物理化学性质。对于疏水性、亲水性、分子量、pK1、pK2和等电点等性质,我们可以计算每个氨基酸的这些性质,然后通过PseAAC方法构建特征向量。 步骤: 1. 准备7000个蛋白质序列(FASTA格式)。 2. 为每种物理化学性质(疏水性、亲水性、分子量、pK1、pK2、等电点)定义每个氨基酸的数值。通常这些数值有标准表。 3. 计算每个序列的伪氨基酸组成特征。PseAAC特征包括两部分:一部分是20个标准氨基酸的组成(20维),另一部分是序列相关因子(λ个相关因子,λ是用户定义的参数)。但这里我们更关注物理化学性质,因此我们可以使用基于物理化学性质的PseAAC变体(如Am-PseAAC,即Amphiphilic Pseudo Amino Acid Composition,它考虑疏水和亲水性质)或者自定义一组性质。 具体到每个性质(疏水性、亲水性、分子量、pK1、pK2、等电点),我们可以分别计算每个序列的PseAAC特征。但是,通常PseAAC一次使用一组性质(例如疏水和亲水一起用于Am-PseAAC)。我们可以分别对每个性质进行特征提取,也可以组合多个性质。 然而,传统的PseAAC一次只能使用一种物理化学性质(或者一组相关的性质,如疏水和亲水作为两亲性)。如果要同时考虑多个独立的物理化学性质,我们可以为每个性质单独计算PseAAC特征,然后将这些特征拼接起来。 另一种方法是使用扩展的PseAAC(如PseAAC-General),它允许同时使用多个性质。具体步骤: a. 选择一组物理化学性质(这里包括疏水性、亲水性、分子量、pK1、pK2、等电点)。 b. 为每个性质定义一个归一化的数值表(即每个氨基酸在该性质下的数值)。 c. 对于每个性质,计算序列的相关因子(自相关函数),然后组合成特征向量。 由于计算量较大,我们需要一个工具来帮助计算。常用的工具有: - `protr` (R包) - `iFeature` (Python) - 或者自己编写脚本 这里我们使用Python,因为更通用。我们可以使用`biopython`来读取序列,然后自己实现PseAAC计算,或者使用`iFeature`工具。 使用`iFeature`: `iFeature`是一个强大的Python工具,可以计算多种特征,包括PseAAC及其变体。 安装:`pip install iFeature` 但是,`iFeature`目前支持的PseAAC类型包括: - 标准的PseAAC(使用λ个序列顺序相关因子) - Am-PseAAC(两亲性伪氨基酸组成,使用疏水和亲水性质) - PseAAC-General(可以自定义多个物理化学性质) 因此,我们可以使用PseAAC-General来同时计算多个性质。 步骤: 1. 准备数据:将7000个蛋白质序列保存到一个FASTA文件中,例如`sequences.fasta`。 2. 编写Python脚本,使用`iFeature`计算PseAAC-General特征。 注意:PseAAC-General需要指定以下参数: - 使用的物理化学性质列表(例如我们选择:hydrophobicity, hydrophilicity, molecular_weight, pK1, pK2, pI) - λ:表示相关因子的阶数(通常取小于序列最小长度的一半,比如5,10,15等) - w:权重因子(默认0.05) 但是,`iFeature`的PseAAC-General默认支持的物理化学性质只有20多种(在`AAIndex`数据库中)。我们需要确认这些性质在AAIndex中的编号。例如: - 疏水性:AAIndex中有多个,常用的是HOPT810101(来自Hope等) - 亲水性:HOPT810101(注意:同一个编号可能代表不同性质,需要确认) - 分子量:通常使用标准分子量(在`iFeature`中可能已经内置) - pK1, pK2, pI:这些是解离常数和等电点,AAIndex中可能没有直接对应的,需要自定义。 如果AAIndex中没有我们需要的性质,我们可以自定义每个氨基酸的这些性质值。 在`iFeature`中,我们可以通过提供自定义的属性文件来使用自定义性质。自定义属性文件是一个文本文件,每行一个氨基酸字母,后面跟着该性质的数值(多个性质则多列)。 例如,我们可以创建一个名为`custom_properties.txt`的文件,包含以下内容(示例数据,实际数值需要从文献或数据库中获取): A 0.62 0.5 89.09 2.34 9.69 6.00 R -2.53 0.5 174.20 2.17 9.04 10.76 ... (其他氨基酸) 其中,每列分别代表:氨基酸字母、疏水性(自定义标度)、亲水性(自定义标度)、分子量、pK1、pK2、等电点。 注意:由于不同文献给出的数值可能有差异,需要确保使用一致的数值。 然后,在`iFeature`中,我们可以指定自定义属性文件,并选择使用哪些列(即哪些性质)。 具体步骤: 1. 准备自定义属性文件(如上所述)。 2. 使用`iFeature`的`PseAAC-General`描述器,并指定自定义属性文件路径,以及选择使用哪些列(列索引从0开始,比如我们想用第1列(疏水性)、第2列(亲水性)、第3列(分子量)等,则指定`cols=0,1,2,3,4,5`,因为索引从0开始,且第一列是氨基酸字母,后面6列是6个性质,但注意:我们自定义文件中,第一列是氨基酸字母,后面6列是6个性质值,所以实际性质值从第1列到第6列(索引1到6)?不对,在自定义文件中,第一列是氨基酸字母,然后第二列是第一个性质,第三列是第二个性质,以此类推。所以我们在指定列的时候,应该指定从第二列开始(索引1)到第七列(索引6)?不对,`iFeature`要求自定义属性文件中,第一列是氨基酸,后面每列是一个性质。在调用时,我们通过`cols`参数指定使用哪些性质列(列索引从0开始,但0是氨基酸列,所以性质列从1开始)。例如,我们想使用第2列(索引1)到第7列(索引6),则`cols=1,2,3,4,5,6`。 3. 设置λ(自相关因子阶数)和权重w。 但是,由于`iFeature`的文档中关于自定义属性文件的使用说明不够详细,我们可以参考其GitHub上的例子:https://github.com/Superzchen/iFeature/ 另外,我们也可以自己编写代码计算。但考虑到7000条序列,自己编写可能效率较低,使用现成工具更好。 如果使用`iFeature`,我们可以这样写脚本: ```python from iFeature import iFeature # 创建一个iFeature对象 ifr = iFeature.iFeature() # 设置输入文件(FASTA格式) input_file = "sequences.fasta" # 设置输出文件 output_file = "PseAAC_features.tsv" # 自定义属性文件路径 custom_prop_file = "custom_properties.txt" # 使用PseAAC-General,参数: # - lambda: 自相关阶数,例如5 # - w: 权重因子,例如0.05 # - custom: 自定义属性文件 # - cols: 自定义属性文件中使用的列(从0开始索引,但0是氨基酸列,所以性质列从1开始) # 注意:cols参数是一个字符串,用逗号分隔的列索引,例如我们使用第2列(索引1)到第7列(索引6),则"1,2,3,4,5,6" # 另外,我们还需要指定AAIndex类型为'custom' params = { "type": "PseAAC-General", "method": "auto", "lambda": 5, # 阶数,可以根据需要调整 "w": 0.05, # 权重 "custom": custom_prop_file, "cols": "1,2,3,4,5,6" # 使用自定义文件中的第2列到第7列(共6个性质) } # 运行 ifr.get_feature(input_file, output_file, params) ``` 注意:由于`iFeature`的PseAAC-General计算多个性质时,会为每个性质计算λ个自相关因子,然后再加上20个氨基酸组成(如果设置了,默认是包含的),所以特征维度为:20 + 6 * λ(因为6个性质,每个性质计算λ个自相关因子)。 另外,需要确保自定义属性文件中的氨基酸字母是标准的20种氨基酸(不包括特殊字母),如果有非标准氨基酸,需要处理。 如果序列中有非标准氨基酸,我们需要在自定义属性文件中为它们指定数值,或者将其映射为标准氨基酸,或者忽略。 由于我们使用的是标准性质,通常只考虑20种标准氨基酸。 步骤总结: 1. 收集7000个蛋白质序列,保存为FASTA文件。 2. 创建自定义属性文件,包含20种氨基酸的6个性质(疏水性、亲水性、分子量、pK1、pK2、等电点)。注意:这些性质的具体数值需要从可靠的生物化学数据源获取。 3. 安装`iFeature`(注意:`iFeature`需要Python3.6+,以及一些依赖包)。 4. 运行上述脚本,得到特征矩阵(每行一个序列,每列一个特征)。 特征提取完成后,我们将得到一个包含每个序列的PseAAC-General特征向量的文件。 注意:由于计算量较大,7000条序列可能需要一些时间(取决于序列长度和λ的大小)。 关于性质数值的参考: - 分子量:标准氨基酸分子量(例如,丙氨酸89.09,精氨酸174.20等) - pK1(α-羧基解离常数)和pK2(α-氨基解离常数)以及等电点(pI)可以从生物化学手册或数据库中获取。例如,我们参考下面的表格: 氨基酸 | pK1 (α-COOH) | pK2 (α-NH3+) | pI ------|-------------|-------------|---- Ala | 2.34 | 9.69 | 6.00 Arg | 2.17 | 9.04 | 10.76 ... | ... | ... | ... 疏水性和亲水性有多种标度,这里我们可以选择常用的标度,例如: - 疏水性:Kyte-Doolittle标度(在AAIndex中为HOPT810101) - 亲水性:Hopp-Woods标度(在AAIndex中为HOPT810101?注意:HOPT810101是疏水性,亲水性有其他的编号) 由于我们需要自定义,所以我们可以直接使用这些标准标度的数值。 下面是一个示例自定义属性文件(部分,注意:数值仅为示例,实际请使用标准数值): ``` A 1.8 -0.5 89.09 2.34 9.69 6.00 R -4.5 3.0 174.20 2.17 9.04 10.76 N -3.5 0.2 132.12 2.02 8.80 5.41 D -3.5 3.0 133.10 1.88 9.60 2.77 C 2.5 -1.0 121.15 1.96 10.28 5.07 Q -3.5 0.2 146.15 2.17 9.13 5.65 E -3.5 3.0 147.13 2.19 9.67 3.22 G -0.4 0.0 75.07 2.34 9.60 5.97 H -3.2 -0.5 155.16 1.82 9.17 7.59 I 4.5 -1.8 131.17 2.36 9.60 6.02 L 3.8 -1.8 131.17 2.36 9.60 5.98 K -3.9 3.0 146.19 2.18 8.95 9.74 M 1.9 -1.3 149.21 2.28 9.21 5.74 F 2.8 -2.5 165.19 1.83 9.13 5.48 P -1.6 0.0 115.13 1.99 10.60 6.30 S -0.8 0.3 105.09 2.21 9.15 5.68 T -0.7 -0.4 119.12 2.09 9.10 5.60 W -0.9 -3.4 204.23 2.83 9.39 5.89 Y -1.3 -2.3 181.19 2.20 9.11 5.66 V 4.2 -1.5 117.15 2.32 9.62 5.96 ``` 注意:这个示例中,第一列是氨基酸,第二列是疏水性(Kyte-Doolittle标度,数值越大越疏水),第三列是亲水性(这里使用了Hopp-Woods标度,数值越大越亲水),后面是分子量、pK1、pK2和pI。 实际数值请根据可靠的来源填写。 如果使用`iFeature`时遇到问题(例如自定义属性文件读取错误),可以查阅其文档或源码。 另外,`iFeature`还支持直接使用AAIndex中的性质,但我们这里需要多个自定义性质,所以使用自定义文件更合适。 如果不想使用`iFeature`,也可以考虑使用`protr`(R包)或自己编写代码。自己编写代码的思路如下: 1. 读取FASTA文件,得到序列列表。 2. 定义20种氨基酸的6个性质的字典(从自定义属性文件加载)。 3. 对于每条序列,计算: a) 20个氨基酸的组成(20维)。 b) 对于每个物理化学性质,计算序列的λ阶自相关函数(autocorrelation)。自相关函数定义为: AC(k, property) = (1/(L-k)) * Σ_{i=1}^{L-k} (P_i * P_{i+k}) 其中,k=1,2,...,λ;P_i表示序列中第i个氨基酸在该性质下的数值。 c) 将20个氨基酸组成和6*λ个自相关值组合成一个向量(20+6*λ维),并乘以权重因子(具体公式参考PseAAC原始论文)。 但是,自己实现需要处理归一化等问题,所以推荐使用现成工具。 最后,由于7000条序列较多,建议在服务器或高性能计算机上运行。 相关问题: 1. 如何获取每个氨基酸的疏水性、亲水性、分子量、pK1、pK2和等电点等物理化学性质的数值? 2. 在PseAAC-General中,如何选择λ和权重因子w? 3. 除了PseAAC-General,还有哪些特征提取方法可以用于蛋白质序列的物理化学性质分析? 4. 如何处理序列中的非标准氨基酸(如X、U等)? 5. 如何验证PseAAC特征的有效性?</think>### PseAAC特征提取方法详解 PseAAC(伪氨基酸组成)是一种将蛋白质序列的氨基酸组成与物理化学性质相结合的生物信息学方法。以下是对7000个蛋白质序列进行疏水性、亲水性、分子量、pK1(α-CO₂H)、pK2(NH₃)和等电点(pI)量化分析的完整流程: --- ### 步骤1:数据准备 1. **序列格式**:确保7000个蛋白质序列为FASTA格式 ```fasta >Protein1 MALWMRLLPL... >Protein2 MKFFIVFTNV... ``` 2. **氨基酸性质表**:创建物理化学性质参考表(示例如下)[^1]: | 氨基酸 | 疏水性 | 亲水性 | 分子量 | pK1(α-CO₂H) | pK2(NH₃) | 等电点(pI) | |--------|--------|--------|--------|-------------|-----------|------------| | A | 0.62 | -0.5 | 89.09 | 2.34 | 9.69 | 6.00 | | R | -2.53 | 3.00 | 174.20 | 2.17 | 9.04 | 10.76 | | ... | ... | ... | ... | ... | ... | ... | > 完整数据需参考权威数据库如[AAIndex](https://www.genome.jp/aaindex/)[^2] --- ### 步骤2:特征计算原理 PseAAC特征向量公式: $$ \mathbf{P} = [p_1, p_2, ..., p_{20}, p_{21}, ..., p_{20+\lambda}]^T $$ 其中: - 前20维:标准氨基酸频率 $ p_u = \frac{f_u}{\sum f_i + w \sum \tau_k} $ - 后$\lambda$维:序列相关因子 $ \tau_k = \frac{1}{L-k} \sum_{i=1}^{L-k} \Psi (R_i, R_{i+k}) $ $\Psi$函数整合六种性质: $$ \Psi (R_i, R_j) = \frac{1}{6} \sum_{m=1}^{6} [P_m(R_i) - P_m(R_j)]^2 $$ $P_m$分别对应:疏水性、亲水性、分子量、pK1、pK2、pI --- ### 步骤3:Python实现方案 ```python import numpy as np from Bio import SeqIO from sklearn.preprocessing import StandardScaler # 加载性质字典(示例片段) PROPERTY_DICT = { 'A': {'hydro': 0.62, 'philic': -0.5, 'mw': 89.09, 'pK1': 2.34, 'pK2': 9.69, 'pI': 6.00}, 'R': {'hydro': -2.53, 'philic': 3.0, 'mw': 174.2, 'pK1': 2.17, 'pK2': 9.04, 'pI': 10.76}, # 补充其他18种氨基酸... } def calculate_pseaac(sequence, lambda_val=5, w=0.05): L = len(sequence) # 1. 计算标准氨基酸组成 aa_counts = {aa: sequence.count(aa) for aa in "ACDEFGHIKLMNPQRSTVWY"} total_aa = sum(aa_counts.values()) # 2. 计算序列相关因子 tau = [] for k in range(1, lambda_val + 1): psi_sum = 0 for i in range(L - k): # 获取两个残基的六种性质 prop_i = [PROPERTY_DICT[sequence[i]][key] for key in ['hydro','philic','mw','pK1','pK2','pI']] prop_j = [PROPERTY_DICT[sequence[i+k]][key] for key in ['hydro','philic','mw','pK1','pK2','pI']] # 计算Ψ函数 psi_val = np.mean([(pi - pj)**2 for pi, pj in zip(prop_i, prop_j)]) psi_sum += psi_val tau.append(psi_sum / (L - k)) # 3. 组合特征向量 p_aa = [aa_counts[aa]/total_aa for aa in "ACDEFGHIKLMNPQRSTVWY"] denominator = sum(p_aa) + w * sum(tau) pseaac = [x/denominator for x in p_aa] + [w * t / denominator for t in tau] return np.array(pseaac) # 主处理流程 def process_7000_sequences(fasta_path, output_file): features = [] for record in SeqIO.parse(fasta_path, "fasta"): seq = str(record.seq).upper() # 过滤非标准氨基酸 seq = ''.join(aa for aa in seq if aa in PROPERTY_DICT) if len(seq) > 10: # 忽略过短序列 features.append(calculate_pseaac(seq)) # 标准化特征矩阵 scaler = StandardScaler() scaled_features = scaler.fit_transform(features) # 保存结果 (7000 x (20+lambda) 矩阵) np.save(output_file, scaled_features) # 执行 process_7000_sequences("proteins.fasta", "pseaac_features.npy") ``` --- ### 关键参数说明 1. **λ值选择**:序列相关因子阶数(默认5-30) - 建议:$$ \lambda < \min(L)/10 $$,$L$为序列长度 2. **权重因子w**:平衡组成与序列效应(默认0.05-0.7) 3. **特征维度**:$20 + \lambda$(当$\lambda=10$时,输出30维向量) --- ### 优化建议 1. **并行计算**:使用`joblib`加速7000条序列处理 ```python from joblib import Parallel, delayed results = Parallel(n_jobs=-1)(delayed(calculate_pseaac)(seq) for seq in sequences) ``` 2. **性质标准化**:不同性质量纲差异大,需Z-score标准化 3. **处理缺失值**:对非常见氨基酸(B、X、Z): - 方案1:替换为平均值 - 方案2:使用K-nearest neighbors插补 --- ### 验证方法 1. **特征相关性分析**:计算特征与蛋白质功能的Spearman相关系数 $$ \rho = 1 - \frac{6\sum d_i^2}{n(n^2-1)} $$ 2. **降维可视化**:使用t-SNE观察特征分布 ```python from sklearn.manifold import TSNE tsne = TSNE(n_components=2) vis_features = tsne.fit_transform(features) ``` --- ### 注意事项 1. **性质数据一致性**:必须使用同一来源的物化参数表 2. **序列长度影响**:短序列(<30aa)需增加λ权重 3. **内存优化**:7000×30矩阵约占用: $$ 7000 \times 30 \times 8 \text{ bytes} \approx 1.68 \text{ MB} $$ > 完整实现可参考GitHub仓库:[bio-pseaac-toolkit](https://github.com/bio-tools/pseaac)[^3]
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值