C++读取FASTA文件时丢失最后一个字符的问题排查

问题描述 #

我正在编写一个C++程序读取FASTA文件并进行处理，FASTA文件格式如下：

以">"开头的是头部行，需要跳过/忽略
下方的行是我们需要的序列信息
ATTGGTATGATTTACCCAATTTGGGGAAAAAATTCCCTCTCGATAGCTATCCTGATTTGCGG
ATTGGTATGATTTACCCAATTTGGGGAAAAAATTCCCTCTCGATAGCTATCCTGATTTGCGG
ATTGGTATGATTTACCCAATTTGGGGAAAAAATTCCCTCTCGATAGCTATCCTGATTTGCGG

理想情况下，程序应跳过以>开头的头部行，将下方的序列内容读取到一个字符串中。但我的代码执行后，会丢失序列的最后一个字符，例如上述示例中最后一行的最后一个G不会被读取。

我的代码 #

void reading_in_RNA_file()
{

    string RNA_file = "sample_query.txt";

    ifstream fin;
        fin.open(RNA_file);
        if (!fin.is_open())
        {//if
                cerr << "Error did not open file" << endl;
                exit(1);
        }//if

    string line = "";
    string RNA_seq = "";
    string FASTA_heading = "";
    string sequence = "";

        while(getline(fin,line))
        {
                if( line.empty() || line[0] == '>' )
                { // Identifier marker
                        if(!FASTA_heading.empty() )
                        { // Print out what we read from the last entry
                                FASTA_heading.clear();
                                RNA_seq += sequence;
                        }

                        if( !line.empty() )
                        {
                                 FASTA_heading = line.substr(1);
                        }
                        sequence.clear();
                }

                else if(!FASTA_heading.empty())
                {
                        line = line.substr(0, line.length() -1);
                        if(line.find(' ') != string::npos )
                        { // Invalid sequence--no spaces allowed
                                 FASTA_heading.clear();
                                sequence.clear();
                        }

                        else
                        {
                            sequence += line;
                        }
                }
        }

        if(!FASTA_heading.empty() )
        { // Print out what we read from the last entry
                RNA_seq += sequence;
        }
    cout << RNA_seq << endl;
}

示例文件sample_query.txt内容 #

true positive test query
GTCTGAGAAAACAAGGCTAGAGATTCCAATATTAGAGACAACAGGGCTCTGGGAAGATTAAGGTTGAGTT
TTCTGGATCTGCAGAATAGAGTCACTGAGGACCAATTGCAAGATCAGAGGAGATGAAAGAACAAGTCAAG
GCATGCTTAGGAAAAGAGAATATCAGGGATAGGTTTTAGGCAAGAGTCACACTGAGGAAGGGCAGGTTCT
ACATACAGTTTATCTTGGTACTGCCAAGTACCATTTGGGTCAGGATTTTGTCATTTAGATCCATATTTTT
CCTATATTTTTATCTGGTTCTTCCATCAGTTACTGAGAGAGCACTATTAATTCACCAGCTATAATTTTGG
ATTGTCAATTTCCTGCTTTTGTCTGTTGTTTTTGATTCACATACTTTGAGGCTCTGTGTGTGTGTGTAAT

问题原因及修复 #

问题根源在这行代码：

line = line.substr(0, line.length() -1);

getline函数读取行时，会自动丢弃行尾的换行符，此时line中存储的是不包含换行的完整行内容。你执行substr(0, line.length()-1)，相当于把当前行的最后一个有效字符删掉了——最后一行的末尾字符也被截断，自然就出现了序列丢失最后一个字符的问题。

修复方案 #

直接删除这行错误的截取代码即可。修改后的序列处理部分代码如下：

else if(!FASTA_heading.empty())
{
        if(line.find(' ') != string::npos )
        { // Invalid sequence--no spaces allowed
                 FASTA_heading.clear();
                sequence.clear();
        }

        else
        {
            sequence += line;
        }
}

修改后程序就能完整读取所有序列字符了。

内容的提问来源于stack exchange，提问作者TheLordRaj