v1.2.2 added --save_control CLI option (#42)

Author: Marc Maxmeister

methylprep/cli.py

@@ -176,6 +176,14 @@ def cli_process(cmd_args):
         help="Change the processed beta or m_value data_type output from float64 to float16 or float32, to save disk space.",
     )
 
+    parser.add_argument(
+        '-c', '--save_control',
+        required=False,
+        action='store_true',
+        default=False,
+        help='If specified, saves an additional "control_probes.pkl" file that contains Control and SNP-I probe data in the data_dir.'
+    )
+
     args = parser.parse_args(cmd_args)
 
     array_type = args.array_type
@@ -200,6 +208,7 @@ def cli_process(cmd_args):
         export=args.no_export, # flag flips here
         meta_data_frame=args.no_meta_export, # flag flips here
         bit=args.bit,
+        save_control=args.save_control,
         )
 
 

methylprep/files/idat.py

@@ -39,6 +39,9 @@ class IdatHeaderLocation(IntEnum):
 class IdatSectionCode(IntEnum):
     """Unique IntEnum defining constant values for byte offsets of IDAT headers.
     These values come from the field codes of the Bioconductor illuminaio package.
+
+    MM: refer to https://bioconductor.org/packages/release/bioc/vignettes/illuminaio/inst/doc/EncryptedFormat.pdf
+    and https://bioconductor.org/packages/release/bioc/vignettes/illuminaio/inst/doc/illuminaio.pdf
     """
     ILLUMINA_ID = 102
     STD_DEV = 103
@@ -239,7 +242,7 @@ def seek_to_section(section_code):
             data=probe_records,
             orient='index',
             columns=['mean_value'],
-            dtype='float64',
+            dtype='float32',
         )
 
         data_frame.index.name = 'illumina_id'

methylprep/files/manifests.py

@@ -50,6 +50,7 @@
     'Control_Type',
     'Color',
     'Extended_Type',
+    # control probes don't have 'IlmnID' values set -- these probes are not locii specific
 )
 
 
@@ -80,6 +81,7 @@ def __init__(self, array_type, filepath_or_buffer=None, on_lambda=False):
         with get_file_object(filepath_or_buffer) as manifest_file:
             self.__data_frame = self.read_probes(manifest_file)
             self.__control_data_frame = self.read_control_probes(manifest_file)
+            self.__snp_data_frame = self.read_snp_probes(manifest_file)
 
     @property
     def columns(self):
@@ -93,6 +95,10 @@ def data_frame(self):
     def control_data_frame(self):
         return self.__control_data_frame
 
+    @property
+    def snp_data_frame(self):
+        return self.__snp_data_frame
+
     @staticmethod
     def download_default(array_type, on_lambda=False):
         """Downloads the appropriate manifest file if one does not already exist.
@@ -150,6 +156,13 @@ def read_probes(self, manifest_file):
         )
 
         def get_probe_type(name, infinium_type):
+            """what:
+        returns one of (I, II, SnpI, SnpII, Control)
+
+        how:
+            .from_manifest_values() returns probe type using either
+            the Infinium_Design_Type (I or II) or the name (starts with 'rs' == SnpI)
+            and 'Control' is none of the above."""
             probe_type = ProbeType.from_manifest_values(name, infinium_type)
             return probe_type.value
 
@@ -163,23 +176,37 @@ def get_probe_type(name, infinium_type):
         return data_frame
 
     def read_control_probes(self, manifest_file):
-        LOGGER.info(f'Reading control probes: {Path(manifest_file.name).stem}')
+        """ Unlike other probes, control probes have no IlmnID because they're not locus-specific. """
+        #LOGGER.info(f'Reading control probes: {Path(manifest_file.name).stem}')
 
         self.seek_to_start(manifest_file)
 
-        num_headers = 4
+        #num_headers = 4 -- removed on Jan 21 2020.
+        # Steve Byerly added this =4, but the manifests seem to start controls at the point specified, with no extra columns based on .num_probes stored.
+        num_headers = 0
 
         return pd.read_csv(
             manifest_file,
             comment='[',
             header=None,
-            index_col=0,
+            index_col=0, # illumina_id, not IlmnID here
             names=CONTROL_COLUMNS,
             nrows=self.array_type.num_controls,
             skiprows=self.array_type.num_probes + num_headers,
             usecols=range(len(CONTROL_COLUMNS)),
         )
 
+    def read_snp_probes(self, manifest_file):
+        """ Unlike cpg and control probes, these rs probes are NOT sequential in all arrays. """
+        #LOGGER.info(f'Reading snp probes: {Path(manifest_file.name).stem} --> {snp_df.shape[0]} found')
+        self.seek_to_start(manifest_file)
+        # since these are not sequential, loading everything and filtering by IlmnID.
+        snp_df = pd.read_csv(
+            manifest_file,
+            low_memory=False)
+        snp_df = snp_df[snp_df['IlmnID'].str.match('rs', na=False)]
+        return snp_df
+
     def map_to_genome(self, data_frame):
         genome_df = self.get_genome_data()
         merged_df = inner_join_data(data_frame, genome_df)
