Scribd
Upload
2 views11 pages

Geo Diff

The document outlines a workflow for analyzing kidney data using various R packages, focusing on background modeling and quality control. It involves fitting a Poisson background model, aggregating probes, and applying a background score test to filter targets based on expression levels. The workflow also includes estimating size factors and assessing signal quality to ensure reliable differential expression modeling.

Uploaded by

dnangelsheila
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
2 views11 pages

Geo Diff

The document outlines a workflow for analyzing kidney data using various R packages, focusing on background modeling and quality control. It involves fitting a Poisson background model, aggregating probes, and applying a background score test to filter targets based on expression levels. The workflow also includes estimating size factors and assessing signal quality to ensure reliable differential expression modeling.

Uploaded by

dnangelsheila
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

# workflow WTA

pacman::p_load(GeoDiff,dplyr,ggfortify,NanoStringNCTools,GeomxTools,Biobase,reshap
e2)

data('kidney')

kidney <- updateGeoMxSet(kidney)

kidney

head(pData(kidney))

table(pData(kidney)$`slide name`)

table(pData(kidney)$region)

# subset some ROIs

kidney <- kidney[,which(kidney$`slide name` %in%


c('disease3','normal3'))][,c(1:4,48:51,60:63,115:118)]

table(kidney$region,kidney$`slide name`)

table(kidney$`slide name`,kidney$class)

# background modeling

# Poisson background model using negative probes. The background model works on
the probe level data with all of the negative probes. Please do not use aggregateCounts
from GeomxTools before modeling.

featureType(kidney)

paste('## of negative probes:',sum(fData(kidney)$Negative))

# This model estimates a feature factor for each negative probe and a background size
factor for each ROI.

kidney <- fitPoisBG(kidney)

summary(pData(kidney)$sizefact)

summary(fData(kidney)$featfact[fData(kidney)$Negative])

# After running the model, we can diagnose it and see if there are any issues in the
dataset. One key metric for Poisson model is the dispersion. When dispersion is big, it is
called over-dispersion which often indicates batch effect or large outliers in the data.

[Link](123)

kidney_diag <- diagPoisBG(kidney)

notes(kidney_diag)$disper

# If the dispersion is >2, one of these factors might be present in the data. We can
check for outlier ROIs. People can choose to set outliers to be missing values and rerun
the Poisson Background model. Since the dispersion is within range here, the model will
not get run.

which(assayDataElement(kidney_diag,'low_outlier')==1,[Link] = T)

which(assayDataElement(kidney_diag,'up_outlier')==1,[Link] = T)

# Or if a batch effect is assumed, the poisson model can be adjusted to take different
groups into account. Here we are grouping the ROIs by slide.

kidney <- fitPoisBG(kidney,groupvar='slide name')

# The diagnosis of this model shows that when splitting by slide we similar results as
without splitting in this dataset.

# [Link](123)

# kidney_diag <- diagPoisBG(kidney,split=T)

# notes(kidney_diag)$disper

# Aggregate function

# After subsetting, we have a couple probes with 0 counts in all 16 ROIs so we will
remove them [Link] is a GeoDiff specific function for probe aggregation and
filtering. Probes get filtered based on either correlation and/or the score test within
targets and then aggregated. The negative probes do not get aggregated or filtered.

all0probeidx <- which(rowSums(exprs(kidney))==0)

if(length(all0probeidx)>0){kidney <- kidney[-all0probeidx,]} #减去全为 0 基因

kidney <- aggreprobe(kidney,use='cor')


# target QC

# score test

# Using the background score test, we can determine which targets are expressed
above the background of the negative probes across this dataset. We can then filter the
data to only targets above background, using a suggested pvalue threshold of 1e-3.

kidney <- BGScoreTest(kidney)

sum(fData(kidney)[['pvalues']]<1e-3,[Link] = T)

# For advanced users, there are 3 variables that can be changed in the score test. The
default for all three variables is FALSE. Any combination of these variables can be used.

