Find the overlapping SNPs in 4 different data sets.

    :param prefix: the prefix of all the files.
    :param referencePrefixes: the prefix of the reference population files.
    :param referencePopulations: the name of the reference population (same
                                 order as ``referencePrefixes``)
    :param outPrefix: the prefix of the output files.

    :type prefix: str
    :type referencePrefixes: list
    :type referencePopulations: list
    :type outPrefix: str

    It starts by reading the ``bim`` file of the source data set
    (``prefix.bim``). It finds all the markers (excluding the duplicated ones).
    Then it reads all of the reference population ``bim`` files
    (``referencePrefixes.bim``) and find all the markers that were found in the
    source data set.

    It creates three output files:

    * ``outPrefix.ref_snp_to_extract``: the name of the markers that needs to
      be extracted from the three reference panels.
    * ``outPrefix.source_snp_to_extract``: the name of the markers that needs
      to be extracted from the source panel.
    * ``outPrefix.update_names``: a file (readable by Plink) that will help in
      changing the names of the selected markers in the reference panels, so
      that they become comparable with the source panel.

def findOverlappingSNPsWithReference(prefix, referencePrefixes,
                                     referencePopulations, outPrefix):
    """Find the overlapping SNPs in 4 different data sets.

    :param prefix: the prefix of all the files.
    :param referencePrefixes: the prefix of the reference population files.
    :param referencePopulations: the name of the reference population (same
                                 order as ``referencePrefixes``)
    :param outPrefix: the prefix of the output files.

    :type prefix: str
    :type referencePrefixes: list
    :type referencePopulations: list
    :type outPrefix: str

    It starts by reading the ``bim`` file of the source data set
    (``prefix.bim``). It finds all the markers (excluding the duplicated ones).
    Then it reads all of the reference population ``bim`` files
    (``referencePrefixes.bim``) and find all the markers that were found in the
    source data set.

    It creates three output files:

    * ``outPrefix.ref_snp_to_extract``: the name of the markers that needs to
      be extracted from the three reference panels.
    * ``outPrefix.source_snp_to_extract``: the name of the markers that needs
      to be extracted from the source panel.
    * ``outPrefix.update_names``: a file (readable by Plink) that will help in
      changing the names of the selected markers in the reference panels, so
      that they become comparable with the source panel.

    """
    # Reading the main file
    sourceSnpToExtract = {}
    duplicates = set()
    try:
        with open(prefix + ".bim", "r") as inputFile:
            for line in inputFile:
                row = line.rstrip("\r\n").split("\t")
                chromosome = row[0]
                position = row[3]
                snpName = row[1]

                if (chromosome, position) not in sourceSnpToExtract:
                    sourceSnpToExtract[(chromosome, position)] = snpName
                else:
                    # It's a duplicate
                    duplicates.add((chromosome, position))
    except IOError:
        msg = "%s.bim: no such file" % prefix
        raise ProgramError(msg)

    # Removing duplicates from the list
    for snpID in duplicates:
        del sourceSnpToExtract[snpID]

    # Reading each of the reference markers
    refSnpToExtract = defaultdict(dict)
    refSnp = defaultdict(set)
    for refPop, refPrefix in izip(referencePopulations, referencePrefixes):
        try:
            with open(refPrefix + ".bim", "r") as inputFile:
                for line in inputFile:
                    row = line.rstrip("\r\n").split("\t")
                    chromosome = row[0]
                    position = row[3]
                    snpName = row[1]

                    key = (chromosome, position)
                    if key in sourceSnpToExtract:
                        # We want this SNP
                        refSnpToExtract[refPop][key] = snpName
                        refSnp[refPop].add(key)
                logger.info("  - {:,d} overlaps with {}".format(
                    len(refSnp[refPop]),
                    refPrefix,
                ))

        except IOError:
            msg = "%(refPrefix)s.bim: no such file" % locals()
            raise ProgramError(msg)

    # Creating the intersect of the reference SNP
    refSnpIntersect = refSnp[referencePopulations[0]]
    for refPop in referencePopulations[1:]:
        refSnpIntersect &= refSnp[refPop]
    logger.info("  - {:,d} in common between reference "
                "panels".format(len(refSnpIntersect)))

    # Printing the names of the SNPs to extract in the reference
    refOutputFiles = {}
    try:
        for refPop in referencePopulations:
            refOutputFiles[refPop] = open(
                outPrefix + ".{}_snp_to_extract".format(refPop),
                mode="w",
            )
    except IOError:
        msg = "{}.POP_snp_to_extract: can't write file".format(outPrefix)
        raise ProgramError(msg)

    sourceOutputFile = None
    try:
        sourceOutputFile = open(outPrefix + ".source_snp_to_extract", "w")
    except IOError:
        msg = "{}.source_snp_to_extract: can't write file".format(outPrefix)
        raise ProgramError(msg)

    changeNameOutputFiles = {}
    try:
        for refPop in referencePopulations:
            changeNameOutputFiles[refPop] = open(
                outPrefix + ".{}_update_names".format(refPop),
                mode="w",
            )
    except IOError:
        msg = "%(outPrefix)s.updateNames: can't write file" % locals()
        raise ProgramError(msg)

    # Writing the file containing the SNPs to extract in the source
    for snpID in refSnpIntersect:
        print >>sourceOutputFile, sourceSnpToExtract[snpID]
        for refPop in referencePopulations:
            print >>refOutputFiles[refPop], refSnpToExtract[refPop][snpID]
            print >>changeNameOutputFiles[refPop], "\t".join([
                refSnpToExtract[refPop][snpID],
                sourceSnpToExtract[snpID],
            ])

    # Closing the output file
    sourceOutputFile.close()
    for refOutputFile in refOutputFiles.itervalues():
        refOutputFile.close()
    for changeNameOutputFile in changeNameOutputFiles.itervalues():
        changeNameOutputFile.close()