Extract related individuals according IBS2* ratio.

    :param fileName: the name of the input file.
    :param outPrefix: the prefix of the output files.
    :param ibs2_ratio_threshold: the ibs2 ratio threshold (tells if sample pair
                                 is related or not).

    :type fileName: str
    :type outPrefix: str
    :type ibs2_ratio_threshold: float

    :returns: a :py:class:`numpy.recarray` data set containing (for each
              related sample pair) the ``ibs2 ratio``, ``Z1``, ``Z2`` and the
              type of relatedness.

    Reads a ``genome`` file (provided by :py:func:`runGenome`) and extract
    related sample pairs according to ``IBS2 ratio``.

    A ``genome`` file contains at least the following information for each
    sample pair:

    * **FID1:** the family ID of the first sample in the pair.
    * **IID1:** the individual ID of the first sample in the pair.
    * **FID2:** the family ID of the second sample in the pair.
    * **IID2:** the individual ID of the second sample in the pair.
    * **Z0:** the probability that :math:`IBD = 0`.
    * **Z1:** the probability that :math:`IBD = 1`.
    * **Z2:** the probability that :math:`IBD = 2`.
    * **HOMHOM:** the number of :math:`IBS = 0` SNP pairs used in ``PPC`` test.
    * **HETHET:** the number of :math:`IBS = 2` het/het SNP pairs in ``PPC``
      test.

    The ``IBS2 ratio`` is computed using the following formula:

    .. math::
        \\textrm{IBS2 ratio} = \\frac{\\textrm{HETHET}}
                                     {\\textrm{HOMHOM} + \\textrm{HETHET}}

    If the ``IBS2 ratio`` is higher than the threshold, the samples in the pair
    are related. The following values help in finding the relatedness of the
    sample pair.

    +---------------------------------------+-----------------------+------+
    |           Values                      |        Relation       | Code |
    +=======================================+=======================+======+
    | :math:`0.17 \\leq z_0 \\leq 0.33` and   | Full-sibs             | 1    |
    | :math:`0.40 \\leq z_1 \\leq 0.60`       |                       |      |
    +---------------------------------------+-----------------------+------+
    | :math:`0.40 \\leq z_0 \\leq 0.60` and   | Half-sibs or          | 2    |
    | :math:`0.40 \\leq z_1 \\leq 0.60`       | Grand-parent-Child or |      |
    |                                       | Uncle-Nephew          |      |
    +---------------------------------------+-----------------------+------+
    | :math:`z_0 \\leq 0.05` and             | Parent-Child          | 3    |
    | :math:`z_1 \\geq 0.95` and             |                       |      |
    | :math:`z_2 \\leq 0.05`                 |                       |      |
    +---------------------------------------+-----------------------+------+
    | :math:`z_0 \\leq 0.05` and             | Twins or Duplicated   | 4    |
    | :math:`z_1 \\leq 0.05` and             | samples               |      |
    | :math:`z_2 \\geq 0.95`                 |                       |      |
    +---------------------------------------+-----------------------+------+

