tests/testthat/test-protein_analysis.R

In promor: Proteomics Data Analysis and Modeling Tools

# Pre-processing----------------------------------------------------------------
## No tech-reps-----------------------------------------------------------------
# Create test_rawdf
test_raw_df <- structure(list(H_1 = c(
  27.3134575527802, 25.0475589619685, 25.7489360071628,
  22.3217409656277, 23.9379563661029, 25.768729719086, 27.6466382487042,
  NA, 26.2290168009068, 20.7899474944062, 20.6427110188685, NA,
  23.3404162709489, 23.9196169668074, 23.9141515162468
), H_2 = c(
  27.3965634965896,
  25.4638231809381, 25.7198105464424, 23.7216107357952, 23.7859632470773,
  25.1020545097735, 27.8039664255518, 21.3147336941717, 26.6504781237034,
  20.3212455062761, NA, NA, 23.5160500915308, 24.1067737781079,
  23.7869603409319
), H_3 = c(
  27.0682521470214, 24.5609251894038,
  25.8093128192732, 23.0995044870265, 23.4582634148661, 24.3433337543312,
  27.7661660273555, 20.9903006674839, 26.9776011840784, 21.7031583967283,
  NA, NA, 23.4146749142913, NA, 23.9672804455828
), L_4 = c(
  27.5967623420646,
  26.1093667962031, 26.3575984475757, 21.2096722968756, 24.2172666952823,
  24.7096194080498, 27.915816116319, 22.5738542523838, 27.3090394121738,
  23.0277123297888, 21.4292575700939, NA, 24.4428468344072, 25.0422241309514,
  24.9093979862761
), L_5 = c(
  27.8685009058466, 26.1012931825391,
  26.5974729242791, 23.7620267328653, 24.9532259575982, 25.0495544486124,
  27.9875559091121, 21.4550803438114, 27.2600639899901, 22.5863936035805,
  21.7253396278684, 23.3327894314844, 24.0624994391506, 25.141762588201,
  24.1083700330184
), L_6 = c(
  27.9516039906164, 25.7985259878485,
  26.452446499214, 22.6047857618104, 24.2072461353296, 24.6898717464561,
  27.8163496237955, 21.4407576690532, 27.1474412722278, 22.4058228851568,
  21.2955245730356, NA, 24.2565231471014, 24.651792303055, 24.2982404854359
)), row.names = c(
  "A0AV96;B7Z8Z7;A0AV96-2;D6R9D6", "A0AVT1;A0AVT1-2",
  "H7BXI1;A0FGR8-6;A0FGR8-2;A0FGR8;C9JGI7;A0FGR8-4", "A0JLT2;A0JLT2-2",
  "A0MZ66-4;A0MZ66;A0MZ66-3;B7Z7Z9;A0MZ66-5;A0MZ66-6;A0MZ66-2",
  "A1A528;O43264", "A1L0T0;E9PJS0", "A1L491;Q86XK2;Q86XK2-5;Q86XK2-2;H0YAV3;Q86XK2-3",
  "Q99798;A2A274;F5H2A5;B4DW08;B4DEC3", "Q86X10-2;A2A2F0;Q86X10-3;A2A2E9;Q86X10;B4E2E8;Q86X10-4",
  "Q9Y672;A2A2G4", "Q9Y312;A2A2Q9", "A2A2V2;P42696;P42696-2;Q5TCT4",
  "P35611-2;E7EV99;Q86XM2;A2A3N8;E7ENY0;P35611;P35611-3;H0Y9H2;Q96D30",
  "A2RRF3;Q9UBC2;Q9UBC2-2"
), class = "data.frame")

# Check if the function works
test_that("create_df works for data without tech reps", {
  expect_equal(
    create_df(
      "./testdata/test_protgroups.txt",
      "./testdata/test_expdesign.txt"
    ),
    test_raw_df
  )
})


