<jupyter_start><jupyter_text>Data processing and statistics demo<NAME>, PhD Bioinformatician/Computational Biologist Bioinformatics Core Research Facility, Center for Biotechnology Data Life Science Core, NCIBC [<EMAIL>] BEAD E206, UNL, 68588-0665 * Assess the quality of data: missing values* Assess feature presence/absence: missing value imputation* High-dimensional data processing    + Log2 transformation    + Normalization    + Scaling* Univariate analysis    + T-test    + Wilcoxon rank sum test* Multivariate analysis    + PCA    + PLS-DA**Goal: feature selection/extraction** and visualization.<jupyter_code>suppressPackageStartupMessages(library("tidyverse"))<jupyter_output>Warning message:
"package 'tidyverse' was built under R version 4.0.5"
Warning message:
"package 'ggplot2' was built under R version 4.0.5"
Warning message:
"package 'tibble' was built under R version 4.0.5"
Warning message:
"package 'dplyr' was built under R version 4.0.5"
Warning message:
"package 'forcats' was built under R version 4.0.5"<jupyter_text>Set working directory where the data is.<jupyter_code>setwd("./demo/")<jupyter_output><empty_output><jupyter_text>Import data and data wrangling.<jupyter_code>fh = read.csv(file = "malaria_feature_table.csv", header = TRUE, row.names = 1, stringsAsFactors = FALSE) %>%
    as.data.frame()

head(fh)
dim(fh)<jupyter_output><empty_output><jupyter_text>This dataset has 12 samples and 5114 features. If your samples were not organized by group, this string of code can take care of it. We also subset the labels of each sample out as a seperated matrix for future use.<jupyter_code>fh <- fh[,order(fh["Label",])] ## order samples by group
labels_d1 <- as.matrix(fh[1,])<jupyter_output><empty_output><jupyter_text>Check what we have in our global environment.<jupyter_code>ls()<jupyter_output><empty_output><jupyter_text>Check missing data.<jupyter_code>d1 <- fh[-1,] # removing label from the dataset, just data left. 

sum(d1==0) ## Total missing values
sum(d1 ==0)/(dim(d1)[1]*dim(d1)[2])<jupyter_output><empty_output><jupyter_text>The dataset contains 10448 missing values, missing data percentage is 17%, which is acceptable for raw peak intensity table. We want to convert all zeros to NA.<jupyter_code>d1[d1 == 0] <- NA ## Change zeros to NA
sum(is.na(d1))<jupyter_output><empty_output><jupyter_text>Access feature presence/absence. As we stated, a feature is deemed present when missing values of a feature accounts for less than 50% of sample space, all sample or within each group.<jupyter_code>d1$nacount <- as.matrix(apply(d1, 1, function(x) sum(is.na(x)))) #count the missing signal values by feature
d1$nacount_Native = as.matrix(apply(d1[, 1:6], 1, function(x) sum(is.na(x))))
d1$nacount_Semi = as.matrix(apply(d1[, 7:12], 1, function(x) sum(is.na(x))))

head(d1)<jupyter_output><empty_output><jupyter_text>Scenario 1, keep a feature when it has na less than 50% of all samples.<jupyter_code>d2 <- d1[which(d1$nacount < round((length(d1)-3)*.5)), ] %>%
        select(-c(nacount, nacount_Semi, nacount_Native))
sum(is.na(d2)) # number of missing values
sum(is.na(d2))/(dim(d2)[1]*dim(d2)[2]) # percentage missing data<jupyter_output><empty_output><jupyter_text>Scenario 2, keep a feature when it has na less than 50% of either group.<jupyter_code>d2_2 = d1[which(d1$nacount_Native < round(6*.5) & d1$nacount_Semi < round(6*.5)), ] %>%
        select(-c(nacount, nacount_Semi, nacount_Native))
sum(is.na(d2_2)) # number of missing values
sum(is.na(d2_2))/(dim(d2_2)[1]*dim(d2_2)[2]) # percentage missing data<jupyter_output><empty_output><jupyter_text>Missing data imputation.We are gonna impute both datasets, because their missing value percentages are both low enough. We will use k-nearest neighbors for imputation. This is how we decide k, square root of number of features.<jupyter_code># decide K (how many neighbor to use)
k <- d2_2 %>% length() %>% sqrt() # square root of number of features
if (k>=3){
        if (is.integer(k)) k = k else k = ceiling(k)
} else {k = 3}

