CommandParameter prqual("rqfile", "InputTypes", "", "", "none", "none", "qfileGroup","",false,false,true); parameters.push_back(prqual);
CommandParameter pfile("file", "InputTypes", "", "", "FastaFastqFile", "FastaFastqFile", "none","fasta-qfile",false,false,true); parameters.push_back(pfile);
CommandParameter poligos("oligos", "InputTypes", "", "", "none", "none", "none","group",false,false,true); parameters.push_back(poligos);
+ CommandParameter pfindex("findex", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(pfindex);
+ CommandParameter prindex("rindex", "InputTypes", "", "", "none", "none", "none","",false,false,true); parameters.push_back(prindex);
CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs);
CommandParameter pbdiffs("bdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(pbdiffs);
CommandParameter ptdiffs("tdiffs", "Number", "", "0", "", "", "","",false,false); parameters.push_back(ptdiffs);
-
+ CommandParameter preorient("checkorient", "Boolean", "", "F", "", "", "","",false,false,true); parameters.push_back(preorient);
CommandParameter palign("align", "Multiple", "needleman-gotoh", "needleman", "", "", "","",false,false); parameters.push_back(palign);
CommandParameter pallfiles("allfiles", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pallfiles);
CommandParameter ptrimoverlap("trimoverlap", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(ptrimoverlap);
string MakeContigsCommand::getHelpString(){
try {
string helpString = "";
- helpString += "The make.contigs command reads a file, forward fastq file and a reverse fastq file or forward fasta and reverse fasta files and outputs new fasta. It will also provide new quality files if the fastq or file parameter is used.\n";
+ helpString += "The make.contigs command reads a file, forward fastq file and a reverse fastq file or forward fasta and reverse fasta files and outputs new fasta. \n";
helpString += "If an oligos file is provided barcodes and primers will be trimmed, and a group file will be created.\n";
- helpString += "The make.contigs command parameters are file, ffastq, rfastq, ffasta, rfasta, fqfile, rqfile, oligos, format, tdiffs, bdiffs, pdiffs, align, match, mismatch, gapopen, gapextend, insert, deltaq, allfiles and processors.\n";
+ helpString += "If a forward index or reverse index file is provided barcodes be trimmed, and a group file will be created. The oligos parameter is required if an index file is given.\n";
+ helpString += "The make.contigs command parameters are file, ffastq, rfastq, ffasta, rfasta, fqfile, rqfile, oligos, findex, rindex, format, tdiffs, bdiffs, pdiffs, align, match, mismatch, gapopen, gapextend, insert, deltaq, allfiles and processors.\n";
helpString += "The ffastq and rfastq, file, or ffasta and rfasta parameters are required.\n";
- helpString += "The file parameter is 2 or 3 column file containing the forward fastq files in the first column and their matching reverse fastq files in the second column, or a groupName then forward fastq file and reverse fastq file. Mothur will process each pair and create a combined fasta and report file with all the sequences.\n";
+ helpString += "The file parameter is 2, 3 or 4 column file containing the forward fastq files in the first column and their matching reverse fastq files in the second column, or a groupName then forward fastq file and reverse fastq file, or forward fastq file then reverse fastq then forward index and reverse index file. If you only have one index file add 'none' for the other one. Mothur will process each pair and create a combined fasta and report file with all the sequences.\n";
helpString += "The ffastq and rfastq parameters are used to provide a forward fastq and reverse fastq file to process. If you provide one, you must provide the other.\n";
helpString += "The ffasta and rfasta parameters are used to provide a forward fasta and reverse fasta file to process. If you provide one, you must provide the other.\n";
helpString += "The fqfile and rqfile parameters are used to provide a forward quality and reverse quality files to process with the ffasta and rfasta parameters. If you provide one, you must provide the other.\n";
helpString += "The format parameter is used to indicate whether your sequences are sanger, solexa, illumina1.8+ or illumina, default=illumina1.8+.\n";
+ helpString += "The findex and rindex parameters are used to provide a forward index and reverse index files to process. \n";
helpString += "The align parameter allows you to specify the alignment method to use. Your options are: gotoh and needleman. The default is needleman.\n";
helpString += "The tdiffs parameter is used to specify the total number of differences allowed in the sequence. The default is pdiffs + bdiffs + sdiffs + ldiffs.\n";
helpString += "The bdiffs parameter is used to specify the number of differences allowed in the barcode. The default is 0.\n";
//helpString += "The sdiffs parameter is used to specify the number of differences allowed in the spacer. The default is 0.\n";
helpString += "The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n";
helpString += "The mistmatch parameter allows you to specify the penalty for having different bases. The default is -1.0.\n";
+ helpString += "The checkorient parameter will check look for the reverse compliment of the barcode or primer in the sequence. If found the sequence is flipped. The default is false.\n";
helpString += "The deltaq parameter allows you to specify the delta allowed between quality scores of a mismatched base. For example in the overlap, if deltaq=5 and in the alignment seqA, pos 200 has a quality score of 30 and the same position in seqB has a quality score of 20, you take the base from seqA (30-20 >= 5). If the quality score in seqB is 28 then the base in the consensus will be an N (30-28<5) The default is 6.\n";
helpString += "The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -2.0.\n";
helpString += "The gapextend parameter allows you to specify the penalty for extending a gap in an alignment. The default is -1.0.\n";
//if the user changes the input directory command factory will send this info to us in the output parameter
- string inputDir = validParameter.validFile(parameters, "inputdir", false);
+ inputDir = validParameter.validFile(parameters, "inputdir", false);
if (inputDir == "not found"){ inputDir = ""; }
else {
string path;
//if the user has not given a path then, add inputdir. else leave path alone.
if (path == "") { parameters["oligos"] = inputDir + it->second; }
}
+
+ it = parameters.find("findex");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["findex"] = inputDir + it->second; }
+ }
+
+ it = parameters.find("rindex");
+ //user has given a template file
+ if(it != parameters.end()){
+ path = m->hasPath(it->second);
+ //if the user has not given a path then, add inputdir. else leave path alone.
+ if (path == "") { parameters["rindex"] = inputDir + it->second; }
+ }
}
ffastqfile = validParameter.validFile(parameters, "ffastq", true);
else if(oligosfile == "not open") { abort = true; }
else { m->setOligosFile(oligosfile); }
+ findexfile = validParameter.validFile(parameters, "findex", true);
+ if (findexfile == "not found") { findexfile = ""; }
+ else if(findexfile == "not open") { abort = true; }
+
+ rindexfile = validParameter.validFile(parameters, "rindex", true);
+ if (rindexfile == "not found") { rindexfile = ""; }
+ else if(rindexfile == "not open") { abort = true; }
+
+ if ((rindexfile != "") || (findexfile != "")) {
+ if (oligosfile == ""){
+ oligosfile = m->getOligosFile();
+ if (oligosfile != "") { m->mothurOut("Using " + oligosfile + " as input file for the oligos parameter.\n"); }
+ else {
+ m->mothurOut("You need to provide an oligos file if you are going to use an index file.\n"); abort = true;
+ }
+ }
+
+ //can only use an index file with the fastq parameters not fasta and qual
+ if ((ffastafile != "") || (rfastafile != "")) {
+ m->mothurOut("[ERROR]: You can only use an index file with the fastq parameters or the file option.\n"); abort = true;
+ }
+ }
+
//if the user changes the output directory command factory will send this info to us in the output parameter
outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
outputDir = "";
abort=true;
}
+ temp = validParameter.validFile(parameters, "checkorient", false); if (temp == "not found") { temp = "F"; }
+ reorient = m->isTrue(temp);
+
//fill convert table - goes from solexa to sanger. Used fq_all2std.pl as a reference.
for (int i = -64; i < 65; i++) {
char temp = (char) ((int)(33 + 10*log(1+pow(10,(i/10.0)))/log(10)+0.499));
groupCounts.clear();
groupMap.clear();
vector<vector<string> > fastaFileNames;
+ map<string, string> uniqueFastaNames;// so we don't add the same groupfile multiple times
createOligosGroup = false;
+ oligos = new Oligos();
+ numBarcodes = 0; numFPrimers= 0; numLinkers= 0; numSpacers = 0; numRPrimers = 0;
string outputGroupFileName;
map<string, string> variables;
string thisOutputDir = outputDir;
if (outputDir == "") { thisOutputDir = m->hasPath(filesToProcess[l][0][0]); }
variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(filesToProcess[l][0][0]));
variables["[tag]"] = "";
- if(oligosfile != ""){ createOligosGroup = getOligos(fastaFileNames, variables["[filename]"]); }
+ if(oligosfile != ""){ createOligosGroup = getOligos(fastaFileNames, variables["[filename]"], uniqueFastaNames); }
if (createOligosGroup || createFileGroup) {
outputGroupFileName = getOutputFileName("group",variables);
}
m->mothurOut("Making contigs...\n");
createProcesses(filesToProcess[l], outFastaFile, outScrapFastaFile, outMisMatchFile, fastaFileNames, l);
- m->mothurOut("Here...\n");
//remove temp fasta and qual files
for (int i = 0; i < processors; i++) { for(int j = 0; j < filesToProcess[l][i].size(); j++) { m->mothurRemove(filesToProcess[l][i][j]); } }
- if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
+ if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } delete oligos; return 0; }
if(allFiles){
- map<string, string> uniqueFastaNames;// so we don't add the same groupfile multiple times
+ // so we don't add the same groupfile multiple times
map<string, string>::iterator it;
set<string> namesToRemove;
for(int i=0;i<fastaFileNames.size();i++){
if(m->isBlank(fastaFileNames[i][j])){
m->mothurRemove(fastaFileNames[i][j]);
namesToRemove.insert(fastaFileNames[i][j]);
- }else{
- it = uniqueFastaNames.find(fastaFileNames[i][j]);
- if (it == uniqueFastaNames.end()) {
- uniqueFastaNames[fastaFileNames[i][j]] = barcodeNameVector[i];
- }
+ uniqueFastaNames.erase(fastaFileNames[i][j]); //remove from list for group file print
}
}
}
m->openInputFile(it->first, in);
ofstream out;
- string thisGroupName = thisOutputDir + m->getRootName(m->getSimpleName(it->first));
- thisGroupName += getOutputFileName("group",variables); outputNames.push_back(thisGroupName); outputTypes["group"].push_back(thisGroupName);
+ variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(it->first));
+ string thisGroupName = getOutputFileName("group",variables); outputNames.push_back(thisGroupName); outputTypes["group"].push_back(thisGroupName);
m->openOutputFile(thisGroupName, out);
while (!in.eof()){
}
}
m->mothurOut("Done.\n");
+ delete oligos;
}
m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process " + toString(numReads) + " sequences.\n");
vector< vector< vector<string> > > filesToProcess;
if (ffastqfile != "") { //reading one file
- vector< vector<string> > files = readFastqFiles(numReads, ffastqfile, rfastqfile);
+ vector< vector<string> > files = readFastqFiles(numReads, ffastqfile, rfastqfile, findexfile, rindexfile);
//adjust for really large processors or really small files
if (numReads == 0) { m->mothurOut("[ERROR]: no good reads.\n"); m->control_pressed = true; }
if (numReads < processors) {
if (m->control_pressed) { for (int l = 0; l < filesToProcess.size(); l++) { for (int k = 0; k < filesToProcess[l].size(); k++) { for(int j = 0; j < filesToProcess[l][k].size(); j++) { m->mothurRemove(filesToProcess[l][k][j]); } filesToProcess[l][k].clear(); } return filesToProcess; } }
unsigned long int thisFilesReads;
- vector< vector<string> > files = readFastqFiles(thisFilesReads, filePairsToProcess[i][0], filePairsToProcess[i][1]);
+ vector< vector<string> > files = readFastqFiles(thisFilesReads, filePairsToProcess[i][0], filePairsToProcess[i][1], filePairsToProcess[i][2], filePairsToProcess[i][3]);
//adjust for really large processors or really small files
if (thisFilesReads < processors) {
m->mothurOut("[ERROR]: " + filePairsToProcess[i][0] + " has less than " + toString(processors) + " good reads, skipping\n");
for (int k = 0; k < files.size(); k++) { for(int j = 0; j < files[k].size(); j++) { m->mothurRemove(files[k][j]); } files[k].clear(); }
+ //remove from file2Group if necassary
+ map<int, string> cFile2Group;
+ for (map<int, string>::iterator it = file2Group.begin(); it != file2Group.end(); it++) {
+ if ((it->first) < i) { cFile2Group[it->first] = it->second; }
+ else if ((it->first) == i) { } //do nothing, we removed files for i
+ else { cFile2Group[(it->first-1)] = it->second; } //adjust files because i was removed
+ }
+ file2Group = cFile2Group;
}else {
filesToProcess.push_back(files);
numReads += thisFilesReads;
//loop through and create all the processes you want
while (process != processors-1) {
- int pid = fork();
+ pid_t pid = fork();
if (pid > 0) {
processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
for(int i=0;i<tempFASTAFileNames.size();i++){
for(int j=0;j<tempFASTAFileNames[i].size();j++){
if (tempFASTAFileNames[i][j] != "") {
- tempFASTAFileNames[i][j] += toString(getpid()) + ".temp";
+ tempFASTAFileNames[i][j] += m->mothurGetpid(process) + ".temp";
m->openOutputFile(tempFASTAFileNames[i][j], temp); temp.close();
}
}
}
num = driver(files[process],
- outputFasta + toString(getpid()) + ".temp",
- outputScrapFasta + toString(getpid()) + ".temp",
- outputMisMatches + toString(getpid()) + ".temp",
+ outputFasta + m->mothurGetpid(process) + ".temp",
+ outputScrapFasta + m->mothurGetpid(process) + ".temp",
+ outputMisMatches + m->mothurGetpid(process) + ".temp",
tempFASTAFileNames, process, group);
//pass groupCounts to parent
ofstream out;
- string tempFile = toString(getpid()) + ".num.temp";
+ string tempFile = m->mothurGetpid(process) + ".num.temp";
m->openOutputFile(tempFile, out);
out << num << endl;
if (createFileGroup || createOligosGroup) {
}
}
- contigsData* tempcontig = new contigsData(group, files[h], (outputFasta + extension), (outputScrapFasta + extension), (outputMisMatches + extension), align, m, match, misMatch, gapOpen, gapExtend, insert, deltaq, barcodes, primers, tempFASTAFileNames, barcodeNameVector, primerNameVector, pdiffs, bdiffs, tdiffs, createOligosGroup, createFileGroup, allFiles, trimOverlap, h);
+ contigsData* tempcontig = new contigsData(group, files[h], (outputFasta + extension), (outputScrapFasta + extension), (outputMisMatches + extension), align, m, match, misMatch, gapOpen, gapExtend, insert, deltaq, tempFASTAFileNames, oligosfile, reorient, pdiffs, bdiffs, tdiffs, createOligosGroup, createFileGroup, allFiles, trimOverlap, h);
pDataArray.push_back(tempcontig);
hThreadArray[h] = CreateThread(NULL, 0, MyContigsThreadFunction, pDataArray[h], 0, &dwThreadIdArray[h]);
string thisfqualfile = files[1];
string thisrfastafile = files[2];
string thisrqualfile = files[3];
+ string thisfindexfile = files[4];
+ string thisrindexfile = files[5];
- if (m->debug) { m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n"); }
+ if (m->debug) { m->mothurOut("[DEBUG]: ffasta = " + thisffastafile + ".\n[DEBUG]: fqual = " + thisfqualfile + ".\n[DEBUG]: rfasta = " + thisrfastafile + ".\n[DEBUG]: rqual = " + thisrqualfile + ".\n[DEBUG]: findex = " + thisfindexfile + ".\n[DEBUG]: rindex = " + thisrindexfile + ".\n"); }
- ifstream inFFasta, inRFasta, inFQual, inRQual;
+ ifstream inFFasta, inRFasta, inFQual, inRQual, inFIndex, inRIndex;
ofstream outFasta, outMisMatch, outScrapFasta;
m->openInputFile(thisffastafile, inFFasta);
m->openInputFile(thisrfastafile, inRFasta);
m->openInputFile(thisfqualfile, inFQual);
m->openInputFile(thisrqualfile, inRQual);
}
+
+ if (thisfindexfile != "") { m->openInputFile(thisfindexfile, inFIndex); }
+ if (thisrindexfile != "") { m->openInputFile(thisrindexfile, inRIndex); }
+
m->openOutputFile(outputFasta, outFasta);
m->openOutputFile(outputScrapFasta, outScrapFasta);
m->openOutputFile(outputMisMatches, outMisMatch);
outMisMatch << "Name\tLength\tOverlap_Length\tOverlap_Start\tOverlap_End\tMisMatches\tNum_Ns\n";
- TrimOligos trimOligos(pdiffs, bdiffs, 0, 0, primers, barcodes);
+ TrimOligos trimOligos(pdiffs, bdiffs, 0, 0, oligos->getPairedPrimers(), oligos->getPairedBarcodes());
+
+ TrimOligos* rtrimOligos = NULL;
+ if (reorient) {
+ rtrimOligos = new TrimOligos(pdiffs, bdiffs, 0, 0, oligos->getReorientedPairedPrimers(), oligos->getReorientedPairedBarcodes()); numBarcodes = oligos->getReorientedPairedBarcodes().size();
+ }
while ((!inFFasta.eof()) && (!inRFasta.eof())) {
fQual = new QualityScores(inFQual); m->gobble(inFQual);
rQual = new QualityScores(inRQual); m->gobble(inRQual);
}
+ Sequence findexBarcode("findex", "NONE"); Sequence rindexBarcode("rindex", "NONE");
+ if (thisfindexfile != "") {
+ Sequence temp(inFIndex); m->gobble(inFIndex);
+ findexBarcode.setAligned(temp.getAligned());
+ }
+
+ if (thisrindexfile != "") {
+ Sequence temp(inRIndex); m->gobble(inRIndex);
+ rindexBarcode.setAligned(temp.getAligned());
+ }
int barcodeIndex = 0;
int primerIndex = 0;
-
- if(barcodes.size() != 0){
+ Sequence savedFSeq(fSeq.getName(), fSeq.getAligned()); Sequence savedRSeq(rSeq.getName(), rSeq.getAligned());
+ Sequence savedFindex(findexBarcode.getName(), findexBarcode.getAligned()); Sequence savedRIndex(rindexBarcode.getName(), rindexBarcode.getAligned());
+ QualityScores* savedFQual = NULL; QualityScores* savedRQual = NULL;
+ if (thisfqualfile != "") {
+ savedFQual = new QualityScores(fQual->getName(), fQual->getQualityScores());
+ savedRQual = new QualityScores(rQual->getName(), rQual->getQualityScores());
+ }
+
+ if(numBarcodes != 0){
if (thisfqualfile != "") {
- success = trimOligos.stripBarcode(fSeq, rSeq, *fQual, *rQual, barcodeIndex);
+ if ((thisfindexfile != "") || (thisrindexfile != "")) {
+ success = trimOligos.stripBarcode(findexBarcode, rindexBarcode, *fQual, *rQual, barcodeIndex);
+ }else {
+ success = trimOligos.stripBarcode(fSeq, rSeq, *fQual, *rQual, barcodeIndex);
+ }
}else {
success = trimOligos.stripBarcode(fSeq, rSeq, barcodeIndex);
}
else{ currentSeqsDiffs += success; }
}
- if(primers.size() != 0){
+ if(numFPrimers != 0){
if (thisfqualfile != "") {
success = trimOligos.stripForward(fSeq, rSeq, *fQual, *rQual, primerIndex);
}else {
if (currentSeqsDiffs > tdiffs) { trashCode += 't'; }
+ if (reorient && (trashCode != "")) { //if you failed and want to check the reverse
+ int thisSuccess = 0;
+ string thisTrashCode = "";
+ int thisCurrentSeqsDiffs = 0;
+
+ int thisBarcodeIndex = 0;
+ int thisPrimerIndex = 0;
+
+ if(numBarcodes != 0){
+ if (thisfqualfile != "") {
+ if ((thisfindexfile != "") || (thisrindexfile != "")) {
+ thisSuccess = rtrimOligos->stripBarcode(savedFindex, savedRIndex, *savedFQual, *savedRQual, thisBarcodeIndex);
+ }else {
+ thisSuccess = rtrimOligos->stripBarcode(savedFSeq, savedRSeq, *savedFQual, *savedRQual, thisBarcodeIndex);
+ }
+ }else {
+ thisSuccess = rtrimOligos->stripBarcode(savedFSeq, savedRSeq, thisBarcodeIndex);
+ }
+ if(thisSuccess > bdiffs) { thisTrashCode += 'b'; }
+ else{ thisCurrentSeqsDiffs += thisSuccess; }
+ }
+
+ if(numFPrimers != 0){
+ if (thisfqualfile != "") {
+ thisSuccess = rtrimOligos->stripForward(savedFSeq, savedRSeq, *savedFQual, *savedRQual, thisPrimerIndex);
+ }else {
+ thisSuccess = rtrimOligos->stripForward(savedFSeq, savedRSeq, thisPrimerIndex);
+ }
+ if(thisSuccess > pdiffs) { thisTrashCode += 'f'; }
+ else{ thisCurrentSeqsDiffs += thisSuccess; }
+ }
+
+ if (thisCurrentSeqsDiffs > tdiffs) { thisTrashCode += 't'; }
+
+ if (thisTrashCode == "") {
+ trashCode = thisTrashCode;
+ success = thisSuccess;
+ currentSeqsDiffs = thisCurrentSeqsDiffs;
+ barcodeIndex = thisBarcodeIndex;
+ primerIndex = thisPrimerIndex;
+ savedFSeq.reverseComplement();
+ savedRSeq.reverseComplement();
+ fSeq.setAligned(savedFSeq.getAligned());
+ rSeq.setAligned(savedRSeq.getAligned());
+ if(thisfqualfile != ""){
+ savedFQual->flipQScores(); savedRQual->flipQScores();
+ fQual->setScores(savedFQual->getScores()); rQual->setScores(savedRQual->getScores());
+ }
+ }else { trashCode += "(" + thisTrashCode + ")"; }
+ }
+
+
//flip the reverse reads
rSeq.reverseComplement();
if (thisfqualfile != "") { rQual->flipQScores(); }
if (thisfqualfile != "") {
scores1 = fQual->getQualityScores();
scores2 = rQual->getQualityScores();
- delete fQual; delete rQual;
+ delete fQual; delete rQual; delete savedFQual; delete savedRQual;
}
// if (num < 5) { cout << fSeq.getStartPos() << '\t' << fSeq.getEndPos() << '\t' << rSeq.getStartPos() << '\t' << rSeq.getEndPos() << endl; }
if (seq2End < overlapEnd) { overlapEnd = seq2End; } //smallest end position is where overlapping ends
int oStart = contig.length();
+ //cout << fSeq.getAligned() << endl; cout << rSeq.getAligned() << endl;
for (int i = overlapStart; i < overlapEnd; i++) {
+ //cout << seq1[i] << ' ' << seq2[i] << ' ' << scores1[ABaseMap[i]] << ' ' << scores2[BBaseMap[i]] << endl;
if (seq1[i] == seq2[i]) { //match, add base and choose highest score
contig += seq1[i];
}else if (((seq1[i] == '.') || (seq1[i] == '-')) && ((seq2[i] != '-') && (seq2[i] != '.'))) { //seq1 is a gap and seq2 is a base, choose seq2, unless quality score for base is below insert. In that case eliminate base
}else { //seq2 ends before seq1 so take from overlap to length from seq1
for (int i = overlapEnd; i < length; i++) { contig += seq1[i]; }
}
-
+ //cout << contig << endl;
+ //exit(1);
if (trimOverlap) { contig = contig.substr(overlapStart-1, oend-oStart); if (contig.length() == 0) { trashCode += "l"; } }
if(trashCode.length() == 0){
if (m->debug) { m->mothurOut(fSeq.getName()); }
if (createOligosGroup) {
- if(barcodes.size() != 0){
- string thisGroup = barcodeNameVector[barcodeIndex];
- if (primers.size() != 0) {
- if (primerNameVector[primerIndex] != "") {
- if(thisGroup != "") {
- thisGroup += "." + primerNameVector[primerIndex];
- }else {
- thisGroup = primerNameVector[primerIndex];
- }
- }
- }
-
+ string thisGroup = oligos->getGroupName(barcodeIndex, primerIndex);
if (m->debug) { m->mothurOut(", group= " + thisGroup + "\n"); }
int pos = thisGroup.find("ignore");
if (it == groupCounts.end()) { groupCounts[thisGroup] = 1; }
else { groupCounts[it->first] ++; }
}else { ignore = true; }
-
- }
}else if (createFileGroup) {
int pos = group.find("ignore");
if (pos == string::npos) {
}else { ignore = true; }
}
if (m->debug) { m->mothurOut("\n"); }
-
+
if(allFiles && !ignore){
ofstream output;
m->openOutputFileAppend(fastaFileNames[barcodeIndex][primerIndex], output);
num++;
//report progress
- if((num) % 1000 == 0){ m->mothurOut(toString(num)); m->mothurOutEndLine(); }
+ if((num) % 1000 == 0){ m->mothurOutJustToScreen(toString(num)); m->mothurOutEndLine(); }
}
//report progress
inRQual.close();
}
delete alignment;
+ if (reorient) { delete rtrimOligos; }
if (m->control_pressed) { m->mothurRemove(outputFasta); m->mothurRemove(outputScrapFasta);m->mothurRemove(outputMisMatches); }
}
}
//**********************************************************************************************************************
-vector< vector<string> > MakeContigsCommand::readFastqFiles(unsigned long int& count, string ffastq, string rfastq){
+vector< vector<string> > MakeContigsCommand::readFastqFiles(unsigned long int& count, string ffastq, string rfastq, string findex, string rindex){
try {
vector< vector<string> > files;
//maps processors number to file pointer
- map<int, vector<ofstream*> > tempfiles; //tempfiles[0] = forwardFasta, [1] = forwardQual, [2] = reverseFasta, [3] = reverseQual
+ map<int, vector<ofstream*> > tempfiles; //tempfiles[0] = forwardFasta, [1] = forwardQual, [2] = reverseFasta, [3] = reverseQual, tempfiles[4] = forwardIndex, [4] = forwardReverse
map<int, vector<ofstream*> >::iterator it;
//create files to write to
ofstream* outFQ = new ofstream; temp.push_back(outFQ);
ofstream* outRF = new ofstream; temp.push_back(outRF);
ofstream* outRQ = new ofstream; temp.push_back(outRQ);
+ ofstream* outFI = new ofstream; temp.push_back(outFI);
+ ofstream* outRI = new ofstream; temp.push_back(outRI);
tempfiles[i] = temp;
vector<string> names;
string rfastafilename = thisOutputDir + m->getRootName(m->getSimpleName(rfastq)) + toString(i) + "rfastatemp";
string fqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(ffastq)) + toString(i) + "fqualtemp";
string rqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(rfastq)) + toString(i) + "rqualtemp";
+ string findexfilename = ""; string rindexfilename = "";
+ noneOk = false; //flag to oligos file read that its okay to allow for non paired barcodes
+ if (findex != "") { findexfilename = thisOutputDir + m->getRootName(m->getSimpleName(findex)) + toString(i) + "findextemp"; m->openOutputFile(findexfilename, *outFI); noneOk = true; }
+ if (rindex != "") { rindexfilename = thisOutputDir + m->getRootName(m->getSimpleName(rindex)) + toString(i) + "rindextemp"; m->openOutputFile(rindexfilename, *outRI); noneOk = true; }
names.push_back(ffastafilename); names.push_back(fqualfilename);
names.push_back(rfastafilename); names.push_back(rqualfilename);
+ names.push_back(findexfilename); names.push_back(rindexfilename);
files.push_back(names);
m->openOutputFile(ffastafilename, *outFF);
for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
//remove files
for (int i = 0; i < files.size(); i++) {
- for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); }
+ for(int j = 0; j < files[i].size(); j++) { if (files[i][j] != "") { m->mothurRemove(files[i][j]); } }
}
}
ifstream inReverse;
m->openInputFile(rfastq, inReverse);
+ ifstream infIndex, inrIndex;
+ bool findexIsGood = false;
+ bool rindexIsGood = false;
+ if (findex != "") { m->openInputFile(findex, infIndex); findexIsGood = true; }
+ if (rindex != "") { m->openInputFile(rindex, inrIndex); rindexIsGood = true; }
+
count = 0;
map<string, fastqRead> uniques;
+ map<string, fastqRead> iUniques;
+ map<string, pairFastqRead> pairUniques;
map<string, fastqRead>::iterator itUniques;
- while ((!inForward.eof()) || (!inReverse.eof())) {
+ while ((!inForward.eof()) || (!inReverse.eof()) || (findexIsGood) || (rindexIsGood)) {
- if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inForward.close(); inReverse.close(); return files; }
+ if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { if (files[i][j] != "") { m->mothurRemove(files[i][j]); } } } inForward.close(); inReverse.close(); if (findex != "") { infIndex.close(); } if (findex != "") { inrIndex.close(); } return files; }
//get a read from forward and reverse fastq files
- bool ignoref, ignorer;
- fastqRead thisFread, thisRread;
+ bool ignoref, ignorer, ignorefi, ignoreri;
+ fastqRead thisFread, thisRread, thisFIread, thisRIread;
if (!inForward.eof()) { thisFread = readFastq(inForward, ignoref); }
else { ignoref = true; }
if (!inReverse.eof()) { thisRread = readFastq(inReverse, ignorer); }
else { ignorer = true; }
+ if (findexIsGood) { thisFIread = readFastq(infIndex, ignorefi); if (infIndex.eof()) { findexIsGood = false; } }
+ else { ignorefi = true; }
+ if (rindexIsGood) { thisRIread = readFastq(inrIndex, ignoreri); if (inrIndex.eof()) { rindexIsGood = false; } }
+ else { ignoreri = true; }
+
+ bool allowOne = false;
+ if ((findex == "") || (rindex == "")) { allowOne = true; }
+ vector<pairFastqRead> frReads = getReads(ignoref, ignorer, thisFread, thisRread, uniques, false);
+ vector<pairFastqRead> friReads = getReads(ignorefi, ignoreri, thisFIread, thisRIread, iUniques, allowOne);
+
+ //add in index info if provided
+ vector<pairFastqRead> reads = frReads;
+ if ((findex != "") || (rindex != "")) { reads = mergeReads(frReads, friReads, pairUniques); }
- vector<pairFastqRead> reads = getReads(ignoref, ignorer, thisFread, thisRread, uniques);
-
for (int i = 0; i < reads.size(); i++) {
fastqRead fread = reads[i].forward;
fastqRead rread = reads[i].reverse;
+ fastqRead firead = reads[i].findex;
+ fastqRead riread = reads[i].rindex;
- if (m->debug) { m->mothurOut(toString(count) + '\t' + fread.name + '\t' + rread.name + '\n'); }
+ if (m->debug) { m->mothurOut(toString(count) + '\t' + fread.name + '\t' + rread.name + '\n'); if (findex != "") { m->mothurOut(toString(count) + '\t' + firead.name + '\n'); } if (rindex != "") { m->mothurOut(toString(count) + '\t' + riread.name + '\n'); } }
//if (checkReads(fread, rread, ffastq, rfastq)) {
- if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { m->mothurRemove(files[i][j]); } } inForward.close(); inReverse.close(); return files; }
+ if (m->control_pressed) { for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } } for (int i = 0; i < files.size(); i++) { for(int j = 0; j < files[i].size(); j++) { if (files[i][j] != "") { m->mothurRemove(files[i][j]); } } } inForward.close(); inReverse.close(); if (findex != "") { infIndex.close(); } if (findex != "") { inrIndex.close(); } return files; }
//if the reads are okay write to output files
int process = count % processors;
*(tempfiles[process][3]) << ">" << rread.name << endl;
for (int i = 0; i < rread.scores.size(); i++) { *(tempfiles[process][3]) << rread.scores[i] << " "; }
*(tempfiles[process][3]) << endl;
-
+ if (findex != "") { *(tempfiles[process][4]) << ">" << firead.name << endl << firead.sequence << endl; }
+ if (rindex != "") { *(tempfiles[process][5]) << ">" << riread.name << endl << riread.sequence << endl; }
+
count++;
//report progress
- if((count) % 10000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
+ if((count) % 10000 == 0){ m->mothurOutJustToScreen(toString(count)); m->mothurOutEndLine(); }
//}
}
}
if (uniques.size() != 0) {
for (itUniques = uniques.begin(); itUniques != uniques.end(); itUniques++) {
+ if (m->control_pressed) { break; }
m->mothurOut("[WARNING]: did not find paired read for " + itUniques->first + ", ignoring.\n");
}
+ for (map<string, pairFastqRead>:: iterator it = pairUniques.begin(); it != pairUniques.end(); it++) {
+ if (m->control_pressed) { break; }
+ m->mothurOut("[WARNING]: did not find paired read for " + (it->first).substr(1) + ", ignoring.\n");
+ }
m->mothurOutEndLine();
}
for (it = tempfiles.begin(); it!=tempfiles.end(); it++) { for (int i = 0; i < (it->second).size(); i++) { (*(it->second)[i]).close(); delete (it->second)[i]; } }
inForward.close();
inReverse.close();
-
+ if (findex != "") { infIndex.close(); }
+ if (rindex != "") { inrIndex.close(); }
+
return files;
}
catch(exception& e) {
if (fqualfile != "") { fqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(fqualfile)) + toString(i) + "fqual.temp"; m->openOutputFile(fqualfilename, *outFQ); }
string rqualfilename = "";
if (rqualfile != "") { rqualfilename = thisOutputDir + m->getRootName(m->getSimpleName(rqualfile)) + toString(i) + "rqual.temp"; m->openOutputFile(rqualfilename, *outRQ); }
+ string findexfilename = ""; string rindexfilename = "";
names.push_back(ffastafilename); names.push_back(fqualfilename);
names.push_back(rfastafilename); names.push_back(rqualfilename);
+ names.push_back(findexfilename); names.push_back(rindexfilename);
files.push_back(names);
m->openOutputFile(ffastafilename, *outFF);
}else { ignorer = true; }
}
- vector<pairFastqRead> reads = getReads(ignoref, ignorer, thisFread, thisRread, uniques);
+ vector<pairFastqRead> reads = getReads(ignoref, ignorer, thisFread, thisRread, uniques, false);
for (int i = 0; i < reads.size(); i++) {
fastqRead fread = reads[i].forward;
}
}
//**********************************************************************************************************************
-vector<pairFastqRead> MakeContigsCommand::getReads(bool ignoref, bool ignorer, fastqRead forward, fastqRead reverse, map<string, fastqRead>& uniques){
+vector<pairFastqRead> MakeContigsCommand::getReads(bool ignoref, bool ignorer, fastqRead forward, fastqRead reverse, map<string, fastqRead>& uniques, bool allowOne){
try {
vector<pairFastqRead> reads;
map<string, fastqRead>::iterator itUniques;
}
}else if (!ignoref && ignorer) { //ignore reverse keep forward
- //look for forward pair
- itUniques = uniques.find(forward.name);
- if (itUniques != uniques.end()) { //we have the pair for this read
- pairFastqRead temp(forward, itUniques->second);
+ if (allowOne) {
+ fastqRead dummy;
+ pairFastqRead temp(forward, dummy);
reads.push_back(temp);
- uniques.erase(itUniques);
- }else { //save this read for later
- uniques[forward.name] = forward;
+ }else {
+ //look for forward pair
+ itUniques = uniques.find(forward.name);
+ if (itUniques != uniques.end()) { //we have the pair for this read
+ pairFastqRead temp(forward, itUniques->second);
+ reads.push_back(temp);
+ uniques.erase(itUniques);
+ }else { //save this read for later
+ uniques[forward.name] = forward;
+ }
}
-
}else if (ignoref && !ignorer) { //ignore forward keep reverse
- //look for reverse pair
- itUniques = uniques.find(reverse.name);
- if (itUniques != uniques.end()) { //we have the pair for this read
- pairFastqRead temp(itUniques->second, reverse);
+ if (allowOne) {
+ fastqRead dummy;
+ pairFastqRead temp(dummy, reverse);
reads.push_back(temp);
- uniques.erase(itUniques);
- }else { //save this read for later
- uniques[reverse.name] = reverse;
+ }else {
+ //look for reverse pair
+ itUniques = uniques.find(reverse.name);
+ if (itUniques != uniques.end()) { //we have the pair for this read
+ pairFastqRead temp(itUniques->second, reverse);
+ reads.push_back(temp);
+ uniques.erase(itUniques);
+ }else { //save this read for later
+ uniques[reverse.name] = reverse;
+ }
}
}//else ignore both and do nothing
}
}
//**********************************************************************************************************************
+//look through the reads from the forward and reverse files and try to find matching reads from index files.
+vector<pairFastqRead> MakeContigsCommand::mergeReads(vector<pairFastqRead> thisReads, vector<pairFastqRead> indexes, map<string, pairFastqRead>& uniques){
+ try {
+ vector<pairFastqRead> reads;
+ map<string, pairFastqRead>::iterator itUniques;
+
+ set<int> foundIndexes;
+ for (int i = 0; i < thisReads.size(); i++) {
+ bool found = false;
+ for (int j = 0; j < indexes.size(); j++) {
+
+ //incase only one index
+ string indexName = indexes[j].forward.name;
+ if (indexName == "") { indexName = indexes[j].reverse.name; }
+
+ if (thisReads[i].forward.name == indexName){
+ thisReads[i].findex = indexes[j].forward;
+ thisReads[i].rindex = indexes[j].reverse;
+ reads.push_back(thisReads[i]);
+ found = true;
+ foundIndexes.insert(j);
+ }
+ }
+
+ if (!found) {
+ //look for forward pair
+ itUniques = uniques.find('i'+thisReads[i].forward.name);
+ if (itUniques != uniques.end()) { //we have the pair for this read
+ thisReads[i].findex = itUniques->second.forward;
+ thisReads[i].rindex = itUniques->second.reverse;
+ reads.push_back(thisReads[i]);
+ uniques.erase(itUniques);
+ }else { //save this read for later
+ uniques['r'+thisReads[i].forward.name] = thisReads[i];
+ }
+ }
+ }
+
+ if (foundIndexes.size() != indexes.size()) { //if we didnt match all the indexes look for them in uniques
+ for (int j = 0; j < indexes.size(); j++) {
+ if (foundIndexes.count(j) == 0) { //we didnt find this one
+ //incase only one index
+ string indexName = indexes[j].forward.name;
+ if (indexName == "") { indexName = indexes[j].reverse.name; }
+
+ //look for forward pair
+ itUniques = uniques.find('r'+indexName);
+ if (itUniques != uniques.end()) { //we have the pair for this read
+ pairFastqRead temp(itUniques->second.forward, itUniques->second.reverse, indexes[j].forward, indexes[j].reverse);
+ reads.push_back(temp);
+ uniques.erase(itUniques);
+ }else { //save this read for later
+ uniques['i'+indexName] = indexes[j];
+ }
+ }
+ }
+ }
+
+
+ return reads;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "MakeContigsCommand", "mergeReads");
+ exit(1);
+ }
+}
+//**********************************************************************************************************************
fastqRead MakeContigsCommand::readFastq(ifstream& in, bool& ignore){
try {
fastqRead read;
vector<int> qualScores = convertQual(quality);
+ m->checkName(name);
read.name = name;
read.sequence = sequence;
read.scores = qualScores;
+
+ if (m->debug) { m->mothurOut("[DEBUG]: " + read.name + " " + read.sequence + " " + quality + "\n"); }
return read;
}
}
}*/
//***************************************************************************************************************
+//lines can be 2, 3, or 4 columns
+// forward.fastq reverse.fastq -> 2 column
+// groupName forward.fastq reverse.fastq -> 3 column
+// forward.fastq reverse.fastq forward.index.fastq reverse.index.fastq -> 4 column
+// forward.fastq reverse.fastq none reverse.index.fastq -> 4 column
+// forward.fastq reverse.fastq forward.index.fastq none -> 4 column
vector< vector<string> > MakeContigsCommand::readFileNames(string filename){
try {
vector< vector<string> > files;
- string forward, reverse;
+ string forward, reverse, findex, rindex;
ifstream in;
m->openInputFile(filename, in);
if (m->control_pressed) { return files; }
- in >> forward; m->gobble(in);
- in >> reverse;
+ string line = m->getline(in); m->gobble(in);
+ vector<string> pieces = m->splitWhiteSpace(line);
string group = "";
- while (!in.eof()) { //do we have a group assigned to this pair
- char c = in.get();
- if (c == 10 || c == 13 || c == -1){ break; }
- else if (c == 32 || c == 9){;} //space or tab
- else { group += c; }
- }
-
- if (group != "") {
- //line in file look like: group forward reverse
- string temp = forward;
- forward = reverse;
- reverse = group;
- group = temp;
+ if (pieces.size() == 2) {
+ forward = pieces[0];
+ reverse = pieces[1];
+ group = "";
+ findex = "";
+ rindex = "";
+ }else if (pieces.size() == 3) {
+ group = pieces[0];
+ forward = pieces[1];
+ reverse = pieces[2];
+ findex = "";
+ rindex = "";
createFileGroup = true;
+ }else if (pieces.size() == 4) {
+ forward = pieces[0];
+ reverse = pieces[1];
+ findex = pieces[2];
+ rindex = pieces[3];
+ if ((findex == "none") || (findex == "NONE")){ findex = ""; }
+ if ((rindex == "none") || (rindex == "NONE")){ rindex = ""; }
+ }else {
+ m->mothurOut("[ERROR]: file lines can be 2, 3, or 4 columns. The forward fastq files in the first column and their matching reverse fastq files in the second column, or a groupName then forward fastq file and reverse fastq file, or forward fastq file then reverse fastq then forward index and reverse index file. If you only have one index file add 'none' for the other one. \n"); m->control_pressed = true;
}
- m->gobble(in);
+ if (m->debug) { m->mothurOut("[DEBUG]: group = " + group + ", forward = " + forward + ", reverse = " + reverse + ", forwardIndex = " + findex + ", reverseIndex = " + rindex + ".\n"); }
+
+ if (inputDir != "") {
+ string path = m->hasPath(forward);
+ if (path == "") { forward = inputDir + forward; }
+
+ path = m->hasPath(reverse);
+ if (path == "") { reverse = inputDir + reverse; }
+
+ if (findex != "") {
+ path = m->hasPath(findex);
+ if (path == "") { findex = inputDir + findex; }
+ }
+
+ if (rindex != "") {
+ path = m->hasPath(rindex);
+ if (path == "") { rindex = inputDir + rindex; }
+ }
+ }
//check to make sure both are able to be opened
ifstream in2;
m->mothurOut("[WARNING]: can't find " + reverse + ", ignoring pair.\n");
}else{ in3.close(); }
- if ((openForward != 1) && (openReverse != 1)) { //good pair
+ int openFindex = 0;
+ if (findex != "") {
+ ifstream in4;
+ openFindex = m->openInputFile(findex, in4, "noerror"); in4.close();
+
+ //if you can't open it, try default location
+ if (openFindex == 1) {
+ if (m->getDefaultPath() != "") { //default path is set
+ string tryPath = m->getDefaultPath() + m->getSimpleName(findex);
+ m->mothurOut("Unable to open " + findex + ". Trying default " + tryPath); m->mothurOutEndLine();
+ ifstream in5;
+ openFindex = m->openInputFile(tryPath, in5, "noerror");
+ in5.close();
+ findex = tryPath;
+ }
+ }
+
+ //if you can't open it, try output location
+ if (openFindex == 1) {
+ if (m->getOutputDir() != "") { //default path is set
+ string tryPath = m->getOutputDir() + m->getSimpleName(findex);
+ m->mothurOut("Unable to open " + findex + ". Trying output directory " + tryPath); m->mothurOutEndLine();
+ ifstream in6;
+ openFindex = m->openInputFile(tryPath, in6, "noerror");
+ findex = tryPath;
+ in6.close();
+ }
+ }
+
+ if (openFindex == 1) { //can't find it
+ m->mothurOut("[WARNING]: can't find " + findex + ", ignoring pair.\n");
+ }
+ }
+
+ int openRindex = 0;
+ if (rindex != "") {
+ ifstream in7;
+ openRindex = m->openInputFile(rindex, in7, "noerror"); in7.close();
+
+ //if you can't open it, try default location
+ if (openRindex == 1) {
+ if (m->getDefaultPath() != "") { //default path is set
+ string tryPath = m->getDefaultPath() + m->getSimpleName(rindex);
+ m->mothurOut("Unable to open " + rindex + ". Trying default " + tryPath); m->mothurOutEndLine();
+ ifstream in8;
+ openRindex = m->openInputFile(tryPath, in8, "noerror");
+ in8.close();
+ rindex = tryPath;
+ }
+ }
+
+ //if you can't open it, try output location
+ if (openRindex == 1) {
+ if (m->getOutputDir() != "") { //default path is set
+ string tryPath = m->getOutputDir() + m->getSimpleName(rindex);
+ m->mothurOut("Unable to open " + rindex + ". Trying output directory " + tryPath); m->mothurOutEndLine();
+ ifstream in9;
+ openRindex = m->openInputFile(tryPath, in9, "noerror");
+ rindex = tryPath;
+ in9.close();
+ }
+ }
+
+ if (openRindex == 1) { //can't find it
+ m->mothurOut("[WARNING]: can't find " + rindex + ", ignoring pair.\n");
+ }
+ }
+
+
+ if ((openForward != 1) && (openReverse != 1) && (openFindex != 1) && (openRindex != 1)) { //good pair
file2Group[files.size()] = group;
vector<string> pair;
pair.push_back(forward);
pair.push_back(reverse);
+ pair.push_back(findex);
+ pair.push_back(rindex);
+ if (((findex != "") || (rindex != "")) && (oligosfile == "")) { m->mothurOut("[ERROR]: You need to provide an oligos file if you are going to use an index file.\n"); m->control_pressed = true; }
files.push_back(pair);
}
}
return files;
}
catch(exception& e) {
- m->errorOut(e, "MakeContigsCommand", "checkReads");
+ m->errorOut(e, "MakeContigsCommand", "readFileNames");
exit(1);
}
}
//illumina data requires paired forward and reverse data
//BARCODE atgcatgc atgcatgc groupName
//PRIMER atgcatgc atgcatgc groupName
-//PRIMER atgcatgc atgcatgc
-bool MakeContigsCommand::getOligos(vector<vector<string> >& fastaFileNames, string rootname){
+//PRIMER atgcatgc atgcatgc
+bool MakeContigsCommand::getOligos(vector<vector<string> >& fastaFileNames, string rootname, map<string, string>& fastaFile2Group){
try {
- ifstream in;
- m->openInputFile(oligosfile, in);
-
- ofstream test;
-
- string type, foligo, roligo, group;
+ if (m->debug) { m->mothurOut("[DEBUG]: oligosfile = " + oligosfile + "\n"); }
- int indexPrimer = 0;
- int indexBarcode = 0;
- set<string> uniquePrimers;
- set<string> uniqueBarcodes;
-
- while(!in.eof()){
-
- in >> type;
+ bool allBlank = false;
+ oligos->read(oligosfile, false);
+
+ if (m->control_pressed) { return false; } //error in reading oligos
+
+ if (oligos->hasPairedBarcodes()) {
+ numFPrimers = oligos->getPairedPrimers().size();
+ numBarcodes = oligos->getPairedBarcodes().size();
+ }else {
+ m->mothurOut("[ERROR]: make.contigs requires paired barcodes and primers. You can set one end to NONE if you are using an index file.\n"); m->control_pressed = true;
+ }
- if (m->debug) { m->mothurOut("[DEBUG]: reading type - " + type + ".\n"); }
-
- if(type[0] == '#'){
- while (!in.eof()) { char c = in.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there
- m->gobble(in);
- }
- else{
- m->gobble(in);
- //make type case insensitive
- for(int i=0;i<type.length();i++){ type[i] = toupper(type[i]); }
-
- in >> foligo;
-
- if (m->debug) { m->mothurOut("[DEBUG]: reading - " + foligo + ".\n"); }
-
- for(int i=0;i<foligo.length();i++){
- foligo[i] = toupper(foligo[i]);
- if(foligo[i] == 'U') { foligo[i] = 'T'; }
- }
-
- if(type == "PRIMER"){
- m->gobble(in);
-
- in >> roligo;
-
- for(int i=0;i<roligo.length();i++){
- roligo[i] = toupper(roligo[i]);
- if(roligo[i] == 'U') { roligo[i] = 'T'; }
- }
- //roligo = reverseOligo(roligo);
-
- group = "";
-
- // get rest of line in case there is a primer name
- while (!in.eof()) {
- char c = in.get();
- if (c == 10 || c == 13){ break; }
- else if (c == 32 || c == 9){;} //space or tab
- else { group += c; }
- }
-
- oligosPair newPrimer(foligo, roligo);
-
- //check for repeat barcodes
- string tempPair = foligo+roligo;
- if (uniquePrimers.count(tempPair) != 0) { m->mothurOut("primer pair " + newPrimer.forward + " " + newPrimer.reverse + " is in your oligos file already."); m->mothurOutEndLine(); }
- else { uniquePrimers.insert(tempPair); }
-
- if (m->debug) { if (group != "") { m->mothurOut("[DEBUG]: reading group " + group + ".\n"); }else{ m->mothurOut("[DEBUG]: no group for primer pair " + newPrimer.forward + " " + newPrimer.reverse + ".\n"); } }
-
- primers[indexPrimer]=newPrimer; indexPrimer++;
- primerNameVector.push_back(group);
- }else if(type == "BARCODE"){
- m->gobble(in);
-
- in >> roligo;
-
- for(int i=0;i<roligo.length();i++){
- roligo[i] = toupper(roligo[i]);
- if(roligo[i] == 'U') { roligo[i] = 'T'; }
- }
- //roligo = reverseOligo(roligo);
-
- oligosPair newPair(foligo, roligo);
-
- group = "";
- while (!in.eof()) {
- char c = in.get();
- if (c == 10 || c == 13){ break; }
- else if (c == 32 || c == 9){;} //space or tab
- else { group += c; }
- }
-
- if (m->debug) { m->mothurOut("[DEBUG]: barcode pair " + newPair.forward + " " + newPair.reverse + ", and group = " + group + ".\n"); }
-
- //check for repeat barcodes
- string tempPair = foligo+roligo;
- if (uniqueBarcodes.count(tempPair) != 0) { m->mothurOut("barcode pair " + newPair.forward + " " + newPair.reverse + " is in your oligos file already, disregarding."); m->mothurOutEndLine(); }
- else { uniqueBarcodes.insert(tempPair); }
-
- barcodes[indexBarcode]=newPair; indexBarcode++;
- barcodeNameVector.push_back(group);
- }else if(type == "LINKER"){
- linker.push_back(foligo);
- m->mothurOut("[WARNING]: make.contigs is not setup to remove linkers, ignoring.\n");
- }else if(type == "SPACER"){
- spacer.push_back(foligo);
- m->mothurOut("[WARNING]: make.contigs is not setup to remove spacers, ignoring.\n");
- }
- else{ m->mothurOut("[WARNING]: " + type + " is not recognized as a valid type. Choices are primer, barcode, linker and spacer. Ignoring " + foligo + "."); m->mothurOutEndLine(); }
- }
- m->gobble(in);
- }
- in.close();
-
- if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0; }
-
- //add in potential combos
- if(barcodeNameVector.size() == 0){
- oligosPair temp("", "");
- barcodes[0] = temp;
- barcodeNameVector.push_back("");
- }
-
- if(primerNameVector.size() == 0){
- oligosPair temp("", "");
- primers[0] = temp;
- primerNameVector.push_back("");
- }
-
- fastaFileNames.resize(barcodeNameVector.size());
+ if (m->control_pressed) { return false; }
+
+ numLinkers = oligos->getLinkers().size();
+ numSpacers = oligos->getSpacers().size();
+ numRPrimers = oligos->getReversePrimers().size();
+ if (numLinkers != 0) { m->mothurOut("[WARNING]: make.contigs is not setup to remove linkers, ignoring.\n"); }
+ if (numSpacers != 0) { m->mothurOut("[WARNING]: make.contigs is not setup to remove spacers, ignoring.\n"); }
+
+ vector<string> groupNames = oligos->getGroupNames();
+ if (groupNames.size() == 0) { allFiles = 0; allBlank = true; }
+
+
+ fastaFileNames.resize(oligos->getBarcodeNames().size());
for(int i=0;i<fastaFileNames.size();i++){
- fastaFileNames[i].assign(primerNameVector.size(), "");
+ for(int j=0;j<oligos->getPrimerNames().size();j++){ fastaFileNames[i].push_back(""); }
}
-
- if(allFiles){
- set<string> uniqueNames; //used to cleanup outputFileNames
- for(map<int, oligosPair>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
- for(map<int, oligosPair>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
-
- string primerName = primerNameVector[itPrimer->first];
- string barcodeName = barcodeNameVector[itBar->first];
+
+ if (allFiles) {
+ set<string> uniqueNames; //used to cleanup outputFileNames
+ map<int, oligosPair> barcodes = oligos->getPairedBarcodes();
+ map<int, oligosPair> primers = oligos->getPairedPrimers();
+ for(map<int, oligosPair>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
+ for(map<int, oligosPair>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
- if ((primerName == "ignore") || (barcodeName == "ignore")) { } //do nothing
- else {
+ string primerName = oligos->getPrimerName(itPrimer->first);
+ string barcodeName = oligos->getBarcodeName(itBar->first);
+
+ if ((primerName == "ignore") || (barcodeName == "ignore")) { } //do nothing
+ else if ((primerName == "") && (barcodeName == "")) { } //do nothing
+ else {
string comboGroupName = "";
string fastaFileName = "";
string qualFileName = "";
string countFileName = "";
if(primerName == ""){
- comboGroupName = barcodeNameVector[itBar->first];
- }
- else{
+ comboGroupName = barcodeName;
+ }else{
if(barcodeName == ""){
- comboGroupName = primerNameVector[itPrimer->first];
+ comboGroupName = primerName;
}
else{
- comboGroupName = barcodeNameVector[itBar->first] + "." + primerNameVector[itPrimer->first];
+ comboGroupName = barcodeName + "." + primerName;
}
}
ofstream temp;
- fastaFileName = rootname + comboGroupName + ".fasta";
+ map<string, string> variables;
+ variables["[filename]"] = rootname;
+ variables["[tag]"] = comboGroupName;
+ fastaFileName = getOutputFileName("fasta", variables);
if (uniqueNames.count(fastaFileName) == 0) {
outputNames.push_back(fastaFileName);
outputTypes["fasta"].push_back(fastaFileName);
uniqueNames.insert(fastaFileName);
+ fastaFile2Group[fastaFileName] = comboGroupName;
}
fastaFileNames[itBar->first][itPrimer->first] = fastaFileName;
m->openOutputFile(fastaFileName, temp); temp.close();
+ cout << fastaFileName << endl;
}
- }
- }
- }
-
- bool allBlank = true;
- for (int i = 0; i < barcodeNameVector.size(); i++) {
- if (barcodeNameVector[i] != "") {
- allBlank = false;
- break;
- }
- }
- for (int i = 0; i < primerNameVector.size(); i++) {
- if (primerNameVector[i] != "") {
- allBlank = false;
- break;
- }
- }
-
- if (allBlank) {
- m->mothurOut("[WARNING]: your oligos file does not contain any group names. mothur will not create a groupfile."); m->mothurOutEndLine();
- allFiles = false;
- return false;
- }
-
- return true;
-
+ }
+ }
+ }
+
+ if (allBlank) {
+ m->mothurOut("[WARNING]: your oligos file does not contain any group names. mothur will not create a groupfile."); m->mothurOutEndLine();
+ allFiles = false;
+ return false;
+ }
+
+ return true;
+
}
catch(exception& e) {
m->errorOut(e, "MakeContigsCommand", "getOligos");
exit(1);
}
}
-//********************************************************************/
-string MakeContigsCommand::reverseOligo(string oligo){
- try {
- string reverse = "";
-
- for(int i=oligo.length()-1;i>=0;i--){
-
- if(oligo[i] == 'A') { reverse += 'T'; }
- else if(oligo[i] == 'T'){ reverse += 'A'; }
- else if(oligo[i] == 'U'){ reverse += 'A'; }
-
- else if(oligo[i] == 'G'){ reverse += 'C'; }
- else if(oligo[i] == 'C'){ reverse += 'G'; }
-
- else if(oligo[i] == 'R'){ reverse += 'Y'; }
- else if(oligo[i] == 'Y'){ reverse += 'R'; }
-
- else if(oligo[i] == 'M'){ reverse += 'K'; }
- else if(oligo[i] == 'K'){ reverse += 'M'; }
-
- else if(oligo[i] == 'W'){ reverse += 'W'; }
- else if(oligo[i] == 'S'){ reverse += 'S'; }
-
- else if(oligo[i] == 'B'){ reverse += 'V'; }
- else if(oligo[i] == 'V'){ reverse += 'B'; }
-
- else if(oligo[i] == 'D'){ reverse += 'H'; }
- else if(oligo[i] == 'H'){ reverse += 'D'; }
-
- else { reverse += 'N'; }
- }
-
-
- return reverse;
- }
- catch(exception& e) {
- m->errorOut(e, "MakeContigsCommand", "reverseOligo");
- exit(1);
- }
-}
//**********************************************************************************************************************
vector<int> MakeContigsCommand::convertQual(string qual) {
try {
vector<int> qualScores;
bool negativeScores = false;
- for (int i = 0; i < qual.length(); i++) {
+ for (int i = 0; i < qual.length(); i++) {
int temp = 0;
temp = int(qual[i]);
projects / mothur.git / blobdiff