# Filter out rows with high levels of NA in each group
# Create test_rawdf_filtered
test_raw_df_filt <- structure(list(H_1 = c(
  27.3134575527802, 25.0475589619685, 25.7489360071628,
  22.3217409656277, 23.9379563661029, 25.768729719086, 27.6466382487042,
  26.2290168009068, 20.7899474944062, 23.3404162709489, 23.9141515162468
), H_2 = c(
  27.3965634965896, 25.4638231809381, 25.7198105464424,
  23.7216107357952, 23.7859632470773, 25.1020545097735, 27.8039664255518,
  26.6504781237034, 20.3212455062761, 23.5160500915308, 23.7869603409319
), H_3 = c(
  27.0682521470214, 24.5609251894038, 25.8093128192732,
  23.0995044870265, 23.4582634148661, 24.3433337543312, 27.7661660273555,
  26.9776011840784, 21.7031583967283, 23.4146749142913, 23.9672804455828
), L_4 = c(
  27.5967623420646, 26.1093667962031, 26.3575984475757,
  21.2096722968756, 24.2172666952823, 24.7096194080498, 27.915816116319,
  27.3090394121738, 23.0277123297888, 24.4428468344072, 24.9093979862761
), L_5 = c(
  27.8685009058466, 26.1012931825391, 26.5974729242791,
  23.7620267328653, 24.9532259575982, 25.0495544486124, 27.9875559091121,
  27.2600639899901, 22.5863936035805, 24.0624994391506, 24.1083700330184
), L_6 = c(
  27.9516039906164, 25.7985259878485, 26.452446499214,
  22.6047857618104, 24.2072461353296, 24.6898717464561, 27.8163496237955,
  27.1474412722278, 22.4058228851568, 24.2565231471014, 24.2982404854359
)), row.names = c(
  "A0AV96;B7Z8Z7;A0AV96-2;D6R9D6", "A0AVT1;A0AVT1-2",
  "H7BXI1;A0FGR8-6;A0FGR8-2;A0FGR8;C9JGI7;A0FGR8-4", "A0JLT2;A0JLT2-2",
  "A0MZ66-4;A0MZ66;A0MZ66-3;B7Z7Z9;A0MZ66-5;A0MZ66-6;A0MZ66-2",
  "A1A528;O43264", "A1L0T0;E9PJS0", "Q99798;A2A274;F5H2A5;B4DW08;B4DEC3",
  "Q86X10-2;A2A2F0;Q86X10-3;A2A2E9;Q86X10;B4E2E8;Q86X10-4",
  "A2A2V2;P42696;P42696-2;Q5TCT4",
  "A2RRF3;Q9UBC2;Q9UBC2-2"
), class = "data.frame")

# Check if the function works
test_that("filterbygroup_na works ", {
  expect_equal(
    filterbygroup_na(test_raw_df, set_na = 0.33),
    test_raw_df_filt
  )
})