def extractRelatedIndividuals(fileName, outPrefix, ibs2_ratio_threshold):
    """Extract related individuals according IBS2* ratio.

    :param fileName: the name of the input file.
    :param outPrefix: the prefix of the output files.
    :param ibs2_ratio_threshold: the ibs2 ratio threshold (tells if sample pair
                                 is related or not).

    :type fileName: str
    :type outPrefix: str
    :type ibs2_ratio_threshold: float

    :returns: a :py:class:`numpy.recarray` data set containing (for each
              related sample pair) the ``ibs2 ratio``, ``Z1``, ``Z2`` and the
              type of relatedness.

    Reads a ``genome`` file (provided by :py:func:`runGenome`) and extract
    related sample pairs according to ``IBS2 ratio``.

    A ``genome`` file contains at least the following information for each
    sample pair:

    * **FID1:** the family ID of the first sample in the pair.
    * **IID1:** the individual ID of the first sample in the pair.
    * **FID2:** the family ID of the second sample in the pair.
    * **IID2:** the individual ID of the second sample in the pair.
    * **Z0:** the probability that :math:`IBD = 0`.
    * **Z1:** the probability that :math:`IBD = 1`.
    * **Z2:** the probability that :math:`IBD = 2`.
    * **HOMHOM:** the number of :math:`IBS = 0` SNP pairs used in ``PPC`` test.
    * **HETHET:** the number of :math:`IBS = 2` het/het SNP pairs in ``PPC``
      test.

    The ``IBS2 ratio`` is computed using the following formula:

    .. math::
        \\textrm{IBS2 ratio} = \\frac{\\textrm{HETHET}}
                                     {\\textrm{HOMHOM} + \\textrm{HETHET}}

    If the ``IBS2 ratio`` is higher than the threshold, the samples in the pair
    are related. The following values help in finding the relatedness of the
    sample pair.

    +---------------------------------------+-----------------------+------+
    |           Values                      |        Relation       | Code |
    +=======================================+=======================+======+
    | :math:`0.17 \\leq z_0 \\leq 0.33` and   | Full-sibs             | 1    |
    | :math:`0.40 \\leq z_1 \\leq 0.60`       |                       |      |
    +---------------------------------------+-----------------------+------+
    | :math:`0.40 \\leq z_0 \\leq 0.60` and   | Half-sibs or          | 2    |
    | :math:`0.40 \\leq z_1 \\leq 0.60`       | Grand-parent-Child or |      |
    |                                       | Uncle-Nephew          |      |
    +---------------------------------------+-----------------------+------+
    | :math:`z_0 \\leq 0.05` and             | Parent-Child          | 3    |
    | :math:`z_1 \\geq 0.95` and             |                       |      |
    | :math:`z_2 \\leq 0.05`                 |                       |      |
    +---------------------------------------+-----------------------+------+
    | :math:`z_0 \\leq 0.05` and             | Twins or Duplicated   | 4    |
    | :math:`z_1 \\leq 0.05` and             | samples               |      |
    | :math:`z_2 \\geq 0.95`                 |                       |      |
    +---------------------------------------+-----------------------+------+

    """
    # The output file
    outputFile = None
    try:
        outputFile = open(outPrefix + ".related_individuals", "w")
    except IOError:
        msg = "%(outPrefix)s.related_individuals: can't write file" % locals()
        raise ProgramError(msg)

    # The input file
    inputFile = None
    try:
        if fileName.endswith(".gz"):
            inputFile = gzip.open(fileName, 'rb')
        else:
            inputFile = open(fileName, 'r')
    except IOError:
        msg = "{}: no such file".format(fileName)
        raise ProgramError(msg)

    headerIndex = None
    data = []
    for i, line in enumerate(inputFile):
        row = createRowFromPlinkSpacedOutput(line)

        if i == 0:
            # This is the header
            headerIndex = dict([(colName, j) for j, colName in enumerate(row)])

            # Checking columns
            for columnName in ["FID1", "IID1", "FID2", "IID2", "Z0", "Z1",
                               "Z2", "HOMHOM", "HETHET"]:
                if columnName not in headerIndex:
                    msg = "{}: no culumn named {}".format(fileName, columnName)
                    raise ProgramError(msg)

            # Writing header
            print >>outputFile, "\t".join(row + ["IBS2_ratio", "status",
                                                 "code"])

        else:
            # This is data
            homhom = row[headerIndex["HOMHOM"]]
            hethet = row[headerIndex["HETHET"]]
            try:
                homhom = float(homhom)
                hethet = float(hethet)
            except ValueError:
                msg = "{}: invalid HOMHOM or HETHET".format(fileName)
                raise ProgramError(float)

            # Computing IBS2* ratio
            ibs2_ratio = hethet / (homhom + hethet)

            if ibs2_ratio > ibs2_ratio_threshold:
                # Those pairs might be related
                # Finding the status
                status = "unknown"
                code = "5"
                z0 = row[headerIndex["Z0"]]
                z1 = row[headerIndex["Z1"]]
                z2 = row[headerIndex["Z2"]]
                try:
                    z0 = float(z0)
                    z1 = float(z1)
                    z2 = float(z2)
                except ValueError:
                    msg = "{}: invalid value for Z0, Z1 or Z2".format(fileName)
                    raise ProgramError(msg)

                if (z0 >= 0.17 and z0 <= 0.33) and (z1 >= 0.40 and z1 <= 0.60):
                    # Full sibs
                    status = "full-sibs"
                    code = "1"

                elif (z0 >= 0.4 and z0 <= 0.6) and (z1 >= 0.4 and z1 <= 0.6):
                    # half sibs, grand-parent child, uncle nephew
                    status = ";".join(["half-sibs", "grand-parent-child",
                                       "uncle-nephew"])
                    code = "2"

                elif (z0 <= 0.05) and (z1 >= 0.95) and (z2 <= 0.05):
                    # parent child
                    status = "parent-child"
                    code = "3"

                elif (z0 <= 0.05) and (z1 <= 0.05) and (z2 >= 0.95):
                    # twin
                    status = "twins"
                    code = "4"

                # Printing to file
                print >>outputFile, "\t".join(row + [str(ibs2_ratio), status,
                                                     code])
                data.append((ibs2_ratio, z1, z2, code))

    # Closing the output and input files
    inputFile.close()
    outputFile.close()

    # Merging the related individuals
    merge_related_samples(outPrefix + ".related_individuals", outPrefix, False)

    # If there are no related samples, we return nothing
    if len(data) == 0:
        return None

    # Creating the numpy array if there are related samples
    data = np.array(data, dtype=[
        ("IBS2_RATIO", float),
        ("Z1", float),
        ("Z2", float),
        ("CODE", "S{}".format(max([len(i[3]) for i in data]))),
    ])

    return data