python/pymethylprocess/PreProcessDataTypes.py
def preprocessMeffil(self, n_cores=6, n_pcs=4, qc_report_fname="qc/report.html", normalization_report_fname='norm/report.html', pc_plot_fname='qc/pc_plot.pdf', useCache=True, qc_only=True, qc_parameters={'p.beadnum.samples':0.1,'p.detection.samples':0.1,'p.detection.cpgs':0.1,'p.beadnum.cpgs':0.1}, rm_sex=False):
        """Run meffil preprocessing pipeline with functional normalization.

        Parameters
        ----------
        n_cores
            Number of CPUs to use.
        n_pcs
            Number of principal components to use for functional normalization, set to -1 to autoselect via kneedle algorithm.
        qc_report_fname
            HTML filename to store QC report.
        normalization_report_fname
            HTML filename to store normalization report
        pc_plot_fname
            PDF file to store principal components plot.
        useCache
            Use saved QC objects instead of running through QC again.
        qc_only
            Perform QC, then save and quit before normalization.
        qc_parameters
            Python dictionary with parameters for qc.
        rm_sex
            Remove non-autosomal cpgs?"""
        from pymethylprocess.meffil_functions import load_detection_p_values_beadnum, set_missing, remove_sex
        self.pheno = self.meffil.meffil_read_samplesheet(self.idat_dir, verbose=True)
        cache_storage_path = os.path.join(self.idat_dir,'QCObjects.rds')
        qc_parameters=robjects.r("""function(p.beadnum.samples,p.detection.samples,p.detection.cpgs,p.beadnum.cpgs,sex.outlier.sd){
                        qc.parameters <- meffil.qc.parameters(
                            	beadnum.samples.threshold             = p.beadnum.samples,
                            	detectionp.samples.threshold          = p.detection.samples,
                            	detectionp.cpgs.threshold             = p.detection.cpgs,
                            	beadnum.cpgs.threshold                = p.beadnum.cpgs,
                            	sex.outlier.sd                        = sex.outlier.sd,
                            	snp.concordance.threshold             = 0.95,
                            	sample.genotype.concordance.threshold = 0.8
                            )
                        return(qc.parameters)
                        }""")(qc_parameters['p.beadnum.samples'], qc_parameters['p.detection.samples'], qc_parameters['p.detection.cpgs'], qc_parameters['p.beadnum.cpgs'], qc_parameters['sex.outlier.sd'])
        if useCache:
            qc_list = robjects.r('readRDS')(cache_storage_path)
            qc_list = robjects.r("""function(qc.list,qc.parameters, qc.report.fname) {
                                qc.list$qc.summary <- meffil.qc.summary(qc.list$qc.objects,parameters=qc.parameters,verbose=F)
                                meffil.qc.report(qc.list$qc.summary, output.file=qc.report.fname)
                                return(qc.list)
                                }""")(qc_list,qc_parameters, qc_report_fname)
        else:
            qc_list = robjects.r("""function(samplesheet,n.cores,qc.parameters,qc.report.fname){
                qc.objects<-meffil.qc(samplesheet,mc.cores=n.cores,detection.threshold=0.000001,verbose=F)
                qc.summary<-meffil.qc.summary(qc.objects,parameters=qc.parameters,verbose=F)
                meffil.qc.report(qc.summary, output.file=qc.report.fname)
                return(list(qc.objects=qc.objects,qc.summary=qc.summary))
                }""")(self.pheno, n_cores, qc_parameters, qc_report_fname) # p.values=meffil.load.detection.pvalues(qc.objects)
        pc_df = robjects.r("""function(qc.list,pc.plot.fname){
            y <- meffil.plot.pc.fit(qc.list$qc.objects)
            print(y)
            ggsave(y$plot,filename=pc.plot.fname,height=6,width=6)
            return(y$data)
            }""")(qc_list,pc_plot_fname)
        print("Check QC report and select number of PCs. Will add option in future to adjust thresholds.")
        if n_pcs == -1:
            from kneed import KneeLocator
            pc_df=pandas2ri.ri2py(robjects.r['as'](pc_df,'data.frame'))
            pc_df['B'] = pc_df['U']+pc_df['M']
            knee=KneeLocator(pc_df['n'].values, pc_df['B'].values, S=1.0, curve='convex', direction='decreasing').knee
            if knee != None:
                n_pcs=int(knee)
            else:
                n_pcs=2#pc_df['n'].max() # try 2
            with open(pc_plot_fname.replace('.pdf','.txt'),'w') as f:
                f.write('pcs_selected:{}'.format(n_pcs))
        if qc_only:
            robjects.r('saveRDS')(qc_list,cache_storage_path)
            exit()
        if not os.path.exists(cache_storage_path):
            robjects.r('saveRDS')(qc_list,cache_storage_path)
        pval_beadnum=load_detection_p_values_beadnum(qc_list,n_cores)
        #print(pval_beadnum)
        self.beta_final = robjects.r("""function(qc.list, n.pcs, norm.report.fname,mc.cores) {
            options(mc.cores=mc.cores)
            qc.objects = qc.list$qc.objects
            qc.summary = qc.list$qc.summary
            outlier <- qc.summary$bad.samples
            if (nrow(outlier) > 0) {
            table(outlier$issue)
            index <- outlier$issue %in% c("Control probe (dye.bias)",
                                          "Methylated vs Unmethylated",
                                          "X-Y ratio outlier",
                                          "Low bead numbers",
                                          "Detection p-value",
                                          "Sex mismatch",
                                          "Genotype mismatch",
                                          "Control probe (bisulfite1)",
                                          "Control probe (bisulfite2)")
            outlier <- outlier[index,]
            if (nrow(outlier) > 0) {
                qc.objects <- meffil.remove.samples(qc.objects, outlier$sample.name)
            }
            }
            norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=n.pcs, verbose=F)
            norm <- meffil.normalize.samples(norm.objects, just.beta=F, cpglist.remove=qc.summary$bad.cpgs$name)
            beta <- meffil.get.beta(norm$M, norm$U)
            pcs <- meffil.methylation.pcs(beta)
            norm.summary <- meffil.normalization.summary(norm.objects, pcs=pcs)
            meffil.normalization.report(norm.summary, output.file=norm.report.fname)
            return(beta)}""")(qc_list, n_pcs, normalization_report_fname,n_cores)

        self.beta_final = set_missing(self.beta_final, pval_beadnum)

        if rm_sex:
            self.beta_final = remove_sex(self.beta_final)

        try:
            grdevice = importr("grDevices")
            geneplotter = importr("geneplotter")
            qr_report_fname=os.path.abspath(qc_report_fname).split('/')
            qr_report_fname[-1]='beta_dist.jpg'
            grdevice.jpeg('/'.join(qr_report_fname),height=900,width=600)
            base.par(mfrow=robjects.vectors.IntVector([1,2]))
            self.enmix.multidensity(self.beta_final, main="Multidensity")
            self.enmix.multifreqpoly(self.beta_final, xlab="Beta value")
            grdevice.dev_off()
        except:
            pass