## Data with tech-reps----------------------------------------------------------
# Create test_rawdf_tr
test_raw_df_tr <- structure(list(WT_4_1 = c(
  28.2936438898012, 29.6186292576641,
  NA, 30.7945158173451, 29.0678932740327, 28.0406156470249, 29.1721214197861,
  27.119910145335, 30.6453846495587, 30.0410074885123, 28.5316673138427
), WT_4_2 = c(
  28.2302727014579, 29.732064069791, 23.9675443129893,
  30.8093894910794, 29.1898754836933, 27.7491918268356, 29.0156859984869,
  27.2171557372336, 30.7830833879877, 30.109610898189, 28.1675349098489
), WT_4_3 = c(
  28.2168783047708, 29.7297001102491, NA, 30.8191699938329,
  29.1923705191376, 28.0976311997323, 28.8926709332398, 27.2714185722088,
  30.8046281169856, 29.9632008573523, 28.6045595631519
), WT_5_1 = c(
  28.0692770067938,
  29.4347383877629, NA, 30.5958380616431, 28.9967423270321, 27.914619166411,
  29.4056800832666, 26.9129646012963, 30.9832533443557, 29.3775145533344,
  28.3533388778899
), WT_5_2 = c(
  28.3002932737853, 29.6902082063487,
  NA, 30.7136658360682, 28.8854739844503, 27.6736421690804, 28.9897497517868,
  27.2529555521063, 31.0995511375747, 29.3940094928478, 28.1271622395638
), WT_5_3 = c(
  28.1671520416985, 29.7997429229766, NA, 30.747912274609,
  28.9950175821295, 27.6961748124338, 28.5973308358325, 27.2422722265493,
  30.9876643253623, 29.4029389567881, 28.3268456315426
), WT_6_1 = c(
  28.3297564494538,
  29.8846446110208, NA, 30.5892442985113, 29.1034458987105, 28.0054423458266,
  28.9629803081492, 27.3427405583237, 30.9307229924469, 29.4434002727089,
  28.4333659551936
), WT_6_2 = c(
  28.3335999942551, 29.9338210210556,
  NA, 30.7303643531446, 29.1101836284433, 27.9655971554024, 29.1076918404721,
  27.2578182701616, 30.9116367246489, 29.5180706258109, 28.1363842174365
), WT_6_3 = c(
  28.309776709179, 29.9616839002322, 24.3187242869872,
  30.6472732188454, 28.9547130090708, 27.9666425667018, 28.9621805345338,
  27.3910001879615, 30.8879447265738, 29.5763081393039, 28.2454039333088
), D8_1_1 = c(
  27.8445776263161, 29.448430767204, NA, 30.6236609447517,
  29.1416310525644, 27.7035201554364, 29.2664228386558, 27.2927228106225,
  31.0073151073851, 30.2687444379105, 28.2735543165528
), D8_1_2 = c(
  28.1898048066551,
  29.8917877583825, NA, 30.6470158322001, 28.9819070074915, 27.7271537144422,
  29.4102977482881, 27.5540597213853, 31.0176384315078, 30.3827796811565,
  28.2186806628685
), D8_1_3 = c(
  28.2679412301097, 29.7679187969981,
  NA, 30.718573244987, 29.0732414760752, 27.7022654943766, 29.1479264237409,
  27.5512016651033, 31.0933932835392, 30.3665094859859, 28.4562278107371
), D8_2_1 = c(
  28.2224629079324, 29.8526210142944, NA, 30.5359364692179,
  28.940669358688, 27.8083924396305, 29.1452317705679, 27.4432421691426,
  31.2667355793252, 29.237994874927, 27.9604139464238
), D8_2_2 = c(
  28.2595253359549,
  29.9418169204995, NA, 30.52738548788, 28.9021634262318, 27.5600505637081,
  29.2013030632228, 27.414658786775, 31.2676845555708, 29.2900879108747,
  28.2321956846597
), D8_2_3 = c(
  28.3436302608953, 29.9901015395258,
  NA, 30.5929013221194, 28.8589613663302, 27.5899228860983, 29.2357829739292,
  27.4522265330772, 31.3103511059622, 29.144405373359, 28.0186070260911
), D8_3_1 = c(
  28.0400407535959, 30.20559721838, NA, 30.6619537896456,
  29.0941864816023, 27.6380991218623, 29.4514603177377, 27.5613623034931,
  31.0582085194794, 29.3470357135931, 28.1787840509819
), D8_3_2 = c(
  28.3079544669147,
  30.0972293208228, NA, 30.6279276018739, 28.8933076419283, 27.5586646543679,
  29.3034323171483, 27.5390468442083, 30.8756217150601, 29.3415110187243,
  28.0666621741898
), D8_3_3 = c(
  28.286390947795, 30.218241834318,
  NA, 30.6085834258799, 29.0433563051252, 27.4601398881949, 29.2906811830317,
  27.4296598378651, 30.9773964738913, 29.3815012993247, 28.1019443425373
)), row.names = c(
  "gi|118496616", "gi|118496617", "gi|118496621",
  "gi|118496635", "gi|118496636", "gi|118496637", "gi|118496639",
  "gi|118496647", "gi|118496648", "gi|118496649", "gi|118496650"
), class = "data.frame")

# Check if the function works
test_that("create_df works for data with tech reps", {
  expect_equal(
    create_df("./testdata/test_protgroups_tr.txt",
      "./testdata/test_expdesign_tr.txt",
      tech_reps = TRUE
    ),
    test_raw_df_tr
  )
})

# Average across tech reps
# Create test_rawdf_avg_tr