k<jupyter_output><empty_output><jupyter_text>Scenario 1:<jupyter_code>imputed_meta_d1 <- impute::impute.knn(as.matrix(d2), k = k, rowmax = 0.5,
                              colmax = 0.8, maxp = 1500)

imputed_d1 <- as.data.frame(imputed_meta_d1[["data"]], colnames = TRUE)
sum(imputed_d1 == 0)
sum(imputed_d1 < 0)<jupyter_output>Cluster size 4291 broken into 39 4252 
Done cluster 39 
Cluster size 4252 broken into 196 4056 
Done cluster 196 
Cluster size 4056 broken into 340 3716 
Done cluster 340 
Cluster size 3716 broken into 3115 601 
Cluster size 3115 broken into 827 2288 
Done cluster 827 
Cluster size 2288 broken into 3 2285 
Done cluster 3 
Cluster size 2285 broken into 1855 430 
Cluster size 1855 broken into 501 1354 
Done cluster 501 
Done cluster 1354 
Done cluster 1855 
Done cluster 430 
Done cluster 2285 
Done cluster 2288 
Done cluster 3115 
Done cluster 601 
Done cluster 3716 
Done cluster 4056 
Done cluster 4252<jupyter_text>There is 8 zeros after imputation. This usually means there were too much missing values to begin with, and we are pushing the limit of KNN. Scenatio 2:<jupyter_code>imputed_meta_d1 <- impute::impute.knn(as.matrix(d2_2), k = k, rowmax = 0.5,
                              colmax = 0.8, maxp = 1500)