Run meffil preprocessing pipeline with functional normalization.

        Parameters
        ----------
        n_cores
            Number of CPUs to use.
        n_pcs
            Number of principal components to use for functional normalization, set to -1 to autoselect via kneedle algorithm.
        qc_report_fname
            HTML filename to store QC report.
        normalization_report_fname
            HTML filename to store normalization report
        pc_plot_fname
            PDF file to store principal components plot.
        useCache
            Use saved QC objects instead of running through QC again.
        qc_only
            Perform QC, then save and quit before normalization.
        qc_parameters
            Python dictionary with parameters for qc.
        rm_sex
            Remove non-autosomal cpgs?
Run meffil preprocessing pipeline with functional normalization.
Parameters


Run meffil preprocessing pipeline with functional normalization . Parameters
def preprocessMeffil(self, n_cores=6, n_pcs=4, qc_report_fname="qc/report.html", normalization_report_fname='norm/report.html', pc_plot_fname='qc/pc_plot.pdf', useCache=True, qc_only=True, qc_parameters={'p.beadnum.samples':0.1,'p.detection.samples':0.1,'p.detection.cpgs':0.1,'p.beadnum.cpgs':0.1}, rm_sex=False):
        from pymethylprocess.meffil_functions import load_detection_p_values_beadnum, set_missing, remove_sex
        self.pheno = self.meffil.meffil_read_samplesheet(self.idat_dir, verbose=True)
        cache_storage_path = os.path.join(self.idat_dir,'QCObjects.rds')
        qc_parameters=robjects.r("""function(p.beadnum.samples,p.detection.samples,p.detection.cpgs,p.beadnum.cpgs,sex.outlier.sd){
                        qc.parameters <- meffil.qc.parameters(
                            	beadnum.samples.threshold             = p.beadnum.samples,
                            	detectionp.samples.threshold          = p.detection.samples,
                            	detectionp.cpgs.threshold             = p.detection.cpgs,
                            	beadnum.cpgs.threshold                = p.beadnum.cpgs,
                            	sex.outlier.sd                        = sex.outlier.sd,
                            	snp.concordance.threshold             = 0.95,
                            	sample.genotype.concordance.threshold = 0.8
                            )
                        return(qc.parameters)
                        }""")(qc_parameters['p.beadnum.samples'], qc_parameters['p.detection.samples'], qc_parameters['p.detection.cpgs'], qc_parameters['p.beadnum.cpgs'], qc_parameters['sex.outlier.sd'])
        if useCache:
            qc_list = robjects.r('readRDS')(cache_storage_path)
            qc_list = robjects.r("""function(qc.list,qc.parameters, qc.report.fname) {
                                qc.list$qc.summary <- meffil.qc.summary(qc.list$qc.objects,parameters=qc.parameters,verbose=F)
                                meffil.qc.report(qc.list$qc.summary, output.file=qc.report.fname)
                                return(qc.list)
                                }""")(qc_list,qc_parameters, qc_report_fname)
        else:
            qc_list = robjects.r("""function(samplesheet,n.cores,qc.parameters,qc.report.fname){
                qc.objects<-meffil.qc(samplesheet,mc.cores=n.cores,detection.threshold=0.000001,verbose=F)
                qc.summary<-meffil.qc.summary(qc.objects,parameters=qc.parameters,verbose=F)
                meffil.qc.report(qc.summary, output.file=qc.report.fname)
                return(list(qc.objects=qc.objects,qc.summary=qc.summary))
                }""")(self.pheno, n_cores, qc_parameters, qc_report_fname) 
        pc_df = robjects.r("""function(qc.list,pc.plot.fname){
            y <- meffil.plot.pc.fit(qc.list$qc.objects)
            print(y)
            ggsave(y$plot,filename=pc.plot.fname,height=6,width=6)
            return(y$data)
            }""")(qc_list,pc_plot_fname)
        print("Check QC report and select number of PCs. Will add option in future to adjust thresholds.")
        if n_pcs == -1:
            from kneed import KneeLocator
            pc_df=pandas2ri.ri2py(robjects.r['as'](pc_df,'data.frame'))
            pc_df['B'] = pc_df['U']+pc_df['M']
            knee=KneeLocator(pc_df['n'].values, pc_df['B'].values, S=1.0, curve='convex', direction='decreasing').knee
            if knee != None:
                n_pcs=int(knee)
            else:
                n_pcs=2
            with open(pc_plot_fname.replace('.pdf','.txt'),'w') as f:
                f.write('pcs_selected:{}'.format(n_pcs))
        if qc_only:
            robjects.r('saveRDS')(qc_list,cache_storage_path)
            exit()
        if not os.path.exists(cache_storage_path):
            robjects.r('saveRDS')(qc_list,cache_storage_path)
        pval_beadnum=load_detection_p_values_beadnum(qc_list,n_cores)
        self.beta_final = robjects.r("""function(qc.list, n.pcs, norm.report.fname,mc.cores) {
            options(mc.cores=mc.cores)
            qc.objects = qc.list$qc.objects
            qc.summary = qc.list$qc.summary
            outlier <- qc.summary$bad.samples
            if (nrow(outlier) > 0) {
            table(outlier$issue)
            index <- outlier$issue %in% c("Control probe (dye.bias)",
                                          "Methylated vs Unmethylated",
                                          "X-Y ratio outlier",
                                          "Low bead numbers",
                                          "Detection p-value",
                                          "Sex mismatch",
                                          "Genotype mismatch",
                                          "Control probe (bisulfite1)",
                                          "Control probe (bisulfite2)")
            outlier <- outlier[index,]
            if (nrow(outlier) > 0) {
                qc.objects <- meffil.remove.samples(qc.objects, outlier$sample.name)
            }
            }
            norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=n.pcs, verbose=F)
            norm <- meffil.normalize.samples(norm.objects, just.beta=F, cpglist.remove=qc.summary$bad.cpgs$name)
            beta <- meffil.get.beta(norm$M, norm$U)
            pcs <- meffil.methylation.pcs(beta)
            norm.summary <- meffil.normalization.summary(norm.objects, pcs=pcs)
            meffil.normalization.report(norm.summary, output.file=norm.report.fname)
            return(beta)}""")(qc_list, n_pcs, normalization_report_fname,n_cores)
        self.beta_final = set_missing(self.beta_final, pval_beadnum)
        if rm_sex:
            self.beta_final = remove_sex(self.beta_final)
        try:
            grdevice = importr("grDevices")
            geneplotter = importr("geneplotter")
            qr_report_fname=os.path.abspath(qc_report_fname).split('/')
            qr_report_fname[-1]='beta_dist.jpg'
            grdevice.jpeg('/'.join(qr_report_fname),height=900,width=600)
            base.par(mfrow=robjects.vectors.IntVector([1,2]))
            self.enmix.multidensity(self.beta_final, main="Multidensity")
            self.enmix.multifreqpoly(self.beta_final, xlab="Beta value")
            grdevice.dev_off()
        except:
            pass