# split - should the poisson background values split by group be used

# removeoutlier - should outlier negatives be removed

# useprior - use prior that the expression level of background follows a Beta
distribution, this will lead to a more conservative test but is prone to influence by
outliers

# kidneySplit <- BGScoreTest(kidney,split=T,removeoutlier=F,useprior=F) #16265

# kidneyOutliers <- BGScoreTest(kidney,split=F,removeoutlier=T,useprior=F) #16449

# kidneyPrior <- BGScoreTest(kidney,split=F,removeoutlier=F,useprior=T) #15662

# Estimate the size factor

# To estimate the signal size factor, we use the fit negative binomial threshold function.
This size factor represents technical variation between ROIs like sequencing depth. The
feature_high_fitNBth labeled genes are ones well above background that will be used in
later steps.

[Link](123)

kidney <- fitNBth(kidney,split=T)

features_high <- rownames(fData(kidney))[fData(kidney)$feature_high_fitNBth==1]

length(features_high)

# We can compare this threshold to the mean of the background as a sanity check.

bgMean <- mean(fData(kidney)$featfact,[Link]=T)


notes(kidney)[['threshold']]

bgMean

# This is a sanity check to see that the signal size factor and background size factor are
correlated but not redundant.

cor(kidney$sizefact,kidney$sizefact_fitNBth)

plot(kidney$sizefact, kidney$sizefact_fitNBth, xlab = "Background Size Factor",ylab =


"Signal Size Factor")

abline(a = 0, b = 1)

# In this dataset, this size factor correlate well with different quantiles, including 75%
quantile which is used in Q3 normalization.

# get only biological probes

posdat <- kidney[-which(fData(kidney)$CodeClass=='Negative'),]

posdat <- exprs(posdat)

quan <- sapply(c(0.75, 0.8, 0.9, 0.95), function(y) apply(posdat, 2, function(x) quantile(x,
probs = y)))

corrs <- apply(quan, 2, function(x) cor(x, kidney$sizefact_fitNBth))

names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs

quan75 <- apply(posdat,2,function(x) quantile(x,probs=.75))

# Quantile range (quantile - background size factor scaled by the mean feature factor
of negative probes) has better correlation with the signal size factor.

kidney <- QuanRange(kidney, split = FALSE, probs = c(0.75, 0.8, 0.9, 0.95))

corrs <- apply(pData(kidney)[, [Link](c(0.75, 0.8, 0.9, 0.95))], 2, function(x) cor(x,


kidney$sizefact_fitNBth))

names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs
# sample QC

# To filter out poor quality ROIs, we only keep those which have a high enough signal
in comparison to the background. In this dataset, all ROIs remain.

ROIs_hgh <- sampleNames(kidney)[which((quantile(fData(kidney)[['para']][,1],probs


= .9,[Link] = T) - notes(kidney)[['threshold']])*kidney$sizefact_fitNBth>2)]

features_all <- rownames(posdat)

# DE modeling

# Fixed Effect Model

# Running the DE model with default values.

# workflow WTA

pacman::p_load(GeoDiff,dplyr,ggfortify,NanoStringNCTools,GeomxTools,Biobase,reshap
e2)

data('kidney')

kidney <- updateGeoMxSet(kidney)

kidney

head(pData(kidney))

table(pData(kidney)$`slide name`)

table(pData(kidney)$region)

# subset some ROIs

kidney <- kidney[,which(kidney$`slide name` %in%


c('disease3','normal3'))][,c(1:4,48:51,60:63,115:118)]

table(kidney$region,kidney$`slide name`)

table(kidney$`slide name`,kidney$class)

# background modeling
# Poisson background model using negative probes. The background model works on
the probe level data with all of the negative probes. Please do not use aggregateCounts
from GeomxTools before modeling.

featureType(kidney)

paste('## of negative probes:',sum(fData(kidney)$Negative))

