2  GA-adjusted Z-score

The gestational windows of blood sampling were ~12wkGA, ~20wkGA, ~28wkGA and ~36wkGA. However, at each time point, there was a degree of variation in the exact gestational age. Therefore we calculated a GA-adjusted z score for each protein concentration within each time window, referent to the whole of the random sub-cohort (or a comparator group).

2.1 Data Import

#####################
# Meta Information ##
#####################
dt.sample<-readxl::read_excel("OLINK-SAMPLEMANIFEST.xlsx") %>% data.table 
setnames(dt.sample,"Unique_Sample_ID","SampleID") # Unique_Sample_ID => SampleID
dt.sample[,Sampletakenat:=gsub(" weeks","wk",Sampletakenat)] 

## phenotype information
dt.meta<-openxlsx::read.xlsx("PHEN_META_INFO.xlsx") %>% data.table
dt.GA<-dt.meta[,.(POPSID,Comparator,`12wk`=pc_ga1,`20wk`=pc_ga2,`28wk`=pc_ga3,`36wk`=pc_ga4)]

######################
# NPX raw data file ##
######################
dt.olink.raw<-fread("OLINK_RAW_DATA.csv.gz")

# Now, update the input table after QC
dt.olink.input<-dt.olink.raw[Sample_Type=="SAMPLE" & QC_Warning!="EXCLUDED" & Assay_Warning!="EXCLUDED"]

## add POPSID, and GA info and split by GA
dt.olink<-merge(dt.olink.input, dt.sample[,.(SampleID,POPSID,GA=Sampletakenat)])
dl.olink<-split(dt.olink, dt.olink$GA) # by each GA

2.2 Normalisation

################
# for each GA ##
################
dl.npxz.GA<-lapply(dl.olink, function(dt.olink.GA){
            my.GA<-dt.olink.GA[,unique(GA)]
            this.col<-c("POPSID",my.GA)
            dt.com.GA<-dt.GA[Comparator==1,..this.col] # get the GA of Comparator group of this GA
            setnames(dt.com.GA,c("POPSID","GA")) # make the name of a given GA neutral

            # get the NPX of Comparator for a given GA
            dt.com.npx<-merge(dt.com.GA,dt.olink.GA[,.(POPSID,OlinkID,NPX)]) 
            dl.com.npx<-split(dt.com.npx, dt.com.npx$OlinkID) # split by OlinkID
            dl.com.npx  %>% length # should be 2904

            # Now, for each of the NPX data.table of Comparator group for a given OlinkID, GA, and disease type
            dt.npxz.GA <-lapply(dl.com.npx, function(dt.com.npx.GA) {
              my.olink <- dt.com.npx.GA[,unique(OlinkID)] # set this olinkid
              my.fit<-lm(NPX ~ GA , data=dt.com.npx.GA) # a linear regress model of NPX by the function of GA

              # get the NPX data.table of a target disease for a given OlinkID and GA 
              dt.this.olink<-dt.olink.GA[OlinkID==my.olink][,.(POPSID,OlinkID,NPX)] 

              # check the number of outcome from this data.table and compared with that of dl.npx - should be the same
              # get the exact GA for this NPX data.table
              dt.this.olink.GA<-merge(dt.this.olink, dt.GA[,..this.col])
              setnames(dt.this.olink.GA,my.GA, "GA") # make the name of a given GA neutral

              # get the residual (i.e. observed - predicted) of this fitted model 
              # the following equation is the formular to calculate standardised residual.
              # NB, `my.fit` from the Comparator group, and `my.target.fit` from the target disease
              # 1) `summary(my.fit)$sigma` is the RSE (residual standard error of the model)
              # 2) lvg: the leverage of this observation (see https://en.wikipedia.org/wiki/Leverage_(statistics))
              resid<-dt.this.olink.GA$NPX - predict(my.fit, newdata=dt.this.olink.GA) 
              my.target.fit <-lm(NPX ~ GA , data=dt.this.olink.GA, na.action="na.exclude")
              lvg<-hatvalues(my.target.fit) 
              std.resid=resid/(summary(my.fit)$sigma*(sqrt(1-lvg)))

              dt.this.olink.GA[,NPXZ:=std.resid][,.(POPSID,OlinkID,NPX,NPXZ)]
            }) %>% rbindlist # end of for each OlinkID
            # long to wide format
            dt.npxz.GA %>% dcast.data.table(POPSID~OlinkID,value.var="NPXZ")
          }) # end of for each GA

2.3 Winsorisation

In order to limit the influence of extreme values, all z scores were winsorised -5 and +5.

my.sd<-5
dl.npxz5<-lapply(dl.npxz.GA, function(dt.mat){
  # wide-to-long
  dt.bar<- dt.mat %>% melt.data.table(id.var="POPSID",variable.name="OlinkID",value.name="NPXZ")
  dt.bar[,NPXZ5:=ifelse(NPXZ > my.sd,my.sd,ifelse(NPXZ < -my.sd, -my.sd, NPXZ))]
  dcast.data.table(dt.bar, POPSID~OlinkID, value.var="NPXZ5")
})