test_raw_df_tr_avg <- structure(list(WT_4 = c(
  28.24693163201, 29.6934644792347, 23.9675443129893,
  30.8076917674191, 29.1500464256212, 27.9624795578643, 29.0268261171709,
  27.2028281515925, 30.744365384844, 30.0379397480179, 28.4345872622812
), WT_5 = c(
  28.1789074407592, 29.6415631723627, NaN, 30.6858053907734,
  28.9590779645373, 27.7614787159751, 28.9975868902953, 27.1360641266506,
  31.0234896024309, 29.3914876676568, 28.2691155829988
), WT_6 = c(
  28.3243777176293,
  29.9267165107695, 24.3187242869872, 30.6556272901671, 29.0561141787415,
  27.9792273559769, 29.010950894385, 27.3305196721489, 30.9101014812232,
  29.5125930126079, 28.271718035313
), D8_1 = c(
  28.1007745543603,
  29.7027124408615, NaN, 30.6630833406463, 29.0655931787104, 27.7109797880851,
  29.2748823368949, 27.4659947323704, 31.0394489408107, 30.3393445350176,
  28.3161542633861
), D8_2 = c(
  28.2752061682609, 29.9281798247732,
  NaN, 30.5520744264058, 28.9005980504167, 27.6527886298123, 29.1941059359066,
  27.4367091629983, 31.2815904136194, 29.2241627197202, 28.0704055523916
), D8_3 = c(
  28.2114620561019, 30.1736894578403, NaN, 30.6328216057998,
  29.0102834762186, 27.552301221475, 29.3485246059726, 27.5100229951888,
  30.9704089028103, 29.3566826772141, 28.115796855903
)), class = "data.frame", row.names = c(
  "gi|118496616",
  "gi|118496617", "gi|118496621", "gi|118496635", "gi|118496636",
  "gi|118496637", "gi|118496639", "gi|118496647", "gi|118496648",
  "gi|118496649", "gi|118496650"
))

# Check if the function works
test_that("aver_techreps works", {
  expect_equal(
    aver_techreps(test_raw_df_tr),
    test_raw_df_tr_avg
  )
})

## Removing samples-------------------------------------------------------------
# Make a test data frame
df <- data.frame(
  a = c(12, 34, 11, 83, NA, 67),
  b = c(NA, NA, 39, 29, 20, NA),
  c = c(82, 40, 22, NA, NA, 70)
)
# Create df_rem
df_rem <- structure(list(a = c(12, 34, 11, 83, NA, 67), c = c(
  82, 40, 22,
  NA, NA, 70
)), class = "data.frame", row.names = c(NA, -6L))


# Check if the function works
test_that("rem_sample works", {
  expect_equal(
    rem_sample(df, "b"),
    df_rem
  )
})

# Missing data imputation-------------------------------------------------------

# minprob results
df_mp <- structure(
  list(
    a = c(12, 34, 11, 83, 9.80950652319248, 67),
    b = c(12.6260927833319, 29.8391137373841, 39, 29, 20, -1.37591651207944),
    c = c(82, 40, 22, 14.7098046749082, 22.299106363566, 70)
  ),
  row.names = c(NA, -6L),
  class = "data.frame"
)

# rf - results
df_rf <- structure(
  list(
    a = c(12, 34, 11, 83, 62.1633333333333, 67),
    b = c(26.6, 32.15, 39, 29, 20, 22.3),
    c = c(82, 40, 22, 53.125, 66.025, 70)
  ),
  row.names = c(NA, -6L),
  class = "data.frame"
)

# knn - results
df_knn <- structure(
  list(
    a = c(12, 34, 11, 83, 34, 67),
    b = c(39, 39, 39, 29, 20, 39),
    c = c(82, 40, 22, 70, 70, 70)
  ),
  row.names = c("1", "2", "3", "4", "5", "6"),
  class = "data.frame"
)

# svd - results
df_svd <- structure(
  list(
    a = c(12, 34, 11, 83, 40.5316879411941, 67),
    b = c(15.8432402677882, 34.3000432486029, 39, 29, 20, 24.3678401876791),
    c = c(82, 40, 22, 59.6698946596492, 73.7190158521945, 70)
  ),
  class = "data.frame",
  row.names = c(NA, -6L)
)