# This model estimates a feature factor for each negative probe and a background size
factor for each ROI.

kidney <- fitPoisBG(kidney)

summary(pData(kidney)$sizefact)

summary(fData(kidney)$featfact[fData(kidney)$Negative])

# After running the model, we can diagnose it and see if there are any issues in the
dataset. One key metric for Poisson model is the dispersion. When dispersion is big, it is
called over-dispersion which often indicates batch effect or large outliers in the data.

[Link](123)

kidney_diag <- diagPoisBG(kidney)

notes(kidney_diag)$disper

# If the dispersion is >2, one of these factors might be present in the data. We can
check for outlier ROIs. People can choose to set outliers to be missing values and rerun
the Poisson Background model. Since the dispersion is within range here, the model will
not get run.

which(assayDataElement(kidney_diag,'low_outlier')==1,[Link] = T)

which(assayDataElement(kidney_diag,'up_outlier')==1,[Link] = T)

# Or if a batch effect is assumed, the poisson model can be adjusted to take different
groups into account. Here we are grouping the ROIs by slide.

kidney <- fitPoisBG(kidney,groupvar='slide name')

# The diagnosis of this model shows that when splitting by slide we similar results as
without splitting in this dataset.

# [Link](123)

# kidney_diag <- diagPoisBG(kidney,split=T)

# notes(kidney_diag)$disper
# Aggregate function

# After subsetting, we have a couple probes with 0 counts in all 16 ROIs so we will
remove them [Link] is a GeoDiff specific function for probe aggregation and
filtering. Probes get filtered based on either correlation and/or the score test within
targets and then aggregated. The negative probes do not get aggregated or filtered.

all0probeidx <- which(rowSums(exprs(kidney))==0)

if(length(all0probeidx)>0){kidney <- kidney[-all0probeidx,]} #减去全为 0 基因

kidney <- aggreprobe(kidney,use='cor')

# target QC

# score test

# Using the background score test, we can determine which targets are expressed
above the background of the negative probes across this dataset. We can then filter the
data to only targets above background, using a suggested pvalue threshold of 1e-3.

kidney <- BGScoreTest(kidney)

sum(fData(kidney)[['pvalues']]<1e-3,[Link] = T)

# For advanced users, there are 3 variables that can be changed in the score test. The
default for all three variables is FALSE. Any combination of these variables can be used.

# split - should the poisson background values split by group be used

# removeoutlier - should outlier negatives be removed

# useprior - use prior that the expression level of background follows a Beta
distribution, this will lead to a more conservative test but is prone to influence by
outliers

# kidneySplit <- BGScoreTest(kidney,split=T,removeoutlier=F,useprior=F) #16265

# kidneyOutliers <- BGScoreTest(kidney,split=F,removeoutlier=T,useprior=F) #16449

# kidneyPrior <- BGScoreTest(kidney,split=F,removeoutlier=F,useprior=T) #15662

# Estimate the size factor


# To estimate the signal size factor, we use the fit negative binomial threshold function.
This size factor represents technical variation between ROIs like sequencing depth. The
feature_high_fitNBth labeled genes are ones well above background that will be used in
later steps.

[Link](123)

kidney <- fitNBth(kidney,split=T)

features_high <- rownames(fData(kidney))[fData(kidney)$feature_high_fitNBth==1]

length(features_high)

# We can compare this threshold to the mean of the background as a sanity check.

bgMean <- mean(fData(kidney)$featfact,[Link]=T)

notes(kidney)[['threshold']]

bgMean

# This is a sanity check to see that the signal size factor and background size factor are
correlated but not redundant.

cor(kidney$sizefact,kidney$sizefact_fitNBth)

plot(kidney$sizefact, kidney$sizefact_fitNBth, xlab = "Background Size Factor",ylab =


"Signal Size Factor")

abline(a = 0, b = 1)

# In this dataset, this size factor correlate well with different quantiles, including 75%
quantile which is used in Q3 normalization.

