C++ DNASequence类代码示例

// diagnose
//   SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS.....................................................
//   ..........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX......................
//   ...............................IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII......................
//   .................................JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ......................
//   LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL....................................................
//   !"#$%&'()*+,-./0123456789:;<=>[email protected][\]^_`abcdefghijklmnopqrstuvwxyz{|}~
//   |                         |    |        |                              |                     |
//  33                        59   64       73                            104                   126 <- maxValue is value from here
// S 0........................26...31.......40
// X                          -5....0........9.............................40
// I                                0........9.............................40
// J                                   3.....9.............................40
// L 0.2......................26...31........41
//
//  S - Sanger        Phred+33,  raw reads typically (0, 40)
//  X - Solexa        Solexa+64, raw reads typically (-5, 40)
//  I - Illumina 1.3+ Phred+64,  raw reads typically (0, 40)
//  J - Illumina 1.5+ Phred+64,  raw reads typically (3, 40) with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator (bold)
//  L - Illumina 1.8+ Phred+33,  raw reads typically (0, 41)
DNAQualityType FastqQualityTrimTask::detectQualityType(){
    int maxValue = 33;
    int minValue = 126;
    FASTQIterator iter_qual(settings.inputUrl, stateInfo);
    CHECK(!stateInfo.hasError(), DNAQualityType_Sanger);

    int counter = 0;
    while (iter_qual.hasNext()) {
        CHECK(!stateInfo.isCoR(), DNAQualityType_Sanger);

        if (counter > 1000) {   // check only first 1000 reads in file
            break;
        }

        DNASequence dna = iter_qual.next();
        int seqLen = dna.length();
        if (seqLen > dna.quality.qualCodes.length()) {
            continue;
        } else {
            for (int pos = 0; pos <= seqLen - 1; pos++) {
                maxValue = qMax(static_cast<int>(dna.quality.qualCodes.at(pos)), maxValue);
                minValue = qMin(static_cast<int>(dna.quality.qualCodes.at(pos)), minValue);
            }
        }
        counter++;
    }
    return DNAQuality::detectTypeByMinMaxQualityValues(minValue, maxValue);
}

void CASAVAFilterTask::runStep(){
    int ncount = 0;
    int ycount = 0;

    QScopedPointer<IOAdapter> io  (IOAdapterUtils::open(settings.outDir + settings.outName, stateInfo, IOAdapterMode_Append));

    //1:N:0:TAAGGG
    QRegExp pattern (":Y:[^:]:");
    FASTQIterator iter(settings.inputUrl);
    while(iter.hasNext()){
        if(stateInfo.isCoR()){
            return;
        }
        DNASequence seq = iter.next();
        QString comment = DNAInfo::getFastqComment(seq.info);
        if(pattern.indexIn(comment) != -1){
            ycount++;
        }else{
            FastqFormat::writeEntry(seq.getName() + " " + comment, seq, io.data(), "Writing error", stateInfo, false);
            ncount++;
        }
    }

    algoLog.info(QString("Discarded by CASAVA filter %1").arg(ycount));
    algoLog.info(QString("Accepted by CASAVA filter %1").arg(ncount));
    algoLog.info(QString("Total by CASAVA FILTER: %1").arg(ncount + ycount));
}

//Test DNASequence Allocate(DNALength)
TEST_F(DNASequenceTest, Allocate) {
    dnaOne.Allocate(0);
    EXPECT_EQ(dnaOne.length, 0);

    DNASequence dnaTwo;
    dnaTwo.Allocate(100);
    EXPECT_EQ(dnaTwo.length, 100);
}

