Skip to content

Commit 41a7406

Browse files
h-2h-2
authored
Merge pull request #207 from sarahet/bs_mode
Add parameters for bisulfite profiles
2 parents df1c74d + 817446a commit 41a7406

File tree

2 files changed

+164
-126
lines changed

2 files changed

+164
-126
lines changed

‎src/search_options.hpp‎

Lines changed: 85 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -240,6 +240,11 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
240240
.validator = sharg::arithmetic_range_validator{0, 100}
241241
});
242242

243+
if (options.domain == domain_t::bisulfite)
244+
{
245+
options.minBitScore = 68;
246+
}
247+
243248
parser.add_option(options.minBitScore,
244249
sharg::config{
245250
.short_id = '\0',
@@ -282,18 +287,9 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
282287

283288
std::string samBamSeqDescr;
284289

285-
//TODONE: these are TODOs for you
286-
if (options.domain != domain_t::protein) //TODO add separate for BS
290+
if (options.domain != domain_t::protein)
287291
{
288-
samBamSeqDescr = "Write matching DNA subsequence into SAM/BAM file.";
289-
options.searchOpts0.seedLength = 14;
290-
options.searchOpts0.seedOffset = 9;
291-
options.searchOpts0.maxSeedDist = 0;
292-
options.searchOpts.seedLength = 14;
293-
options.searchOpts.seedOffset = 7;
294-
options.searchOpts.maxSeedDist = 1;
295-
options.preScoringThresh = 1.4;
296-
292+
samBamSeqDescr = "Write matching DNA subsequence into SAM/BAM file.";
297293
options.gapOpen = -5;
298294
options.gapExtend = -2;
299295
}
@@ -306,18 +302,40 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
306302
"there is no DNA sequence so a \"*\" is written to the SEQ column. The matching protein sequence "
307303
"can be "
308304
"written as an optional tag, see --sam-bam-tags.";
309-
310-
options.searchOpts0.seedLength = 10;
311-
options.searchOpts0.seedOffset = 5;
312-
options.searchOpts0.maxSeedDist = 0;
313-
options.searchOpts.seedLength = 11;
314-
options.searchOpts.seedOffset = 3;
315-
options.searchOpts.maxSeedDist = 1;
316-
317305
options.gapOpen = -11;
318306
options.gapExtend = -1;
319307
}
320308

309+
switch (options.domain)
310+
{
311+
case domain_t::protein:
312+
options.searchOpts0.seedLength = 10;
313+
options.searchOpts0.seedOffset = 5;
314+
options.searchOpts0.maxSeedDist = 0;
315+
options.searchOpts.seedLength = 11;
316+
options.searchOpts.seedOffset = 3;
317+
options.searchOpts.maxSeedDist = 1;
318+
break;
319+
case domain_t::nucleotide:
320+
options.searchOpts0.seedLength = 14;
321+
options.searchOpts0.seedOffset = 9;
322+
options.searchOpts0.maxSeedDist = 0;
323+
options.searchOpts.seedLength = 14;
324+
options.searchOpts.seedOffset = 7;
325+
options.searchOpts.maxSeedDist = 1;
326+
options.preScoringThresh = 1.4;
327+
break;
328+
case domain_t::bisulfite:
329+
options.searchOpts0.seedLength = 17;
330+
options.searchOpts0.seedOffset = 10;
331+
options.searchOpts0.maxSeedDist = 0;
332+
options.searchOpts.seedLength = 17;
333+
options.searchOpts.seedOffset = 10;
334+
options.searchOpts.maxSeedDist = 1;
335+
options.preScoringThresh = 1.5;
336+
break;
337+
}
338+
321339
std::string samBamSeqDescrTmp = "uniq";
322340
parser.add_option(samBamSeqDescrTmp,
323341
sharg::config{
@@ -545,22 +563,29 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
545563
"arguments!",
546564
true);
547565

548-
if (options.domain != domain_t::protein) //TODO add separate for BS
566+
switch (options.domain)
549567
{
550-
parser.add_line("\"fast\"", false);
551-
parser.add_line("--seed-length 14 --seed-offset 9 --seed-delta 0", true);
552-
parser.add_line("\"sensitive\"", false);
553-
parser.add_line("--seed-length0 14 --seed-offset0 3 --seed-length 14 --seed-offset 3", true);
554-
}
555-
else
556-
{
557-
parser.add_line("\"fast\"", false);
558-
parser.add_line("--seed-length0 12 --seed-offset0 8 --seed-length 10 --seed-offset 5 --seed-delta 0", true);
559-
parser.add_line("\"sensitive\"", false);
560-
parser.add_line(
561-
"--seed-length0 9 --seed-offset0 4 --seed-length 8 --seed-offset 3 --pre-scoring 3 "
562-
"--pre-scoring-threshold 1.9",
563-
true);
568+
case domain_t::protein:
569+
parser.add_line("\"fast\"", false);
570+
parser.add_line("--seed-length0 12 --seed-offset0 8 --seed-length 10 --seed-offset 5 --seed-delta 0", true);
571+
parser.add_line("\"sensitive\"", false);
572+
parser.add_line(
573+
"--seed-length0 9 --seed-offset0 4 --seed-length 8 --seed-offset 3 --pre-scoring 3 "
574+
"--pre-scoring-threshold 1.9",
575+
true);
576+
break;
577+
case domain_t::nucleotide:
578+
parser.add_line("\"fast\"", false);
579+
parser.add_line("--seed-length 14 --seed-offset 9 --seed-delta 0", true);
580+
parser.add_line("\"sensitive\"", false);
581+
parser.add_line("--seed-length0 14 --seed-offset0 3 --seed-length 14 --seed-offset 3", true);
582+
break;
583+
case domain_t::bisulfite:
584+
parser.add_line("\"fast\"", false);
585+
parser.add_line("--seed-length 17 --seed-offset 10 --seed-delta 0", true);
586+
parser.add_line("\"sensitive\"", false);
587+
parser.add_line("--seed-length0 16 --seed-offset0 8 --seed-length 15 --seed-offset 10", true);
588+
break;
564589
}
565590
parser.add_line("For further information see the wiki: <https://github.com/seqan/lambda/wiki>", false);
566591

@@ -584,11 +609,15 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
584609

585610
if (options.profile == "fast")
586611
{
587-
if (options.domain != domain_t::protein) //TODO add separate for BS
612+
if (options.domain != domain_t::protein)
588613
{
589614
options.iterativeSearch = false;
590-
options.searchOpts.seedOffset = 9;
591615
options.searchOpts.maxSeedDist = 0;
616+
617+
if (options.domain == domain_t::nucleotide)
618+
{
619+
options.searchOpts.seedOffset = 9;
620+
}
592621
}
593622
else
594623
{
@@ -601,20 +630,27 @@ void parseCommandLine(LambdaOptions & options, int argc, char const ** argv)
601630
}
602631
else if (options.profile == "sensitive")
603632
{
604-
if (options.domain != domain_t::protein) //TODO add separate for BS
633+
switch (options.domain)
605634
{
606-
options.searchOpts0.seedOffset = 3;
607-
options.searchOpts.seedOffset = 3;
608-
}
609-
else
610-
{
611-
options.searchOpts0.seedLength = 9;
612-
options.searchOpts0.seedOffset = 4;
613-
options.searchOpts.seedLength = 8;
614-
options.searchOpts.seedOffset = 3;
615-
616-
options.preScoring = 3;
617-
options.preScoringThresh = 1.9;
635+
case domain_t::protein:
636+
options.searchOpts0.seedLength = 9;
637+
options.searchOpts0.seedOffset = 4;
638+
options.searchOpts.seedLength = 8;
639+
options.searchOpts.seedOffset = 3;
640+
641+
options.preScoring = 3;
642+
options.preScoringThresh = 1.9;
643+
break;
644+
case domain_t::nucleotide:
645+
options.searchOpts0.seedOffset = 3;
646+
options.searchOpts.seedOffset = 3;
647+
break;
648+
case domain_t::bisulfite:
649+
options.searchOpts0.seedLength = 16;
650+
options.searchOpts0.seedOffset = 8;
651+
options.searchOpts.seedLength = 15;
652+
options.searchOpts.seedOffset = 10;
653+
break;
618654
}
619655
}
620656

0 commit comments

Comments
 (0)