imputed_d1 <- as.data.frame(imputed_meta_d1[["data"]], colnames = TRUE)
sum(imputed_d1 == 0)
sum(imputed_d1 < 0)<jupyter_output>Cluster size 3826 broken into 3776 50 
Cluster size 3776 broken into 3599 177 
Cluster size 3599 broken into 3266 333 
Cluster size 3266 broken into 572 2694 
Done cluster 572 
Cluster size 2694 broken into 2083 611 
Cluster size 2083 broken into 712 1371 
Done cluster 712 
Done cluster 1371 
Done cluster 2083 
Done cluster 611 
Done cluster 2694 
Done cluster 3266 
Done cluster 333 
Done cluster 3599 
Done cluster 177 
Done cluster 3776 
Done cluster 50<jupyter_text>There is no missing values after imputation for scenario 2, so this is the data we will continue with. Check the data distribution after imputation.<jupyter_code>MASS::truehist(unlist(imputed_d1), main = "Imputed", nbins = 50) # check data distribution<jupyter_output><empty_output><jupyter_text>This distribution is very typical for mass spec-generated data. Next step is data processing.  Log2 transformation.<jupyter_code>log2_d1 <- log2(type.convert(imputed_d1))
MASS::truehist(unlist(log2_d1), main = "log2 transformed", nbins = 50) # check data distribution
car::qqPlot(unlist(log2_d1), main = "log2 transformed",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")<jupyter_output><empty_output><jupyter_text>This already improves the data distribution quite a lot. But for the sake of demonstration, we will proceed with eigenMS normalization, which is implemented using package `ProteoMM`. Again, `ProteoMM` is a bioconductor package.  EigenMS normalization<jupyter_code>if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("ProteoMM", quietly = TRUE))
    BiocManager::install("ProteoMM")

suppressPackageStartupMessages(library("ProteoMM"))
d4 <- cbind.data.frame(row.names(log2_d1), log2_d1) # check data structure, feature by sample
head(d4)<jupyter_output><empty_output><jupyter_text>This following section of code is based on `ProteoMM` documentation. You can browse package vignette by typing this line of code in R console:`browseVignettes("ProteoMM")`<jupyter_code>rowmeta <- d4[1,]
rowmeta2 <- as.data.frame(t(colnames(d4)))
colmeta <- d4[-1,1]

intsCols <- 2:dim(d4)[2]
m_logInts = ProteoMM::make_intencities(d4, intsCols)

metaCols <- 1
m_mets.info = as.matrix(ProteoMM::make_meta(d4, metaCols))

grps = as.factor(t(labels_d1))
m_nummiss = sum(is.na(m_logInts)) #check zeros

set.seed(42)
hs_m_ints_eig1 = eig_norm1(m=m_logInts,treatment = grps,
                           prot.info=m_mets.info)
hs_m_ints_norm_1bt = eig_norm2(rv=hs_m_ints_eig1)<jupyter_output>Warning message in eig_norm1(m = m_logInts, treatment = grps, prot.info = m_mets.info):
"This function uses random namber generator. For reproducibility use
    set.seed(12345) with your choce of parameter"
Data dimentions: 382612

Treatment groups: NaiveNaiveNaiveNaiveNaiveNaiveSemi_immueSemi_immueSemi_immueSemi_immueSemi_immueSemi_immue

Selecting complete peptides

Got 2+ treatment grps

Computing SVD, estimating Eigentrends...

Number of treatments: 2

Number of complete peptides (and samples) used in SVD

382612

Number of treatment groups (in svd.id): 2

Starting Bootstrap.....

Iteration 50

Iteration 100

Iteration 150

Iteration 200

Iteration 250

Iteration 300

Iteration 350

Iteration 400

Iteration 450

Iteration 500

Number of bias trends automatically detected 5

Preparing to plot...

Unique number of treatment combinations:2

Normalizing...

Processing peptide 100

Processing peptide 200

Processing peptide 300

Processing peptide 400

Processing peptide 500

Processing[...]<jupyter_text>The dataframe after eigenMS normalization is named `norm_d1`, and we will check its distribution as well.<jupyter_code>norm_d1 <- as.data.frame(hs_m_ints_norm_1bt[["norm_m"]])
par(mfrow=c(1,1)) # after last step, the R graphic space are devided to 3 by 2 grid, this step is simply to restore. 
MASS::truehist(unlist(t(norm_d1)), main = "EigenMS normalization", nbins = 50)
car::qqPlot(unlist(t(norm_d1)),main = "log2-EigenMS",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")<jupyter_output><empty_output><jupyter_text>This data distribution looks normal based on histogram and qq plot, so we will stop here and proceed with univariate and multivariate analysis. Other normalization methods include rank invariant normalization implemented using `lumi`, cubic spline normalization implemented using `affy`. If normalization step cannot transform the data to a near Gaussian distribution, we will follow up with scaling step, implemented using `MetaboAnalyze`. Univariate analysis, t-test and Wilcoxon rank sum test.<jupyter_code>d5 <- rbind.data.frame(labels_d1, norm_d1)

C <- d5[-1,which(d5[1,]=="Naive")]; P <- d5[-1,which(d5[1,]=="Semi_immue")]

wt <- data.frame(Doubles=double(),Doubles=double(),Doubles=double(),Doubles=double(),Doubles=double())

# Calculate p-values from both t-test and Wilcoson rank sum test. 
# P-values will need to be adjusted with Benjamini-Hochberg correction (i.e., FDR)
for (i in 1:dim(norm_d1)[1]){
        wt[i,1] <- t.test(type.convert(t(P[i,])),type.convert(t(C[i,])))$p.value
        wt[i,3] <- wilcox.test(type.convert(t(P[i,])),type.convert(t(C[i,])))$p.value
        az1 <- 2^(type.convert(C[i,])); az1L <- length(az1)
        az2 <- 2^(type.convert(P[i,])); az2L <- length(az2)
        wt[i,5] <- gtools::foldchange(sum(az2)/az2L, sum(az1)/az1L)
        #print(az1)
}
wt[,2] <- as.data.frame(p.adjust(wt[,1], method = 'BH'))#, n = length(p))
wt[,4] <- as.data.frame(p.adjust(wt[,3], method = 'BH'))#, n = length(p))

colnames(wt) <- c("pT", "BHT", "pW", "BHW", "FC(lin)")
row.names(wt) <- rownames(norm_d1)

# Calculate log2FoldChange
log2FoldChange =gtools::foldchange2logratio(wt$`FC(lin)`,base = 2)
wt = cbind(wt,log2FoldChange)
head(wt)<jupyter_output><empty_output><jupyter_text>This is the statistical result table after univariate analysis, BH-corrected t-test and Wilcoxon rank sum tests. Next step is to screen for features that passed BHT < 0.05 & BHW < 0.05 criteria.<jupyter_code>sum(wt[,1] < 0.05 & wt[,3] < 0.05)
sum(wt[,2] < 0.01 & wt[,4] < 0.01)
sig_wt <- wt[which(wt[,2] < 0.01 & wt[,4] < 0.01),]
head(sig_wt)
dim(sig_wt)<jupyter_output><empty_output><jupyter_text>There are 2448 features that are significant different when comparing the two groups. Multivariate analysis, PCA<jupyter_code>if (!requireNamespace("pcaMethods", quietly = TRUE))
    BiocManager::install("pcaMethods")

suppressPackageStartupMessages(library("pcaMethods"))
pc = pcaMethods::pca(t(norm_d1), nPcs = 9, method = "svd", scale = "pareto", center = TRUE)

pca = rbind.data.frame(labels_d1, norm_d1) %>%
        t(.) %>%
        as.data.frame(.) %>%
        rownames_to_column(.) %>%
        as_tibble(.) %>%
        select(rowname, Label) %>%
        column_to_rownames(.,var = "rowname") %>%
        merge(., pcaMethods::scores(pc), by=0)<jupyter_output><empty_output><jupyter_text>Score plot<jupyter_code>ggplot(pca, aes(PC1, PC2,colour=Label)) + geom_point() + stat_ellipse() +
        xlab(paste("PC1", round(pc@R2[1] * 100, digits = 1), "% of the variance")) +
        ylab(paste("PC2", round(pc@R2[2] * 100, digits = 1), "% of the variance"))<jupyter_output><empty_output><jupyter_text>As is shown in score plot, there are clear separation between these two groups, and the first PC accounts for ~74% of total variance of the dataset. The ellipses represent 95% confidence intervals.  Scree plot<jupyter_code>var_explained_pc = data.frame(PC = paste0("PC", 1:9), var_explained = (pc@sDev)^2/sum((pc@sDev)^2), eigenvalues = pc@sDev)
var_explained_pc %>%
        ggplot(aes(x=PC,y=eigenvalues, group=1))+
        geom_point(size=4)+
        geom_line()+
        labs(title="Scree plot: PCA on scaled data")<jupyter_output><empty_output><jupyter_text>As PC1 accounts for majority of the variance in the dataset, we should take a look at what features count heavier in this PC.  Loading plot<jupyter_code>options(repr.plot.width=15, repr.plot.height=10)
barplot( pc@loadings[,1], main = "PC1 loading plot", las = 2, cex.names=0.6)<jupyter_output><empty_output><jupyter_text>Don't panic! It looks messy because we have 2448 (significant) features here in this dataset. Contribution of each feature to PC1 is ranging from -0.04 to 0.04. We can try to visualize features that contribution heaviest to PC1, say > 0.035.<jupyter_code>barplot( sort( pc@loadings[,1][which(abs(pc@loadings[,1]) > 0.035)], decreasing = FALSE), main = "PC1 loading plot", las = 2, cex.names=0.6)<jupyter_output><empty_output><jupyter_text>How many features here? And their names?<jupyter_code>length(  sort( pc@loadings[,1][which(abs(pc@loadings[,1]) > 0.035)], decreasing = FALSE) )
names( sort( pc@loadings[,1][which(abs(pc@loadings[,1]) > 0.035)], decreasing = FALSE) )<jupyter_output><empty_output><jupyter_text>Multivariate analysis, PLS-DA<jupyter_code>if (!requireNamespace("ropls", quietly = TRUE))
    BiocManager::install("ropls")

suppressPackageStartupMessages(library("ropls"))
classes <- as.character(d5[1, ])
my.plsda <- ropls::opls(t(norm_d1), classes)<jupyter_output>PLS-DA
12 samples x 3826 variables and 1 response
standard scaling of predictors and response(s)
      R2X(cum) R2Y(cum) Q2(cum)  RMSEE pre ort pR2Y  pQ2
Total    0.662    0.997   0.991 0.0281   1   0 0.05 0.05<jupyter_text>PLS-DA was performed with only one component. We can show the distribution of our sample on this P1.<jupyter_code>pls_da_score = my.plsda@scoreMN %>%
        as.data.frame(.) %>%
        rownames_to_column(., var = "name") %>%
        merge(rownames_to_column(as.data.frame(t(labels_d1)), var = "name"), by = "name")
pls_da_score
options(repr.plot.width=15, repr.plot.height=2)
ggplot(pls_da_score, aes(x = p1, y = 0, color = Label )) + 
        geom_point(size = 8) + theme(axis.title.y = element_blank(), text = element_text(size = 20), axis.text.y = element_blank(), axis.ticks.y = element_blank())
vip <- sort(ropls::getVipVn(my.plsda), decreasing = T) %>% as.data.frame(.) 
colnames(vip)="VIP"
head(vip)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12290828901284213, 'fraction_numerical': 0.046958100921001424, 'mean_word_length': 5.290901672786617, 'pattern_counts': {'":': 0, '<': 161, '<?xml version=': 0, '>': 132, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '190856', 'license': 'MIT', 'n_tokens_mistral': 5997, 'n_tokens_neox': 5313, 'n_words': 1497, 'path': 'Data_processing_and_statistics.ipynb', 'repo_name': 'whats-in-the-box/tutorials_and_demos', 'size': '15426'}