int main(int argc, char* argv[]) {
	string refFileName, queryFileName;
	int maxHammingDistance;
	if (argc < 4) {
		cout << "usage: hammer ref query maxHam " << endl;
		exit(1);
	}
	refFileName = argv[1];
	queryFileName = argv[2];
	maxHammingDistance = atoi(argv[3]);

	FASTAReader reader;
	reader.Initialize(refFileName);
	FASTASequence ref, refRC;
	reader.GetNext(ref);
	ref.MakeRC(refRC);
	
	FASTAReader queryReader;
	queryReader.Initialize(queryFileName);
	FASTASequence query;
	queryReader.GetNext(query);
	DNALength p;
	for(p=0; p < ref.length-query.length-1; p++ ){
		DNASequence subseq;
		subseq.seq = &ref.seq[p];
		subseq.length = query.length;
		//		cout << "t "; subseq.PrintSeq(cout);
		//		cout << "q "; ((DNASequence*)&query)->PrintSeq(cout);
		if (HammingDistance(&subseq.seq[0], &query.seq[0], query.length) < maxHammingDistance) {
			cout << ">" << p << endl;
			subseq.PrintSeq(cout);
		}
		int i;
		for (i =0; i < query.length; i++) {
			subseq.seq[i] = toupper(subseq.seq[i]);
		}
	}

	for(p=0; p < ref.length-query.length-1; p++ ){
		DNASequence subseq;
		subseq.seq = &refRC.seq[p];
		subseq.length = query.length;
		if (HammingDistance(&subseq.seq[0], &query.seq[0], query.length) < maxHammingDistance) {
			cout << ">" << p << "rc" << endl;
			subseq.PrintSeq(cout);
		}
		int i;
		for (i =0; i < query.length; i++) {
			subseq.seq[i] = toupper(subseq.seq[i]);
		}
	}

}

bool do_test(string sequence, int __expected) {
    time_t startClock = clock();
    DNASequence *instance = new DNASequence();
    int __result = instance->longestDNASequence(sequence);
    double elapsed = (double)(clock() - startClock) / CLOCKS_PER_SEC;
    delete instance;

    if (__result == __expected) {
        cout << "PASSED!" << " (" << elapsed << " seconds)" << endl;
        return true;
    }
    else {
        cout << "FAILED!" << " (" << elapsed << " seconds)" << endl;
        cout << "           Expected: " << to_string(__expected) << endl;
        cout << "           Received: " << to_string(__result) << endl;
        return false;
    }
}

static void saveSequence(IOAdapter* io, const DNASequence& sequence, U2OpStatus& os) {
    writeHeaderToFile( io, sequence.getName( ), os );
    CHECK_OP( os, );

    const char *seq = sequence.seq.constData( );
    const int len = sequence.seq.length( );
    for ( int i = 0; i < len; i += SAVE_LINE_LEN ) {
        const int chunkSize = qMin( SAVE_LINE_LEN, len - i );
        writeBlockToFile( io, seq + i, chunkSize, os );
        CHECK_OP( os, );
    }
}