# get only biological probes

posdat <- kidney[-which(fData(kidney)$CodeClass=='Negative'),]

posdat <- exprs(posdat)

quan <- sapply(c(0.75, 0.8, 0.9, 0.95), function(y) apply(posdat, 2, function(x) quantile(x,
probs = y)))

corrs <- apply(quan, 2, function(x) cor(x, kidney$sizefact_fitNBth))

names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs
quan75 <- apply(posdat,2,function(x) quantile(x,probs=.75))

# Quantile range (quantile - background size factor scaled by the mean feature factor
of negative probes) has better correlation with the signal size factor.

kidney <- QuanRange(kidney, split = FALSE, probs = c(0.75, 0.8, 0.9, 0.95))

corrs <- apply(pData(kidney)[, [Link](c(0.75, 0.8, 0.9, 0.95))], 2, function(x) cor(x,


kidney$sizefact_fitNBth))

names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs

# sample QC

# To filter out poor quality ROIs, we only keep those which have a high enough signal
in comparison to the background. In this dataset, all ROIs remain.

ROIs_hgh <- sampleNames(kidney)[which((quantile(fData(kidney)[['para']][,1],probs


= .9,[Link] = T) - notes(kidney)[['threshold']])*kidney$sizefact_fitNBth>2)]

features_all <- rownames(posdat)

# DE modeling

# Fixed Effect Model

# Running the DE model with default values.

# Estimate the size factor

# To estimate the signal size factor, we use the fit negative binomial threshold function.
This size factor represents technical variation between ROIs like sequencing depth. The
feature_high_fitNBth labeled genes are ones well above background that will be used in
later steps.

[Link](123)

kidney <- fitNBth(kidney,split=T)

features_high <- rownames(fData(kidney))[fData(kidney)$feature_high_fitNBth==1]


length(features_high)

# We can compare this threshold to the mean of the background as a sanity check.

bgMean <- mean(fData(kidney)$featfact,[Link]=T)

notes(kidney)[['threshold']]

bgMean

# This is a sanity check to see that the signal size factor and background size factor are
correlated but not redundant.

cor(kidney$sizefact,kidney$sizefact_fitNBth)

plot(kidney$sizefact, kidney$sizefact_fitNBth, xlab = "Background Size Factor",ylab =


"Signal Size Factor")

abline(a = 0, b = 1)

# In this dataset, this size factor correlate well with different quantiles, including 75%
quantile which is used in Q3 normalization.

# get only biological probes

posdat <- kidney[-which(fData(kidney)$CodeClass=='Negative'),]

posdat <- exprs(posdat)

quan <- sapply(c(0.75, 0.8, 0.9, 0.95), function(y) apply(posdat, 2, function(x) quantile(x,
probs = y)))

corrs <- apply(quan, 2, function(x) cor(x, kidney$sizefact_fitNBth))

names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs

quan75 <- apply(posdat,2,function(x) quantile(x,probs=.75))

# Quantile range (quantile - background size factor scaled by the mean feature factor
of negative probes) has better correlation with the signal size factor.

kidney <- QuanRange(kidney, split = FALSE, probs = c(0.75, 0.8, 0.9, 0.95))

corrs <- apply(pData(kidney)[, [Link](c(0.75, 0.8, 0.9, 0.95))], 2, function(x) cor(x,


kidney$sizefact_fitNBth))
names(corrs) <- c(0.75, 0.8, 0.9, 0.95)

corrs

# sample QC

# To filter out poor quality ROIs, we only keep those which have a high enough signal
in comparison to the background. In this dataset, all ROIs remain.

ROIs_hgh <- sampleNames(kidney)[which((quantile(fData(kidney)[['para']][,1],probs


= .9,[Link] = T) - notes(kidney)[['threshold']])*kidney$sizefact_fitNBth>2)]

features_all <- rownames(posdat)

# DE modeling

# Fixed Effect Model

# Running the DE model with default values.

You might also like