@@ -221,6 +248,9 @@ def get_loci_names(self):
         return self.data_frame['Name'].unique()
 
     def get_probe_details(self, probe_type, channel=None):
+        """given a probe type (I, II, SnpI, SnpII, Control) and a channel (Channel.RED | Channel.GREEN),
+        This will return info needed to map probes to their names (e.g. cg0031313 or rs00542420),
+        which are NOT in the idat files."""
         if not isinstance(probe_type, ProbeType):
             raise Exception('probe_type is not a valid ProbeType')
 

methylprep/models/__init__.py

@@ -6,6 +6,8 @@
     FG_RED_PROBE_SUBSETS,
     METHYLATED_PROBE_SUBSETS,
     UNMETHYLATED_PROBE_SUBSETS,
+    METHYLATED_SNP_PROBES,
+    UNMETHYLATED_SNP_PROBES,
     Channel,
     Probe,
     ProbeAddress,

methylprep/models/arrays.py

@@ -15,7 +15,7 @@ def __str__(self):
 
     @classmethod
     def from_probe_count(cls, probe_count):
-        """Determines array type using number of probes. Returns array string."""
+        """Determines array type using number of probes counted in raw idat file. Returns array string."""
         if probe_count == 1055583:
             return cls.ILLUMINA_EPIC_PLUS
 
@@ -36,20 +36,28 @@ def num_probes(self):
             ArrayType.ILLUMINA_27K: 27578,
             ArrayType.ILLUMINA_450K: 485578,
             ArrayType.ILLUMINA_EPIC: 865919,
-            ArrayType.ILLUMINA_EPIC_PLUS: 868699,
+            ArrayType.ILLUMINA_EPIC_PLUS: 868698, # was 868699 until Jan 21, 2020. corrected.
         }
-
         return probe_counts.get(self)
 
     @property
     def num_controls(self):
         probe_counts = {
-            ArrayType.ILLUMINA_27K: 144,
+            ArrayType.ILLUMINA_27K: 144, # the manifest does not contain control probe data (illumina's site included)
             ArrayType.ILLUMINA_450K: 850,
             ArrayType.ILLUMINA_EPIC: 635,
             ArrayType.ILLUMINA_EPIC_PLUS: 635,
         }
+        return probe_counts.get(self)
 