# mindet - results
df_md <- structure(
  list(
    a = c(12, 34, 11, 83, 11.04, 67),
    b = c(20.18, 20.18, 39, 29, 20, 20.18),
    c = c(82, 40, 22, 22.54, 22.54, 70)
  ),
  row.names = c(NA, -6L),
  class = "data.frame"
)

# testing minprob
test_that("impute_na works for minprob", {
  expect_equal(
    impute_na(df,
      seed = 327
    ),
    df_mp
  )
})

# testing rf
test_that("impute_na works for rf", {
  suppressWarnings(expect_equal(impute_na(df,
    method = "RF",
    seed = 327
  ), df_rf))
})

# testing knn
test_that("impute_na works for knn", {
  expect_equal(impute_na(df,
    method = "kNN",
    seed = 327
  ), df_knn)
})

# testing svd
test_that("impute_na works for svd", {
  suppressMessages(expect_equal(impute_na(df,
    method = "SVD",
    seed = 327
  ), df_svd))
})

# testing minDet
test_that("impute_na works for minDet", {
  expect_equal(impute_na(df,
    method = "minDet",
    seed = 327
  ), df_md)
})

# Normalization------------------------------------------------------------------
# quantile - results
df_norm_q <- structure(
  c(
    18.0997021211887, 34.3333333333333, 15.208697594444,
    68, 7.71446489534041, 55.6130379124614, 15.208697594444,
    55.6130379124614,
    68, 34.3333333333333, 18.0997021211887, 7.71446489534041, 68,
    34.3333333333333, 15.208697594444, 7.71446489534041,
    18.0997021211887,
    55.6130379124614
  ),
  dim = c(6L, 3L),
  dimnames = list(NULL, c("a", "b", "c"))
)

# scale -results
df_norm_s <- structure(
  c(
    13.5592392536598, 38.4178445520362, 12.4293026491882,
    93.7847381711472, 11.0841204923586, 75.7057524996008,
    13.3932139334231,
    31.6520431717492, 41.3695156820846, 30.7619475584732,
    21.2151362472229,
    -1.45951281342845, 68.4138364231764, 33.3726031332568,
    18.3549317232912,
    12.2726118395859, 18.6044806724393, 58.4020554831994
  ),
  dim = c(6L, 3L),
  dimnames = list(NULL, c("a", "b", "c"))
)

# cyclicloess - results
df_norm_c <- structure(
  c(
    8.87594722638867, 35.9180013309211, 24.1612907491603,
    43.5415649767626, 17.3396297815539, 44.825964199607,
    21.9075705937958,
    38.5487769581834, 24.1612907491603, 43.5415649767626,
    17.3396297815539,
    44.825964199607, 81.9288096106512, 35.9180013309212,
    24.1612907491603,
    43.5415649767626, 17.3396297815539, 44.825964199607
  ),
  dim = c(6L, 3L),
  dimnames = list(NULL, c("a", "b", "c"))
)


## testing - quantile method----------------------------------------------------

test_that("normalization works for quantile", {
  expect_equal(
    normalize_data(df_mp),
    df_norm_q
  )
})

## testing - scale method-------------------------------------------------------

test_that("normalization works for sacle", {
  expect_equal(
    normalize_data(df_mp,
      method = "scale"
    ),
    df_norm_s
  )
})

## testing- cyclicloess method--------------------------------------------------

test_that("normalization works for cyclicloess", {
  expect_equal(
    normalize_data(df_mp,
      method = "cyclicloess"
    ),
    df_norm_c
  )
})

# Differential expression -------------------------------------------------------
# make a fake - norm_df
set.seed(581)
sigma2 <- 0.05 / rchisq(1000, df = 10) * 10
fake_norm_df <- matrix(rnorm(100 * 6, sd = sqrt(sigma2)), 100, 6)
fake_norm_df[1, 4:6] <- fake_norm_df[1, 4:6] + 1
colnames(fake_norm_df) <- c("Ctrl_1", "Ctrl_2", "Ctrl_3", "Case_1", "Case_2", "Case_3")

# Load fit_df
fake_fit_df <- readRDS("./testdata/fake_fit_df.rds")

# Test find_dep works
test_that("find_dep works", {
  expect_equal(
    find_dep(fake_norm_df),
    fake_fit_df
  )
})