C : Unidentified if evaluation issue

Question

这比这更复杂，但简而言之，我正在尝试为序列族创建一个共识序列(仅由 A、C、G 和 T 字符组成的字符串)，并且我无法识别在哪里我做的功能失败了。这是:

SEQUENCE calculate_cons(FAMILY FAM)
{
    SEQUENCE cons;
    int maxlenght=seq_with_max_char(FAM);

    for(int i=0; i<maxlenght; i++)
    {
        int nA=0;
        int nC=0;
        int nG=0;
        int nT=0;

        for(int j=0; j<FAM.size; j++)
        {
            if(FAM.seq[j].c[i]=='A')
            {
                nA++;
            }
            if(FAM.seq[j].c[i]=='C')
            {
                nC++;
            }
            if(FAM.seq[j].c[i]=='G')
            {
                nG++;
            }
            if(FAM.seq[j].c[i]=='T')
            {
                nT++;
            }
        }

        if((nA==nC) || (nA==nG) || (nA==nT))
        {
            cons.c[i]='.';
        }
        else
        {
            if((nA>nC) && (nA>nG) && (nA>nT))
            {
                cons.c[i]='A';
            }
            if((nC>nA) && (nC>nG) && (nC>nT))
            {
                cons.c[i]='C';
            }
            if((nG>nA) && (nG>nC) && (nG>nT))
            {
                cons.c[i]='G';
            }
            if((nT>nA) && (nT>nC) && (nT>nG))
            {
                cons.c[i]='T';
            }
        }
    }

    cons.lenght=maxlenght;
    cons.ispartfam=true;

    return cons;
}

问题:使用此代码，共识序列仅由“A”和“.”组成。例如，如果一个系列包含:

TCCTATGGAATCTTTTTA

TTCTATGGAATCTTTTTA

共识序列将是:

...A...AA.......A

该函数写入“.”当不存在2次A时，否则写“A”。它失败的行可能是 if((nA==nC) || (nA==nG) || (nA==nT)) 因为如果我与 nC 进行比较，共识族将只包含 'C' 和 ' .'.

编辑:您将在下面找到一个最小的繁殖示例。在保持清晰的同时，我不能比这更简约了。

#include <stdio.h>

int main()
{
    int nbseq=2;
    int seqsize=5;
    char tseq[nbseq][seqsize];
    char cons[seqsize];

    tseq[0][0]='A';
    tseq[0][1]='C';
    tseq[0][2]='A';
    tseq[0][3]='T';
    tseq[0][4]='A';
    tseq[1][0]='G';
    tseq[1][1]='A';
    tseq[1][2]='A';
    tseq[1][3]='A';
    tseq[1][4]='G';

    for(int i=0; i<seqsize; i++)
    {
        int nA=0;
        int nC=0;
        int nG=0;
        int nT=0;

        for(int j=0; j<nbseq; j++)
        {
            if(tseq[j][i]=='A')
            {
                nA++;
            }
            if(tseq[j][i]=='C')
            {
                nC++;
            }
            if(tseq[j][i]=='G')
            {
                nG++;
            }
            if(tseq[j][i]=='T')
            {
                nT++;
            }
        }

        if((nA==nC) || (nA==nG) || (nA==nT))
        {
            cons[i]='.';
        }
        else
        {
            if((nA>nC) && (nA>nG) && (nA>nT))
            {
                cons[i]='A';
            }
            if((nC>nA) && (nC>nG) && (nC>nT))
            {
                cons[i]='C';
            }
            if((nG>nA) && (nG>nC) && (nG>nT))
            {
                cons[i]='G';
            }
            if((nT>nA) && (nT>nC) && (nT>nG))
            {
                cons[i]='T';
            }
        }
    }

    for(int i=0; i<nbseq; i++)
    {
        for(int j=0; j<seqsize; j++)
        {
        printf("%c", tseq[i][j]);
        }
        printf("\n\n");
    }

    for(int i=0; i<seqsize; i++)
    {
        printf("%c", cons[i]);
    }

    return 0; 
}

Answer 1

我同意有问题的路线是

       if((nA==nC) || (nA==nG) || (nA==nT))

我认为这种逻辑背后的想法应该是，您需要在给定位置上有一个唯一的最大碱基数才能达成共识，因此必须拒绝平局。不幸的是，这个想法既有缺陷，也有错误的实现。

它是有缺陷的，因为平局在一般情况下并不重要，而只有在给定位置存在最大计数的平局时才是重要的。特别注意，如果给定位置的每个碱基都是 A，那么这些条件都不成立，但如果每个碱基都是其他碱基中的任何一个，那么其中两个都将评估为 true，形式为 0 == 0 .

它的实现不正确，因为您没有测试所有对。有六种不同的方法可以从一组四个不同的元素中形成对，但您只提供了三种。缺少的可能是 nC==nG、nC==nT 和 nG==nT。如果您也包含了这些内容，那么您的输出将仅包含“.”，但您恰好省略了允许共识 A 通过的内容。

更糟糕的是，这个特定的部分毫无意义。您不需要预先检查关系，因为尝试识别唯一最大值的部分无论如何都会排除关系(对于最大值)。那么，更有效的方法就是简单地设置“.”，而不是预先检查关系。当没有满足真正的共识条件时作为后备。有多种方法可以做到这一点，例如:

            # Skip the explicit tie-checking altogether, and just do this ...
            if ((nA > nC) && (nA > nG) && (nA > nT)) {
                cons.c[i] = 'A';
            } else if ((nC > nA) && (nC > nG) && (nC > nT)) {
                cons.c[i] = 'C';
            } else if ((nG > nA) && (nG > nC) && (nG > nT)) {
                cons.c[i] = 'G';
            } else if ((nT > nA) && (nT > nC) && (nT > nG)) {
                cons.c[i] = 'T';
            } else {
                cons.c[i] = '.';
            }