+    @property
+    def num_snps(self):
+        probe_counts = {
+            ArrayType.ILLUMINA_27K: 0,
+            ArrayType.ILLUMINA_450K: 65,
+            ArrayType.ILLUMINA_EPIC: 59,
+            ArrayType.ILLUMINA_EPIC_PLUS: 120,
+        }
         return probe_counts.get(self)
 
 ''' # doesn't appear to be used anywhere. -- and pipeline Array model conflicts with it.

methylprep/models/probes.py

@@ -4,6 +4,8 @@
 
 @unique
 class Channel(Enum):
+    """ idat probes measure either a red or green fluorescence.
+    This specifies which to return within idat.py: red_idat or green_idat."""
     RED = 'Red'
     GREEN = 'Grn'
 
@@ -21,6 +23,13 @@ def is_red(self):
 
 @unique
 class ProbeAddress(Enum):
+    """AddressA_ID and AddressB_ID are columns in the manifest csv that contain internal Illumina probe identifiers.
+
+    Type II probes use AddressA_ID; Type I uses both, because there are two probes, two colors.
+
+    probe intensities in .idat files are keyed to one of these ids, but processed data is always keyed to the
+    IlmnID probe "names" -- so this is used in converting between IDs. It is used to define probe sets below
+    in this probes.py"""
     A = 'A'
     B = 'B'
 
@@ -33,6 +42,11 @@ def header_name(self):
 
 @unique
 class ProbeType(Enum):
+    """ probes can either be type I or type II for CpG or Snp sequences. Control probes are used for background
+    correction in different fluorescence ranges and staining efficiency.
+    Type I probes record EITHER a red or a green value.
+    Type II probes record both values together.
+    NOOB uses the red fluorescence on a green probe and vice versa to calculate background fluorescence."""
     ONE = 'I'
     TWO = 'II'
     SNP_ONE = 'SnpI'
@@ -44,6 +58,9 @@ def __str__(self):
 
     @staticmethod
     def from_manifest_values(name, infinium_type):
+        """ this function determines which of four kinds of probe goes with this name, using either
+        the Infinium_Design_Type (I or II) or the name (starts with 'rs')
+        and decides 'Control' is non of the above."""
         is_snp = name.startswith('rs')
 
         if infinium_type == 'I':
@@ -56,6 +73,7 @@ def from_manifest_values(name, infinium_type):
 
 
 class Probe():
+    """ this doesn't appear to be instantiated anywhere in methylprep """
     __slots__ = [
         'address',
         'illumina_id',
@@ -69,6 +87,9 @@ def __init__(self, address, illumina_id, probe_type):
 
 
 class ProbeSubset():
+    """ used below in probes.py to define sub-sets of probes:
+    foreground-(red|green|all), or (un)methylated probes
+    """
     __slots__ = [
         'data_channel',
         'probe_address',
@@ -104,6 +125,81 @@ def get_probe_details(self, manifest):
         )
 
 
+'''
+# notebook equivalent of SNP_PROBES below.
+#manifest = data_containers[0].manifest.data_frame
+#snpProbesI = manifest[manifest['probe_type']=='SnpI']
+#snpProbesI_Grn = snpProbesI[snpProbesI['Color_Channel']=='Grn']
+#snpProbesI_Red = snpProbesI[snpProbesI['Color_Channel']=='Red']
+#snpProbesII = manifest[manifest['probe_type']=='SnpII']
+
+SNP_PROBES = (
+    # SNP_II_PROBES
+    ProbeSubset(
+        data_channel=Channel.GREEN,
+        probe_address=ProbeAddress.A,
+        probe_channel=None,
+        probe_type=ProbeType.SNP_TWO,
+    ),
+    # SNP_I_GREEN_PROBES
+    ProbeSubset(
+        data_channel=Channel.GREEN,
+        probe_address=ProbeAddress.A,
+        probe_channel=Channel.GREEN,
+        probe_type=ProbeType.SNP_ONE,
+    ),
+    # SNP_I_RED_PROBES
+    ProbeSubset(
+        data_channel=Channel.RED,
+        probe_address=ProbeAddress.B,
+        probe_channel=Channel.RED,
+        probe_type=ProbeType.SNP_ONE,
+    )
+)
+'''
+
+METHYLATED_SNP_PROBES = (
+    ProbeSubset(
+        data_channel=Channel.GREEN,
+        probe_address=ProbeAddress.A,
+        probe_channel=None,
+        probe_type=ProbeType.SNP_TWO,
+    ),
+    ProbeSubset(
+        data_channel=Channel.GREEN,
+        probe_address=ProbeAddress.B,
+        probe_channel=Channel.GREEN,
+        probe_type=ProbeType.SNP_ONE,
+    ),
+    ProbeSubset(
+        data_channel=Channel.RED,
+        probe_address=ProbeAddress.B,
+        probe_channel=Channel.RED,
+        probe_type=ProbeType.SNP_ONE,
+    ),
+)
+
+UNMETHYLATED_SNP_PROBES = (
+    ProbeSubset(
+        data_channel=Channel.RED,
+        probe_address=ProbeAddress.A,
+        probe_channel=None,
+        probe_type=ProbeType.SNP_TWO,
+    ),
+    ProbeSubset(
+        data_channel=Channel.GREEN,
+        probe_address=ProbeAddress.A,
+        probe_channel=Channel.GREEN,
+        probe_type=ProbeType.SNP_ONE,
+    ),
+    ProbeSubset(
+        data_channel=Channel.RED,
+        probe_address=ProbeAddress.A,
+        probe_channel=Channel.RED,
+        probe_type=ProbeType.SNP_ONE,
+    ),
+)
+
 FG_GREEN_PROBE_SUBSETS = (
     ProbeSubset(
         data_channel=Channel.GREEN,

methylprep/processing/meth_dataset.py

@@ -2,7 +2,12 @@
 import logging
 import pandas as pd
 # App
-from ..models import METHYLATED_PROBE_SUBSETS, UNMETHYLATED_PROBE_SUBSETS
+from ..models import (
+    METHYLATED_PROBE_SUBSETS,
+    UNMETHYLATED_PROBE_SUBSETS,
+    METHYLATED_SNP_PROBES,
+    UNMETHYLATED_SNP_PROBES,
+)
 
 
 __all__ = ['MethylationDataset']
@@ -25,7 +30,7 @@ class MethylationDataset():
     __preprocessed = False
 
     def __init__(self, raw_dataset, manifest, probe_subsets):
-        LOGGER.info('Preprocessing methylation dataset: %s', raw_dataset.sample)
+        #LOGGER.info('Preprocessing methylation dataset: %s', raw_dataset.sample)
 
         self.probe_subsets = probe_subsets
         self.raw_dataset = raw_dataset
@@ -39,12 +44,24 @@ def __init__(self, raw_dataset, manifest, probe_subsets):
 
     @classmethod
     def methylated(cls, raw_dataset, manifest):
+        """ convenience method that feeds in a pre-defined list of methylated CpG locii probes """
         return cls(raw_dataset, manifest, METHYLATED_PROBE_SUBSETS)
 
     @classmethod
     def unmethylated(cls, raw_dataset, manifest):
+        """ convenience method that feeds in a pre-defined list of UNmethylated CpG locii probes """
         return cls(raw_dataset, manifest, UNMETHYLATED_PROBE_SUBSETS)
 
+    @classmethod
+    def snp_methylated(cls, raw_dataset, manifest):
+        """ convenience method that feeds in a pre-defined list of methylated Snp locii probes """
+        return cls(raw_dataset, manifest, METHYLATED_SNP_PROBES)
+
+    @classmethod
+    def snp_unmethylated(cls, raw_dataset, manifest):
+        """ convenience method that feeds in a pre-defined list of methylated Snp locii probes """
+        return cls(raw_dataset, manifest, UNMETHYLATED_SNP_PROBES)
+
     def build_data_frame(self):
         return pd.concat(self.data_frames.values())