Task* PWMatrixSearchWorker::tick() {
    while (modelPort->hasMessage()) {
        models << modelPort->get().getData().toMap().value(PWMatrixWorkerFactory::WMATRIX_SLOT.getId()).value<PWMatrix>();
    }
    if (!modelPort->isEnded()) {
        return NULL;
    }

    if (dataPort->hasMessage()) {
        Message inputMessage = getMessageAndSetupScriptValues(dataPort);
        if (inputMessage.isEmpty() || models.isEmpty()) {
            output->transit();
            return NULL;
        }
        QVariantMap map = inputMessage.getData().toMap();
        SharedDbiDataHandler seqId = map.value(BaseSlots::DNA_SEQUENCE_SLOT().getId()).value<SharedDbiDataHandler>();
        QScopedPointer<U2SequenceObject> seqObj(StorageUtils::getSequenceObject(context->getDataStorage(), seqId));
        if (seqObj.isNull()) {
            return NULL;
        }
        U2OpStatusImpl os;
        DNASequence seq = seqObj->getWholeSequence(os);
        CHECK_OP(os, new FailTask(os.getError()));

        if (!seq.isNull() && seq.alphabet->getType() == DNAAlphabet_NUCL) {
            WeightMatrixSearchCfg config(cfg);
            config.complOnly = (strand < 0);
            if (strand <= 0) {
                DNATranslation* compTT = AppContext::getDNATranslationRegistry()->
                    lookupComplementTranslation(seq.alphabet);
                if (compTT  != NULL) {
                    config.complTT = compTT  ;
                }
            }
            QList<Task*> subtasks;
            foreach(PWMatrix model, models) {
                subtasks << new WeightMatrixSingleSearchTask(model, seq.seq, config, 0);
            }

void FastqQualityTrimTask::runStep(){
    int ncount = 0;
    int ycount = 0;

    QScopedPointer<IOAdapter> io(IOAdapterUtils::open(settings.outDir + settings.outName, stateInfo, IOAdapterMode_Append));

    int quality = settings.customParameters.value(QUALITY_ID, 20).toInt();
    int minLen = settings.customParameters.value(LEN_ID, 0).toInt();
    bool bothEnds = settings.customParameters.value(BOTH_ID, false).toInt();
    DNAQualityType qualityType = detectQualityType();
    CHECK_OP(stateInfo, );

    FASTQIterator iter(settings.inputUrl, stateInfo);
    CHECK_OP(stateInfo, );

    while (iter.hasNext()) {
        CHECK_OP(stateInfo, );

        DNASequence dna = iter.next();
        dna.quality.type = qualityType;
        const U2Region acceptedRegion = DNASequenceUtils::trimByQuality(dna, quality, minLen, bothEnds);

        if (0 < acceptedRegion.length) {
            ycount++;
        } else {
            ncount++;
            continue;
        }

        FastqFormat::writeEntry(dna.getName(), dna, io.data(), "Writing error", stateInfo, false);
    }

    algoLog.info(QString("Discarded by trimmer %1").arg(ncount));
    algoLog.info(QString("Accepted by trimmer %1").arg(ycount));
    algoLog.info(QString("Total by trimmer %1").arg(ncount + ycount));
}

GeneByGeneCompareResult GeneByGeneComparator::compareGeneAnnotation(const DNASequence& seq, const QList<SharedAnnotationData> &annData,
    const QString& annName, float identity)
{
    GeneByGeneCompareResult result;

    float maxIdentity = -1.0F;
    foreach (const SharedAnnotationData &adata, annData) {
        if (adata->name == annName) {
            U2Location location = adata->location;
            if (location->isSingleRegion()) {
                int reglen = location->regions.first().length;
                float lenRatio  = reglen * 100 /static_cast<float>(seq.length());
                maxIdentity = qMax(maxIdentity, lenRatio);
                if(lenRatio >= identity){ //check length ratio
                    QString ident = adata->findFirstQualifierValue(BLAST_IDENT);
                    if (!ident.isEmpty()){
                        //create BLAST string  YES/identity/gaps
                        float blastIdent = parseBlastQual(ident);
                        if (blastIdent != -1.0f && blastIdent >= identity){
                            result.identical = true;
                            result.identityString = GeneByGeneCompareResult::IDENTICAL_YES;
                            result.identityString.append(QString("\\%1").arg(blastIdent));
                            QString gaps = adata->findFirstQualifierValue(BLAST_GAPS);
                            if (!gaps.isEmpty()){
                                float blastGaps = parseBlastQual(gaps);
                                if (blastGaps!=1.0f){
                                    result.identityString.append(QString("\\%1").arg(blastGaps));
                                }
                            }else{
                                result.identityString.append(QString("\\0"));
                            }
                        }
                    }else{ //not a blast annotation
                        result.identical = true;
                        result.identityString = GeneByGeneCompareResult::IDENTICAL_YES;
                    }
                }
            }
            break;
        }
    }

    if (result.identical == false && maxIdentity != -1.0f){
        result.identityString.append(QString("\\%1").arg(maxIdentity));
    }

    return result;
}

int main(int argc, char* argv[]) {
    string inFileName, readsFileName;
    DNALength readLength;
    float coverage = 0;
    bool noRandInit = false;
    int numReads = -1;
    CommandLineParser clp;
    int qualityValue = 20;
    bool printFastq = false;
    int stratify = 0;
    string titleType = "pacbio";
    string fastqType = "illumina"; // or "sanger"
    clp.RegisterStringOption("inFile", &inFileName, "Reference sequence", 0);
    clp.RegisterPreviousFlagsAsHidden();
    clp.RegisterIntOption("readLength", (int*) &readLength, "The length of reads to simulate.  The length is fixed.",
                          CommandLineParser::PositiveInteger, "Length of every read.", 0);
    clp.RegisterFloatOption("coverage", &coverage, "Total coverage (from which the number of reads is calculated",
                            CommandLineParser::PositiveFloat, 0);
    clp.RegisterFlagOption("nonRandInit", &noRandInit, "Skip initializing the random number generator with time.");
    clp.RegisterIntOption("nReads", &numReads, "Total number of reads (from which coverage is calculated)", CommandLineParser::PositiveInteger, 0);
    clp.RegisterStringOption("readsFile", &readsFileName, "Reads output file", 0);
    clp.RegisterFlagOption("fastq", &printFastq, "Fake fastq output with constant quality value (20)");
    clp.RegisterIntOption("quality", &qualityValue, "Value to use for fastq quality", CommandLineParser::PositiveInteger);
    clp.RegisterIntOption("stratify", &stratify, "Sample a read every 'stratify' bases, rather than randomly.", CommandLineParser::PositiveInteger);
    clp.RegisterStringOption("titleType", &titleType, "Set the name of the title: 'pacbio'|'illumina'");
    clp.RegisterStringOption("fastqType", &fastqType, "Set the type of fastq: 'illumina'|'sanger'");
    vector<string> leftovers;
    clp.ParseCommandLine(argc, argv, leftovers);

    if (!noRandInit) {
        InitializeRandomGeneratorWithTime();
    }

    FASTAReader inReader;
    inReader.Init(inFileName);
    vector<FASTASequence> reference;

    inReader.ReadAllSequences(reference);
    ofstream readsFile;
    if (readsFileName == "") {
        cout << "ERROR.  You must specify a reads file." << endl;
        exit(0);
    }
    CrucialOpen(readsFileName, readsFile, std::ios::out);

    ofstream sangerFastqFile;
    if (fastqType == "sanger") {
        string sangerFastqFileName = readsFileName + ".fastq";
        CrucialOpen(sangerFastqFileName, sangerFastqFile, std::ios::out);
    }

    DNALength refLength = 0;
    int i;
    for (i = 0; i < reference.size(); i++) {
        refLength += reference[i].length;
    }
    if (numReads == -1 and coverage == 0 and stratify == 0) {
        cout << "Error, you must specify either coverage, nReads, or stratify." << endl;
        exit(1);
    }
    else if (numReads == -1) {
        numReads = (refLength / readLength) * coverage;
    }

    if (stratify) {
        if (!readLength) {
            cout << "ERROR. If you are using stratification, a read length must be specified." << endl;
            exit(1);
        }
    }

    DNASequence sampleSeq;
    sampleSeq.length = readLength;
    int maxRetry = 10000000;
    int retryNumber = 0;
    DNALength seqIndex, seqPos;
    if (stratify) {
        seqIndex = 0;
        seqPos   = 0;
    }
    DNALength origReadLength = readLength;
    for (i = 0; stratify or i < numReads; i++) {
        if (stratify == 0) {
            FindRandomPos(reference, seqIndex, seqPos, readLength );
        }
        else {
            //
            // find the next start pos, or bail if done
            //
            if (seqPos >= reference[seqIndex].length) {
                if (seqIndex == reference.size() - 1) {
                    break;
                }
                else {
                    seqIndex = seqIndex + 1;
                    seqPos   = 0;
                    continue;
                }
            }
            readLength = min(reference[seqIndex].length - seqPos, origReadLength);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	ofstream summaryTable(summaryTableFileName.c_str());
	ofstream outputFile(outputFileName.c_str());

	queryReader.Initialize(queryFileName);

	while (queryReader.GetNext(querySequence)) {
		int i;
		cerr << "searching " << querySequence.title << endl;
		if (querySequence.length < maxLength) {
			continue;
		}

		int nMatches = 0;
		querySequence.ToUpper();
		int localMax;
		for (i = 0; i < querySequence.length - maxLength + 1; i++) {
			if ((i + 1) % 100000 == 0) {
				cerr << "processed: " << i + 1 << endl;
			}

			int lcpLength;
			vector<SAIndex> lcpLeftBounds, lcpRightBounds;
			vector<SAIndex> rclcpLeftBounds, rclcpRightBounds;
			localMax = maxN;
			if (i < prefix) {
				localMax = prefixN;
			}
			if (i >= querySequence.length - suffix) {
				localMax = suffixN;
			}
			if (querySequence.length - i <= maxLength) {
				continue;
			}
			if (querySequence.seq[i] == 'N') {
				continue;
			}
			lcpLength = sarray.StoreLCPBounds(genome.seq, genome.length, // The string which the suffix array is built on.
																				&querySequence.seq[i], querySequence.length-i,
																				true,
																				maxLength,
																				lcpLeftBounds, lcpRightBounds,
																				false);
			if (lcpLength < minLength) {
				continue;
			}
			if (lcpLength < maxLength or 
					lcpRightBounds.size() == 0 or 
					(lcpRightBounds.size() > 0 and 
					 lcpLeftBounds.size() > 0 and  
					 lcpRightBounds[lcpRightBounds.size() - 1] - lcpLeftBounds[lcpLeftBounds.size()-1] <= localMax)) {

				FASTASequence rc;
				DNASequence subseq;
				subseq.ReferenceSubstring(querySequence, i, maxLength);
				subseq.MakeRC(rc);
				int rclcpLength;
				int numForwardMatches;
				if (lcpLength == 0) {
					numForwardMatches = 0;
				}
				else {
					numForwardMatches = lcpRightBounds[lcpRightBounds.size() - 1] - lcpLeftBounds[lcpLeftBounds.size()-1];
				}
				rclcpLength = sarray.StoreLCPBounds(genome.seq, genome.length, // The string which the suffix array is built on.
																						rc.seq, maxLength,
																						true,
																						rclcpLength,
																						rclcpLeftBounds, rclcpRightBounds,
																						false);

				string rcstr((const char*)rc.seq, rc.length);

				if (rclcpLength < maxLength or 
						rclcpRightBounds.size() == 0 or
						(numForwardMatches + 
						 rclcpRightBounds[rclcpRightBounds.size() - 1] -
						 rclcpLeftBounds[rclcpLeftBounds.size()-1] <= localMax)) 
					{
						char* substr = new char[maxLength+1];
						substr[maxLength] = '\0';
						memcpy(substr, &querySequence.seq[i], maxLength);

						//						string substr = string((const char*) querySequence.seq, i, maxLength);
						
						outputFile << querySequence.title << "\t" << substr << "\t" << i << endl;

						++nMatches;
						delete[] substr;
						//					}
					}
				rc.Free();
			}

		}
		summaryTable << querySequence.title << "\t" << nMatches << endl;
		querySequence.Free();
	}
	outputFile.close();
	genome.Free();
}

void SAMAlignmentsToCandidates(SAMAlignment &sam,
                               std::vector<FASTASequence> &referenceSequences,
                               std::map<std::string,int> & refNameToRefListIndex,
                               std::vector<AlignmentCandidate<> > &candidates, 
                               bool parseSmrtTitle,
                               bool keepRefAsForward,
                               bool copyQVs) {
  //
  // First determine how many alignments there are from CIGAR string.
  //
  std::vector<int> lengths;
  std::vector<char> ops;
  sam.cigar.Vectorize(lengths, ops);

  DNASequence querySeq;
  // For now just reference the query sequence.
  querySeq.deleteOnExit = false;
  querySeq.seq = (Nucleotide*) sam.seq.c_str();
  querySeq.length = sam.seq.size();

  DNALength samTEnd = 0;
  DNALength samTStart = sam.pos - 1;

  std::vector<std::string> optionalQVs;
  if (copyQVs) {
      sam.CopyQVs(&optionalQVs);    
  }
  if (keepRefAsForward == false and IsReverseComplement(sam.flag)) {
    ReverseAlignmentOperations(lengths, ops);
    DNASequence rcQuerySeq;
    querySeq.CopyAsRC(rcQuerySeq);
    //
    // Zero out the query seq so that the string memory is not
    // deleted.
    //
    querySeq.seq = NULL;
    querySeq.length = 0;
    querySeq = rcQuerySeq;
    rcQuerySeq.Free();
    samTEnd = GetAlignedReferenceLengthByCIGARSum(ops, lengths);
    
    // We also need to reverse any optional QVs
    if (copyQVs) {
      for(int i=0; i<optionalQVs.size(); i++) {
        std::reverse(optionalQVs[i].begin(), optionalQVs[i].end());
      }
    }
  }


  int i;
  int offset = 0;
  if (ops.size() == 0) {
    return;
  }
  bool alignmentStarted = false;
  bool onFirstMatch = true;
  int  curAlignment;
  
  //
  // Advance past any clipping.  This advances in both query and
  // reference position.
  //
  int cigarPos = 0;
  int qPos = 0; 
  int tPos = 0;

  DNALength queryPosOffset = 0;
  if (parseSmrtTitle) {
    //
    // The aligned sequence is really a subread of a full
    // sequence. The position of the aligments start at 0, the
    // beginning of the query sequence, but in the sam file, they
    // may appear as subreads, and are offset from the start of the
    // subread.  By convention, the subread coordinates are embedded
    // in the title of the query, if it is a smrtTitle. 
    // Two types of smrtTitle are supported:
    // movie/zmw/start_end
    // movie/zmw/start_end/start2_end2
    SMRTTitle stitle = SMRTTitle(sam.qName);

    if (not stitle.isSMRTTitle) {
      std::cout << "ERROR. Could not parse title " << sam.qName << std::endl;
      exit(1);
    }
    queryPosOffset = stitle.start;
  }
  else if (sam.xs) {
    queryPosOffset += sam.xs - 1;
  }


  while (cigarPos < lengths.size()) {
    int numClipped;
    //
    // Sequence clipping becomes offsets into the q/t alignedSeqPos
    //


    int numSoftClipped;
//.........这里部分代码省略.........