User request: constrain analysis to one parent sequence

[ArtPoon]

Rather than screening the entire input alignment (default behaviour), or a second set of sequences (-r option) for potential parent sequences, this new option would enable the user to specify a specific sequence as one of the parents, and any of the other sequences in the input alignment as the second parent. This would reduce the complexity from quadratic to linear.

[ArtPoon]

@WilliamZekaiWang - I think we just have to change how combinations is being generated, based on an additional argument to the TripletGenerator class constructor:

OpenRDP/openrdp/common.py

Lines 209 to 217 in 0a899b0

	if isinstance(self.ref_align, np.ndarray):
	self.combinations = product(
	list(combinations(range(self.alignment.shape[0]), 1)),
	list(combinations(range(self.ref_align.shape[0]), 2))
	) # Example: ((1,), (1, 2))
	else:
	self.combinations = combinations(
	range(self.alignment.shape[0]), # number of "rows" (sequences)
	3) # all combinations of three elements

[WilliamZekaiWang]

I think I've got something? I made a new argument for argsparse -l that lets users specify a sequence name in the reference sequence to always be included:

>openrdp test.fasta -r ref.fasta -l Sequence_7

Method          Start   End     Recombinant     Parent1 Parent2 Pvalue
------------------------------------------------------------------------
Geneconv        3       7       Sequence_2      Sequence_4      -       3.92E-03
Geneconv        12      52      Sequence_3      Sequence_4      -       4.86E-01
Bootscan        100     120     Sequence_5      Sequence_6      Sequence_7      5.74E-06
Bootscan        100     120     Sequence_5      Sequence_7      Sequence_8      2.15E-06
Bootscan        100     120     Sequence_7      Sequence_5      Sequence_6      2.15E-06
Bootscan        100     120     Sequence_8      Sequence_4      Sequence_7      4.14E-29
3Seq            7       52      Sequence_4      Sequence_1      Sequence_2      4.12E-02

3Seq does not use our TripletGenerator and so, I'm not too sure on how to incorporate these changes to it

[GopiGugan]

Lets see if there is an option in 3Seq to specify a sequence - https://mol.ax/content/media/2018/02/3seq_manual.20180209.pdf

[ArtPoon]

It might be worth isolating the core 3Seq algorithm and then hooking that C code into Python so that we can feed inputs directly.

[ArtPoon]

3Seq includes the following command line options:

      ./3seq.macOS -triplet  <seq_file>  <P_name>    <Q_name>      <C_name>    [options]

      ./3seq.macOS -single   <seq_file>                 [options]
      ./3seq.macOS -single   <parent_file> <child_file> [options]

      ./3seq.macOS -full     <seq_file>                 [-ptable ptable_file]  [options]
      ./3seq.macOS -full     <parent_file> <child_file> [-ptable ptable_file]  [options]

We could pass triplets from TripletGenerator individually in -triplet mode, but this might be inefficient

[WilliamZekaiWang]

I looked through the options:

the triplet mode doesn't seem to report a p-value:

Child accession :       C

Parent P
   Accession :          A
   Distance to child :  18nt  (18nt ignoring gaps)

Parent Q
   Accession :          B
   Distance to child :  17nt  (17nt ignoring gaps)

Informative sites of type P :  3
Informative sites of type Q :  4
Maximum Descent :              2

Breakpoint Ranges:
  18-18 & 20-22
  13-14 & 20-22

# Breakpoint Pairs: 9
# Breakpoint Ranges: 2

Minimum length of recombinant segment is 3nt 
(excluding gaps and ambiguous nucleotides)

single is for single-breakpoint recombinants only and goes through the whole file

[ArtPoon]

It might be worth isolating the core 3Seq algorithm and then hooking that C code into Python so that we can feed inputs directly.

Note that if we went this route for 3Seq, we would also want to do this for GENECONV (fortunately the source code is available)

[ArtPoon]

two possible approaches that come to mind for 3Seq:

1. locate a function (if one exists) to which we can pass three strings as inputs from Python
2. run multiple instances of 3Seq in parallel, where each instance loads the same input file but is also passed an integer index that determines which subset of triplets it will process

[GopiGugan]

the triplet mode doesn't seem to report a p-value:

• I took a look at the 3Seq code. The latest version of 3Seq does not allow the user to use the -triplet mode. This mode was removed in relatively recent commits. However, looking back at a previous version of 3Seq, we can use the TripletRun::perform() function:

https://gitlab.com/lamhm/3seq/-/blob/0708fa30b7903475132e67b1462a8beb147b4798/source/run-mode/TripletRun.cpp#L92-120

• We can use the calculatePVal() function in Triplet.cpp to calculate the p value:

https://gitlab.com/lamhm/3seq/-/blob/0708fa30b7903475132e67b1462a8beb147b4798/source/data/bio/Triplet.cpp#L346-362

• Triplet constructor takes in 3 Genome objects:

https://gitlab.com/lamhm/3seq/-/blob/0708fa30b7903475132e67b1462a8beb147b4798/source/data/bio/Triplet.h#L33-36