首页 > 解决方案 > C : 未确定是否评估问题

问题描述

这比那更复杂,但简而言之,我正在尝试为序列族(仅由 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; 
}

标签: cif-statement

解决方案


我同意有问题的路线是

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

我认为这个逻辑背后的想法应该是你需要一个给定位置的唯一最大碱基数来达成共识,因此必须拒绝平局。不幸的是,这个想法既有缺陷,也没有正确实施。

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

因为您没有测试所有的对,所以它被错误地实现了。有六种不同的方法可以从一组四个不同的元素中形成对,但您只提供了三种。缺少的是nC==nG,nC==nTnG==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] = '.';
            }

推荐阅读