def preprocessMeffil ( self , n_cores = 6 , n_pcs = 4 , qc_report_fname = "qc/report.html" , normalization_report_fname = 'norm/report.html' , pc_plot_fname = 'qc/pc_plot.pdf' , useCache = True , qc_only = True , qc_parameters = { 'p.beadnum.samples' : 0.1 , 'p.detection.samples' : 0.1 , 'p.detection.cpgs' : 0.1 , 'p.beadnum.cpgs' : 0.1 } , rm_sex = False ) : from pymethylprocess . meffil_functions import load_detection_p_values_beadnum , set_missing , remove_sex self . pheno = self . meffil . meffil_read_samplesheet ( self . idat_dir , verbose = True ) cache_storage_path = os . path . join ( self . idat_dir , 'QCObjects.rds' ) qc_parameters = robjects . r ( """function(p.beadnum.samples,p.detection.samples,p.detection.cpgs,p.beadnum.cpgs,sex.outlier.sd){
                        qc.parameters <- meffil.qc.parameters(
                            	beadnum.samples.threshold             = p.beadnum.samples,
                            	detectionp.samples.threshold          = p.detection.samples,
                            	detectionp.cpgs.threshold             = p.detection.cpgs,
                            	beadnum.cpgs.threshold                = p.beadnum.cpgs,
                            	sex.outlier.sd                        = sex.outlier.sd,
                            	snp.concordance.threshold             = 0.95,
                            	sample.genotype.concordance.threshold = 0.8
                            )
                        return(qc.parameters)
                        }""" ) ( qc_parameters [ 'p.beadnum.samples' ] , qc_parameters [ 'p.detection.samples' ] , qc_parameters [ 'p.detection.cpgs' ] , qc_parameters [ 'p.beadnum.cpgs' ] , qc_parameters [ 'sex.outlier.sd' ] ) if useCache : qc_list = robjects . r ( 'readRDS' ) ( cache_storage_path ) qc_list = robjects . r ( """function(qc.list,qc.parameters, qc.report.fname) {
                                qc.list$qc.summary <- meffil.qc.summary(qc.list$qc.objects,parameters=qc.parameters,verbose=F)
                                meffil.qc.report(qc.list$qc.summary, output.file=qc.report.fname)
                                return(qc.list)
                                }""" ) ( qc_list , qc_parameters , qc_report_fname ) else : qc_list = robjects . r ( """function(samplesheet,n.cores,qc.parameters,qc.report.fname){
                qc.objects<-meffil.qc(samplesheet,mc.cores=n.cores,detection.threshold=0.000001,verbose=F)
                qc.summary<-meffil.qc.summary(qc.objects,parameters=qc.parameters,verbose=F)
                meffil.qc.report(qc.summary, output.file=qc.report.fname)
                return(list(qc.objects=qc.objects,qc.summary=qc.summary))
                }""" ) ( self . pheno , n_cores , qc_parameters , qc_report_fname ) pc_df = robjects . r ( """function(qc.list,pc.plot.fname){
            y <- meffil.plot.pc.fit(qc.list$qc.objects)
            print(y)
            ggsave(y$plot,filename=pc.plot.fname,height=6,width=6)
            return(y$data)
            }""" ) ( qc_list , pc_plot_fname ) print ( "Check QC report and select number of PCs. Will add option in future to adjust thresholds." ) if n_pcs == - 1 : from kneed import KneeLocator pc_df = pandas2ri . ri2py ( robjects . r [ 'as' ] ( pc_df , 'data.frame' ) ) pc_df [ 'B' ] = pc_df [ 'U' ] + pc_df [ 'M' ] knee = KneeLocator ( pc_df [ 'n' ] . values , pc_df [ 'B' ] . values , S = 1.0 , curve = 'convex' , direction = 'decreasing' ) . knee if knee != None : n_pcs = int ( knee ) else : n_pcs = 2 with open ( pc_plot_fname . replace ( '.pdf' , '.txt' ) , 'w' ) as f : f . write ( 'pcs_selected:{}' . format ( n_pcs ) ) if qc_only : robjects . r ( 'saveRDS' ) ( qc_list , cache_storage_path ) exit ( ) if not os . path . exists ( cache_storage_path ) : robjects . r ( 'saveRDS' ) ( qc_list , cache_storage_path ) pval_beadnum = load_detection_p_values_beadnum ( qc_list , n_cores ) self . beta_final = robjects . r ( """function(qc.list, n.pcs, norm.report.fname,mc.cores) {
            options(mc.cores=mc.cores)
            qc.objects = qc.list$qc.objects
            qc.summary = qc.list$qc.summary
            outlier <- qc.summary$bad.samples
            if (nrow(outlier) > 0) {
            table(outlier$issue)
            index <- outlier$issue %in% c("Control probe (dye.bias)",
                                          "Methylated vs Unmethylated",
                                          "X-Y ratio outlier",
                                          "Low bead numbers",
                                          "Detection p-value",
                                          "Sex mismatch",
                                          "Genotype mismatch",
                                          "Control probe (bisulfite1)",
                                          "Control probe (bisulfite2)")
            outlier <- outlier[index,]
            if (nrow(outlier) > 0) {
                qc.objects <- meffil.remove.samples(qc.objects, outlier$sample.name)
            }
            }
            norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=n.pcs, verbose=F)
            norm <- meffil.normalize.samples(norm.objects, just.beta=F, cpglist.remove=qc.summary$bad.cpgs$name)
            beta <- meffil.get.beta(norm$M, norm$U)
            pcs <- meffil.methylation.pcs(beta)
            norm.summary <- meffil.normalization.summary(norm.objects, pcs=pcs)
            meffil.normalization.report(norm.summary, output.file=norm.report.fname)
            return(beta)}""" ) ( qc_list , n_pcs , normalization_report_fname , n_cores ) self . beta_final = set_missing ( self . beta_final , pval_beadnum ) if rm_sex : self . beta_final = remove_sex ( self . beta_final ) try : grdevice = importr ( "grDevices" ) geneplotter = importr ( "geneplotter" ) qr_report_fname = os . path . abspath ( qc_report_fname ) . split ( '/' ) qr_report_fname [ - 1 ] = 'beta_dist.jpg' grdevice . jpeg ( '/' . join ( qr_report_fname ) , height = 900 , width = 600 ) base . par ( mfrow = robjects . vectors . IntVector ( [ 1 , 2 ] ) ) self . enmix . multidensity ( self . beta_final , main = "Multidensity" ) self . enmix . multifreqpoly ( self . beta_final , xlab = "Beta value" ) grdevice . dev_off ( ) except : pass
Run meffil preprocessing pipeline with functional normalization.
Run meffil preprocessing pipeline with functional normalization .
"""Run meffil preprocessing pipeline with functional normalization.

        Parameters
        ----------
        n_cores
            Number of CPUs to use.
        n_pcs
            Number of principal components to use for functional normalization, set to -1 to autoselect via kneedle algorithm.
        qc_report_fname
            HTML filename to store QC report.
        normalization_report_fname
            HTML filename to store normalization report
        pc_plot_fname
            PDF file to store principal components plot.
        useCache
            Use saved QC objects instead of running through QC again.
        qc_only
            Perform QC, then save and quit before normalization.
        qc_parameters
            Python dictionary with parameters for qc.
        rm_sex
            Remove non-autosomal cpgs?""" # p.values=meffil.load.detection.pvalues(qc.objects) #pc_df['n'].max() # try 2 #print(pval_beadnum)
self
n_cores
n_pcs
qc_report_fname
normalization_report_fname
pc_plot_fname
useCache
qc_only
qc_parameters
rm_sex