Interpreting an ML model with SHAP values
Mark Ziemann
2026-07-14
Source: https://github.com/markziemann/bioinformatics_intro_workshop
Introduction
The Shapley object contains the contribution of each feature to the prediction for the specific observation.
each feature has a corresponding SHAP value, indicating its contribution to the model’s prediction for that observation. Positive values indicate that the feature increased the prediction, while negative values suggest it decreased the prediction.
In this workshop, we will start working on some RNA-seq data. Here is the citation:
Cho JW, Shim HS, Lee CY, Park SY et al. The importance of enhancer methylation for epigenetic regulation of tumorigenesis in squamous lung cancer. Exp Mol Med 2022 Jan;54(1):12-22. PMID: 34987166
suppressPackageStartupMessages({
library("kableExtra")
library("DESeq2")
library("gplots")
library("beeswarm")
library("parallel")
library("kableExtra")
library("xgboost")
library("iml")
library("ggplot2")
library("beeswarm")
})Download gene expression counts and read the file
Download from NCBI.
DATA="GSE158420.counts.txt.gz"
if ( ! file.exists(DATA) ) {
download.file("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE158420&format=file&file=GSE158420%5Fcounts%2Etxt%2Egz",
destfile=DATA)
}Read data. Note that there are a couple of duplicated gene symbols. Here I use aggregate() command to deal with them.
mx <- read.table(DATA,row.names = NULL)
str(mx)## 'data.frame': 18747 obs. of 75 variables:
## $ row.names: chr "OR4F5" "OR4F29" "OR4F16" "LINC00115" ...
## $ X10N : int 6 0 0 14 15 1727 43 8 53 98 ...
## $ X10T : int 30 6 10 54 115 1232 190 227 215 107 ...
## $ X11N : int 10 0 0 66 59 2342 55 57 70 57 ...
## $ X11T : int 2 0 0 27 18 3205 74 120 59 68 ...
## $ X12N : int 1 0 0 22 15 1563 21 6 28 64 ...
## $ X12T : int 2 0 0 31 12 1393 36 20 24 29 ...
## $ X14N : int 2 0 0 41 7 1651 49 13 34 43 ...
## $ X14T : int 0 0 0 34 18 2373 55 45 87 282 ...
## $ X15N : int 2 0 0 37 3 1418 41 7 35 34 ...
## $ X15T : int 0 0 0 18 18 1004 33 11 17 52 ...
## $ X16N : int 0 0 0 30 12 1939 49 3 69 106 ...
## $ X16T : int 0 0 0 44 16 2216 63 25 38 37 ...
## $ X17N : int 1 0 0 31 13 1653 39 6 63 76 ...
## $ X17T : int 0 0 0 15 12 2787 18 11 9 54 ...
## $ X18N : int 0 0 0 18 9 1540 46 9 80 99 ...
## $ X18T : int 0 0 5 89 271 994 420 351 64 95 ...
## $ X19N : int 0 0 0 22 14 1510 29 8 28 30 ...
## $ X19T : int 13 2 0 55 72 2051 149 138 125 70 ...
## $ X20N : int 0 0 0 25 16 1656 44 17 83 51 ...
## $ X20T : int 4 0 0 43 30 1723 60 21 62 29 ...
## $ X21N : int 0 0 0 39 12 1824 47 7 139 705 ...
## $ X21T : int 0 0 0 22 30 3537 47 25 34 112 ...
## $ X22N : int 0 0 0 25 9 1590 54 15 40 97 ...
## $ X22T : int 9 0 0 79 90 5094 90 11 101 95 ...
## $ X23N : int 0 0 0 28 25 1499 44 12 79 52 ...
## $ X23T : int 0 0 0 31 57 1723 64 46 56 25 ...
## $ X24N : int 0 0 0 32 2 1613 44 5 54 86 ...
## $ X24T : int 0 0 0 31 15 4485 124 78 249 90 ...
## $ X25N : int 0 0 0 34 19 1533 47 13 60 39 ...
## $ X25T : int 0 0 0 30 11 2369 59 35 58 49 ...
## $ X26N : int 1 0 0 31 16 1407 40 9 44 57 ...
## $ X26T : int 6 1 3 31 26 1772 45 36 61 23 ...
## $ X27N : int 0 0 0 39 13 2051 74 15 28 48 ...
## $ X27T : int 3 0 0 29 29 1808 52 64 81 110 ...
## $ X28N : int 0 0 0 13 15 1824 73 11 71 74 ...
## $ X28T : int 0 0 0 8 1 3714 35 49 31 60 ...
## $ X29N : int 0 0 0 18 28 1591 39 22 62 19 ...
## $ X29T : int 0 0 0 29 19 3730 90 55 47 92 ...
## $ X2N : int 0 0 0 32 16 1739 54 8 45 118 ...
## $ X2T : int 0 0 0 19 15 3025 93 64 100 108 ...
## $ X30N : int 0 0 0 23 41 2500 51 21 98 106 ...
## $ X30T : int 0 0 1 48 119 5999 231 41 79 44 ...
## $ X31N : int 0 0 0 35 14 1729 66 15 51 28 ...
## $ X31T : int 0 0 0 24 26 4151 105 103 130 67 ...
## $ X32N : int 0 0 0 27 16 1666 59 19 48 45 ...
## $ X32T : int 0 0 0 16 4 1195 50 17 10 57 ...
## $ X33N : int 0 0 0 46 16 1731 122 23 52 52 ...
## $ X33T : int 0 0 0 36 50 2503 49 28 67 48 ...
## $ X34N : int 0 0 0 29 16 1745 97 12 60 62 ...
## $ X34T : int 0 0 0 39 92 1673 56 3 91 64 ...
## $ X35N : int 0 0 0 32 19 1678 31 13 114 689 ...
## $ X35T : int 1 0 0 27 43 1709 91 67 134 698 ...
## $ X36N : int 0 0 0 17 10 1885 26 14 95 63 ...
## $ X36T : int 1 0 0 24 48 1105 142 72 114 269 ...
## $ X37N : int 0 0 0 28 7 2478 39 19 59 59 ...
## $ X37T : int 0 0 0 7 25 1630 32 9 53 52 ...
## $ X38N : int 0 0 0 12 15 1576 24 4 70 110 ...
## $ X38T : int 0 0 0 27 85 2113 137 117 107 70 ...
## $ X39N : int 1 0 0 10 11 1519 38 12 15 42 ...
## $ X39T : int 0 0 0 29 40 2100 80 27 47 33 ...
## $ X3N : int 0 0 0 70 6 1788 53 10 48 47 ...
## $ X3T : int 0 2 0 53 56 3764 149 188 117 48 ...
## $ X4N : int 1 0 0 18 35 2057 47 44 68 44 ...
## $ X4T : int 0 0 0 23 28 2110 58 48 42 19 ...
## $ X5N : int 0 0 0 53 13 1273 43 2 34 47 ...
## $ X5T : int 0 0 0 41 12 2708 62 114 101 382 ...
## $ X6N : int 1 0 0 14 11 1478 26 6 48 48 ...
## $ X6T : int 0 0 0 29 4 2454 47 52 97 378 ...
## $ X7N : int 0 0 0 44 10 1722 49 13 77 75 ...
## $ X7T : int 0 0 0 28 21 1892 54 49 64 240 ...
## $ X8N : int 23 0 0 53 33 1442 110 51 53 38 ...
## $ X8T : int 10 0 0 46 28 1931 72 75 64 83 ...
## $ X9N : int 0 0 0 11 16 1701 23 16 43 55 ...
## $ X9T : int 0 0 0 15 5 2018 25 13 92 15 ...which(duplicated(mx$row.names))## [1] 4765 15632xa <- aggregate(. ~ row.names,mx,sum)
rownames(xa) <- xa$row.names
str(xa)## 'data.frame': 18745 obs. of 75 variables:
## $ row.names: chr "1-Dec" "1-Mar" "1-Sep" "10-Mar" ...
## $ X10N : int 0 431 160 286 3541 0 3164 1 0 1075 ...
## $ X10T : int 31 461 131 217 4311 40 2874 72 94 440 ...
## $ X11N : int 0 478 280 130 2174 6 4758 10 4 816 ...
## $ X11T : int 1 98 131 58 3675 1 1750 12 1 204 ...
## $ X12N : int 0 862 132 2 2856 0 5592 1 0 729 ...
## $ X12T : int 1 635 383 18 1864 0 4185 3 0 619 ...
## $ X14N : int 0 489 381 19 3306 0 3610 0 0 988 ...
## $ X14T : int 0 499 124 39 1546 0 2181 10 1 589 ...
## $ X15N : int 0 543 114 417 3507 0 5148 0 1 828 ...
## $ X15T : int 0 608 217 68 2483 1 3100 1 0 818 ...
## $ X16N : int 0 707 258 12 3568 0 4271 0 0 1339 ...
## $ X16T : int 0 498 886 90 1294 0 2933 0 3 882 ...
## $ X17N : int 0 700 302 9 3658 0 3515 4 1 1192 ...
## $ X17T : int 1 411 162 122 542 2 3093 13 2 317 ...
## $ X18N : int 0 607 133 123 3337 1 3549 0 0 1108 ...
## $ X18T : int 79 446 279 313 1050 38 1356 171 127 554 ...
## $ X19N : int 0 608 138 70 4527 0 3569 0 0 1290 ...
## $ X19T : int 18 463 88 127 1865 18 2033 61 97 318 ...
## $ X20N : int 0 361 239 1481 2656 0 3061 7 2 783 ...
## $ X20T : int 1 488 533 1010 1314 0 2811 9 2 759 ...
## $ X21N : int 0 615 212 97 4449 1 4153 0 3 1301 ...
## $ X21T : int 1 350 53 7 1060 3 1586 2 358 182 ...
## $ X22N : int 0 673 258 94 3196 0 3975 0 1 998 ...
## $ X22T : int 0 1838 126 21 400 3 1693 10 3 414 ...
## $ X23N : int 0 297 197 588 3499 1 4309 3 0 871 ...
## $ X23T : int 1 262 229 694 1709 2 3004 1 1 551 ...
## $ X24N : int 0 762 137 92 3280 0 3403 0 0 927 ...
## $ X24T : int 2 529 236 24 3528 5 2128 36 1 498 ...
## $ X25N : int 1 534 132 87 3463 4 2852 0 0 802 ...
## $ X25T : int 0 860 66 23 2832 3 4208 2 4 375 ...
## $ X26N : int 0 451 98 126 2893 0 4659 3 0 727 ...
## $ X26T : int 3 1110 56 368 3445 3 1820 9 1 351 ...
## $ X27N : int 0 798 573 192 2442 0 3134 0 1 1092 ...
## $ X27T : int 2 380 190 114 3703 4 2245 7 0 480 ...
## $ X28N : int 0 572 211 441 3250 0 5131 1 2 1150 ...
## $ X28T : int 2 336 75 5 5465 0 2009 86 2 156 ...
## $ X29N : int 0 289 79 2242 2028 2 3068 8 3 414 ...
## $ X29T : int 0 514 205 19 1971 1 2151 1 1 389 ...
## $ X2N : int 0 575 151 23 3502 0 4091 0 3 1056 ...
## $ X2T : int 0 257 147 80 2354 0 3502 0 0 404 ...
## $ X30N : int 0 464 183 616 2605 1 4596 1 0 836 ...
## $ X30T : int 0 483 144 24 934 3 3291 5 1 111 ...
## $ X31N : int 0 535 212 589 3291 1 3648 0 1 913 ...
## $ X31T : int 2 478 144 28 4226 5 2889 6 8 200 ...
## $ X32N : int 0 311 302 48 3135 0 1931 5 0 950 ...
## $ X32T : int 0 280 1086 27 1416 0 1278 2 0 155 ...
## $ X33N : int 0 630 231 4 3072 0 5982 1 0 1038 ...
## $ X33T : int 0 437 363 11 1702 1 2752 6 0 752 ...
## $ X34N : int 1 713 185 956 3356 1 3695 1 0 1495 ...
## $ X34T : int 0 334 290 13 2515 0 2415 1 0 708 ...
## $ X35N : int 0 643 125 640 3034 0 4439 2 0 1104 ...
## $ X35T : int 1 789 199 63 2532 1 3526 13 1 948 ...
## $ X36N : int 0 286 61 654 3678 1 2726 1 1 1142 ...
## $ X36T : int 1 339 355 19 502 6 1806 12 4 612 ...
## $ X37N : int 0 488 59 1234 1512 0 2326 0 0 419 ...
## $ X37T : int 0 320 121 2076 1324 6 2052 0 0 413 ...
## $ X38N : int 0 427 92 380 2684 0 3143 0 0 829 ...
## $ X38T : int 4 507 203 185 2542 9 2961 22 12 841 ...
## $ X39N : int 0 695 329 43 1770 1 2849 1 2 671 ...
## $ X39T : int 1 458 174 135 1446 2 2355 10 1 1100 ...
## $ X3N : int 1 628 120 11 2921 0 5548 0 1 797 ...
## $ X3T : int 3 287 114 63 1688 10 2582 25 9 257 ...
## $ X4N : int 0 305 119 323 2692 2 3613 8 3 762 ...
## $ X4T : int 2 255 59 33 2447 3 1558 17 15 344 ...
## $ X5N : int 0 801 144 78 3367 0 5042 1 0 1067 ...
## $ X5T : int 8 364 44 63 4396 1 1563 8 4 88 ...
## $ X6N : int 2 626 110 28 3830 0 4163 0 1 1024 ...
## $ X6T : int 0 320 184 6 2841 0 2520 2 2 216 ...
## $ X7N : int 1 638 181 323 2623 2 3087 0 1 1052 ...
## $ X7T : int 1 412 327 16 3671 1 2413 7 4 629 ...
## $ X8N : int 17 264 101 185 1614 7 2039 37 44 376 ...
## $ X8T : int 2 436 217 168 2102 9 1930 20 11 307 ...
## $ X9N : int 0 504 191 76 3507 0 3858 3 5 973 ...
## $ X9T : int 2 317 48 7 2482 5 2377 0 1 232 ...xa$row.names = NULL
head(xa)## X10N X10T X11N X11T X12N X12T X14N X14T X15N X15T X16N X16T X17N X17T
## 1-Dec 0 31 0 1 0 1 0 0 0 0 0 0 0 1
## 1-Mar 431 461 478 98 862 635 489 499 543 608 707 498 700 411
## 1-Sep 160 131 280 131 132 383 381 124 114 217 258 886 302 162
## 10-Mar 286 217 130 58 2 18 19 39 417 68 12 90 9 122
## 10-Sep 3541 4311 2174 3675 2856 1864 3306 1546 3507 2483 3568 1294 3658 542
## 11-Mar 0 40 6 1 0 0 0 0 0 1 0 0 0 2
## X18N X18T X19N X19T X20N X20T X21N X21T X22N X22T X23N X23T X24N X24T
## 1-Dec 0 79 0 18 0 1 0 1 0 0 0 1 0 2
## 1-Mar 607 446 608 463 361 488 615 350 673 1838 297 262 762 529
## 1-Sep 133 279 138 88 239 533 212 53 258 126 197 229 137 236
## 10-Mar 123 313 70 127 1481 1010 97 7 94 21 588 694 92 24
## 10-Sep 3337 1050 4527 1865 2656 1314 4449 1060 3196 400 3499 1709 3280 3528
## 11-Mar 1 38 0 18 0 0 1 3 0 3 1 2 0 5
## X25N X25T X26N X26T X27N X27T X28N X28T X29N X29T X2N X2T X30N X30T
## 1-Dec 1 0 0 3 0 2 0 2 0 0 0 0 0 0
## 1-Mar 534 860 451 1110 798 380 572 336 289 514 575 257 464 483
## 1-Sep 132 66 98 56 573 190 211 75 79 205 151 147 183 144
## 10-Mar 87 23 126 368 192 114 441 5 2242 19 23 80 616 24
## 10-Sep 3463 2832 2893 3445 2442 3703 3250 5465 2028 1971 3502 2354 2605 934
## 11-Mar 4 3 0 3 0 4 0 0 2 1 0 0 1 3
## X31N X31T X32N X32T X33N X33T X34N X34T X35N X35T X36N X36T X37N X37T
## 1-Dec 0 2 0 0 0 0 1 0 0 1 0 1 0 0
## 1-Mar 535 478 311 280 630 437 713 334 643 789 286 339 488 320
## 1-Sep 212 144 302 1086 231 363 185 290 125 199 61 355 59 121
## 10-Mar 589 28 48 27 4 11 956 13 640 63 654 19 1234 2076
## 10-Sep 3291 4226 3135 1416 3072 1702 3356 2515 3034 2532 3678 502 1512 1324
## 11-Mar 1 5 0 0 0 1 1 0 0 1 1 6 0 6
## X38N X38T X39N X39T X3N X3T X4N X4T X5N X5T X6N X6T X7N X7T
## 1-Dec 0 4 0 1 1 3 0 2 0 8 2 0 1 1
## 1-Mar 427 507 695 458 628 287 305 255 801 364 626 320 638 412
## 1-Sep 92 203 329 174 120 114 119 59 144 44 110 184 181 327
## 10-Mar 380 185 43 135 11 63 323 33 78 63 28 6 323 16
## 10-Sep 2684 2542 1770 1446 2921 1688 2692 2447 3367 4396 3830 2841 2623 3671
## 11-Mar 0 9 1 2 0 10 2 3 0 1 0 0 2 1
## X8N X8T X9N X9T
## 1-Dec 17 2 0 2
## 1-Mar 264 436 504 317
## 1-Sep 101 217 191 48
## 10-Mar 185 168 76 7
## 10-Sep 1614 2102 3507 2482
## 11-Mar 7 9 0 5Basic QC - number of reads per sample
Less than 10 million is not acceptable. More than 20 million is considered good.
mycolsums <- colSums(xa)
mycolsums## X10N X10T X11N X11T X12N X12T X14N X14T
## 25014713 21358077 23719271 23009875 23540231 23439861 22829532 21990640
## X15N X15T X16N X16T X17N X17T X18N X18T
## 25092201 22135465 27933368 24850610 26095454 24379916 24331203 15718196
## X19N X19T X20N X20T X21N X21T X22N X22T
## 24878899 21283611 21487092 22777921 27374413 19265745 24658124 25809845
## X23N X23T X24N X24T X25N X25T X26N X26T
## 22491399 23508958 24018825 31021885 22675122 25197199 21507129 20386730
## X27N X27T X28N X28T X29N X29T X2N X2T
## 29175712 27733951 27749273 26197405 23156956 25636790 24559488 24150017
## X30N X30T X31N X31T X32N X32T X33N X33T
## 30343041 24469122 26249413 26005385 23779411 19207983 25951651 21803146
## X34N X34T X35N X35T X36N X36T X37N X37T
## 27848853 27618533 29607974 26996895 25443403 21856651 27091529 24555047
## X38N X38T X39N X39T X3N X3T X4N X4T
## 21897996 28799121 21262314 25770806 26104237 24647793 23048255 25691004
## X5N X5T X6N X6T X7N X7T X8N X8T
## 24139646 25661963 24030603 24994714 23238440 24211560 18291010 22430032
## X9N X9T
## 22689392 21734211summary(mycolsums)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15718196 22678690 24355560 24319084 25916200 31021885mycolsums <- sort(mycolsums)
head(mycolsums,10)## X18T X8N X32T X21T X26T X39N X19T X10T
## 15718196 18291010 19207983 19265745 20386730 21262314 21283611 21358077
## X20N X26N
## 21487092 21507129tail(mycolsums,10)## X34T X27T X28N X34N X16N X38T X27N X35N
## 27618533 27733951 27749273 27848853 27933368 28799121 29175712 29607974
## X30N X24T
## 30343041 31021885barplot(mycolsums)
barplot(head(mycolsums))
Does this dataset need filtering of lowly expressed genes?
There are ~60k genes in the main annotation sets. We can count the number of rows to see how many genes included in this dataset.
I also run my own threshold, which is average of 10 reads per sample.
dim(xa)## [1] 18745 74hist(rowMeans(xa))
table(rowMeans(xa)>=10)##
## FALSE TRUE
## 2309 16436xf <- xa[which(rowMeans(xa)>=10),]
dim(xa)## [1] 18745 74dim(xf)## [1] 16436 74Principal component analysis
Use PCA to understand the overall variation between samples in the dataset.
mds <- cmdscale(dist(t(scale(xf))),k=10)
barplot(colSums(abs(mds)),names.arg = 1:10)
head(mds)## [,1] [,2] [,3] [,4] [,5] [,6]
## X10N 37.643694 -13.024427 0.8597808 4.097087 6.112696 10.7146140
## X10T -76.198299 -19.565410 -12.7071833 -4.364407 -15.251200 -0.5565757
## X11N 30.689243 -1.139396 -1.9729685 -6.541190 -4.298140 -2.8939500
## X11T -38.586767 -1.714293 34.3064374 -4.932127 16.569860 17.1527763
## X12N 38.975575 -6.187959 -10.0227906 -12.146277 -12.417868 -10.0294316
## X12T 3.294181 45.630846 -38.5169987 -8.184535 -29.904639 12.8409121
## [,7] [,8] [,9] [,10]
## X10N 7.568131 -3.935890 5.9196796 8.1464807
## X10T -24.416829 -14.702973 10.6437095 6.0036817
## X11N -1.875214 3.404213 0.3012157 -0.1634874
## X11T -4.899219 2.537063 -4.3955324 -10.1840018
## X12N -3.669772 1.359340 1.6358922 1.5074445
## X12T -3.015237 2.774796 -1.9950542 5.0620557plot(mds[,1:2],pch=19,col="gray",
main="Dim 1 and 2",xlab="PCA1",ylab="PCA2")
text(mds[,1],mds[,2],labels = colnames(xf))
plot(mds[,c(1,3)],pch=19,col="gray",
main="Dim 1 and 3",xlab="PCA1",ylab="PCA3")
text(mds[,1],mds[,3],labels = colnames(xf))
plot(mds[,2:3],pch=19,col="gray",
main="Dim 2 and 3",xlab="PCA2",ylab="PCA3")
text(mds[,2],mds[,3],labels = colnames(xf))
Curate samplesheet
The samplesheet is the sample level metadata, which we need for the differential expression step.
ss <- as.data.frame(colnames(xa))
head(ss)## colnames(xa)
## 1 X10N
## 2 X10T
## 3 X11N
## 4 X11T
## 5 X12N
## 6 X12Tss$tumor <- as.numeric(grepl("T",ss[,1]))
ss$patient <- gsub("T","" , gsub("N","",ss[,1]))
rownames(ss) <- ss[,1]
head(ss)## colnames(xa) tumor patient
## X10N X10N 0 X10
## X10T X10T 1 X10
## X11N X11N 0 X11
## X11T X11T 1 X11
## X12N X12N 0 X12
## X12T X12T 1 X12DE
We will use DESeq2 for differential expression.
dds <- DESeqDataSetFromMatrix(countData = xf , colData=ss, design = ~ patient + tumor)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() functionres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
head(dge) |> kbl(caption="Top significant genes") |> kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | X10N | X10T | X11N | X11T | X12N | X12T | X14N | X14T | X15N | X15T | X16N | X16T | X17N | X17T | X18N | X18T | X19N | X19T | X20N | X20T | X21N | X21T | X22N | X22T | X23N | X23T | X24N | X24T | X25N | X25T | X26N | X26T | X27N | X27T | X28N | X28T | X29N | X29T | X2N | X2T | X30N | X30T | X31N | X31T | X32N | X32T | X33N | X33T | X34N | X34T | X35N | X35T | X36N | X36T | X37N | X37T | X38N | X38T | X39N | X39T | X3N | X3T | X4N | X4T | X5N | X5T | X6N | X6T | X7N | X7T | X8N | X8T | X9N | X9T | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| PDK1 | 937.9546 | 2.135424 | 0.1258421 | 16.96908 | 0 | 0 | 8.650724 | 10.642692 | 8.967844 | 10.498565 | 8.134786 | 10.591653 | 8.785903 | 11.14347 | 8.485555 | 10.412894 | 8.189842 | 10.296493 | 8.159750 | 10.280791 | 8.159321 | 10.781402 | 8.458943 | 10.425273 | 8.240900 | 10.331761 | 8.397134 | 10.466147 | 8.421420 | 10.163560 | 7.907398 | 9.883056 | 8.210950 | 10.372989 | 8.587317 | 10.701480 | 8.335644 | 10.221055 | 8.471277 | 11.556073 | 8.307837 | 11.93393 | 8.305307 | 11.445798 | 8.478097 | 10.881920 | 8.466341 | 10.564357 | 8.595947 | 10.852931 | 9.207234 | 11.698324 | 8.398851 | 9.424032 | 8.181294 | 8.731565 | 8.342664 | 9.416161 | 8.281758 | 11.181796 | 8.360702 | 8.514954 | 8.423486 | 9.84759 | 8.694980 | 9.129945 | 8.386976 | 11.067119 | 8.327773 | 11.114540 | 8.555494 | 11.065938 | 8.289550 | 10.367968 | 8.561558 | 9.485563 | 9.621530 | 10.95606 | 8.125341 | 10.319693 |
| SPP1 | 17975.1510 | 4.417666 | 0.2693309 | 16.40237 | 0 | 0 | 9.802055 | 14.176006 | 10.870985 | 16.825314 | 8.688407 | 12.657360 | 10.557322 | 14.38728 | 9.522463 | 10.686337 | 10.352215 | 12.227207 | 7.823136 | 14.615519 | 8.604498 | 15.602296 | 7.498959 | 13.301173 | 10.576368 | 14.946618 | 8.665934 | 13.679833 | 8.535472 | 12.638242 | 9.390863 | 16.252230 | 10.644801 | 15.998609 | 13.139803 | 16.993093 | 10.831879 | 17.094131 | 9.481431 | 14.384866 | 11.619320 | 16.09676 | 7.953639 | 11.337801 | 8.943416 | 12.980642 | 9.274685 | 14.442331 | 10.909059 | 12.623616 | 8.333667 | 9.721764 | 9.142822 | 15.083419 | 9.177301 | 12.935540 | 8.667962 | 13.620272 | 9.482104 | 14.573337 | 16.849238 | 15.089669 | 9.714998 | 14.94413 | 9.434579 | 13.164718 | 7.680159 | 10.922952 | 9.831650 | 12.937267 | 8.243920 | 13.900129 | 8.454609 | 16.218752 | 9.710107 | 14.391171 | 14.844814 | 14.09893 | 8.620635 | 11.166524 |
| AGER | 11559.6967 | -5.477084 | 0.3355681 | -16.32183 | 0 | 0 | 15.168655 | 6.893846 | 12.872945 | 7.774126 | 14.211504 | 12.702896 | 14.940987 | 11.01977 | 13.667223 | 12.246874 | 14.571089 | 12.406805 | 15.184985 | 8.875605 | 14.333295 | 9.462969 | 14.876306 | 7.152258 | 13.901056 | 11.208358 | 15.016984 | 7.368619 | 14.376976 | 7.482998 | 14.913960 | 8.644041 | 14.072504 | 6.911631 | 15.072956 | 6.885683 | 13.749445 | 9.719213 | 13.588306 | 8.179431 | 14.250593 | 8.67864 | 13.218267 | 9.261431 | 14.676718 | 9.386520 | 12.926690 | 6.397908 | 14.500176 | 7.615689 | 15.010608 | 7.614820 | 14.180083 | 13.422113 | 14.567144 | 7.235510 | 13.978117 | 13.814287 | 14.741739 | 8.525701 | 10.651539 | 8.721279 | 14.509769 | 12.47292 | 13.486554 | 11.975142 | 14.101910 | 8.018593 | 14.218300 | 6.865988 | 13.891117 | 6.501559 | 14.948956 | 6.406268 | 14.411703 | 10.787851 | 12.133789 | 10.50014 | 15.309568 | 7.681559 |
| KRT6A | 101501.5369 | 6.204478 | 0.3804130 | 16.30985 | 0 | 0 | 9.193072 | 18.281839 | 9.516516 | 16.677201 | 5.958339 | 6.729277 | 12.537963 | 17.66984 | 6.000258 | 9.543121 | 6.498843 | 9.294537 | 6.309422 | 7.986419 | 7.160003 | 12.763140 | 5.583972 | 18.514452 | 6.573442 | 8.585020 | 5.795055 | 10.728299 | 5.594051 | 8.746860 | 7.648745 | 11.048317 | 5.015499 | 8.242924 | 10.509087 | 14.326646 | 7.329671 | 16.666850 | 10.764745 | 19.773745 | 12.810230 | 19.99644 | 8.846682 | 17.016851 | 7.263105 | 16.482515 | 8.695843 | 7.076781 | 8.834410 | 18.518557 | 9.164710 | 15.061089 | 7.451653 | 12.464818 | 7.733001 | 12.179191 | 8.488879 | 13.376917 | 8.754957 | 11.945147 | 12.481967 | 15.021821 | 8.826293 | 16.07453 | 7.053591 | 12.239078 | 7.148210 | 19.568155 | 12.172154 | 19.513793 | 6.482059 | 19.631411 | 10.718135 | 19.547211 | 10.227181 | 17.802126 | 10.890398 | 15.85074 | 8.951693 | 17.814489 |
| MMP12 | 1248.2652 | 4.513220 | 0.2774515 | 16.26670 | 0 | 0 | 8.657790 | 12.099376 | 6.731019 | 9.882017 | 5.567395 | 7.050066 | 6.715331 | 10.78122 | 7.200864 | 9.530651 | 8.358858 | 9.221172 | 5.939377 | 12.950840 | 6.216791 | 10.401415 | 4.961707 | 13.143017 | 6.620989 | 9.432495 | 6.324553 | 11.631126 | 6.749005 | 9.939984 | 8.058865 | 10.756126 | 7.296140 | 9.847195 | 6.777252 | 9.647318 | 6.496149 | 9.749994 | 6.683024 | 13.499608 | 6.954529 | 11.63506 | 7.290949 | 12.360101 | 5.156713 | 10.941399 | 6.538170 | 11.154319 | 6.740462 | 10.804220 | 7.485237 | 12.103772 | 6.245587 | 10.292210 | 5.864011 | 8.131906 | 5.810554 | 9.550398 | 5.351996 | 7.667948 | 10.131743 | 9.428738 | 7.772350 | 10.98200 | 7.189583 | 9.352967 | 5.697359 | 9.779781 | 6.557662 | 10.953639 | 5.408323 | 11.284251 | 5.870229 | 11.361134 | 5.962280 | 12.491766 | 8.308241 | 11.67313 | 5.158645 | 7.796870 |
| KRT16 | 6223.7362 | 5.838103 | 0.3652821 | 15.98245 | 0 | 0 | 5.963193 | 13.704380 | 7.391931 | 12.286157 | 4.995210 | 6.307710 | 7.257980 | 12.00832 | 5.181189 | 6.695150 | 5.090391 | 7.000848 | 5.670905 | 9.215005 | 5.204919 | 11.382030 | 5.150114 | 12.312411 | 5.220550 | 6.160121 | 4.659063 | 8.918967 | 5.695141 | 9.258571 | 5.288214 | 8.678109 | 5.015499 | 10.872856 | 7.620772 | 10.344609 | 5.949783 | 12.881752 | 6.488183 | 14.915779 | 7.121346 | 12.84905 | 6.890291 | 15.232684 | 5.399736 | 9.331258 | 5.818049 | 5.854666 | 5.921749 | 13.815014 | 7.237825 | 11.726095 | 5.990253 | 7.302332 | 4.651784 | 7.213946 | 5.061602 | 9.060504 | 6.179562 | 9.365155 | 7.727509 | 8.901664 | 6.032531 | 13.42783 | 6.005928 | 7.697149 | 5.955840 | 17.031628 | 7.900307 | 11.917893 | 5.012995 | 16.420047 | 6.995940 | 16.074262 | 6.230885 | 11.809046 | 9.548481 | 12.06656 | 6.189654 | 12.323947 |
Heatmap visualisation
Heatmaps are great for showing high dimensional data.
topgenes <- head(rownames(dge),50)
colfunc <- colorRampPalette(c("blue", "white", "red"))
rpm <- apply(xf,2,function(x) {x/sum(x) * 1e6 } )
top <- rpm[which(rownames(rpm) %in% topgenes),]
csc = as.character(ss$tumor)
heatmap.2(as.matrix(top),trace="none",col=colfunc(25),scale="row",
margins = c(6,6), cexRow=0.6, cexCol=0.3, ColSideColors=csc )
You may want to look at a few genes in more detail.
myrows <- c(1,2,3,5,6,38)
par(mfrow=c(2,3))
nullplot <- lapply(myrows,function(i) {
nam=rownames(top)[i]
gdata <- top[i,]
ndat <- gdata[grep("N",names(gdata))]
tdat <- gdata[grep("T",names(gdata))]
plotdat <- list("non-tumor"=ndat,"tumor"=tdat)
boxplot(plotdat,ylab="normalised expression (RPM)",cex=0,col="white",main=nam)
beeswarm(plotdat,add=TRUE,cex=0.9,pch=19)
})
nullplot <- lapply(myrows,function(i) {
nam=rownames(top)[i]
gdata <- log10(top[i,]+0.1)
ndat <- gdata[grep("N",names(gdata))]
tdat <- gdata[grep("T",names(gdata))]
plotdat <- list("non-tumor"=ndat,"tumor"=tdat)
boxplot(plotdat,ylab="log10 (RPM)",cex=0,col="white",main=nam)
beeswarm(plotdat,add=TRUE,cex=0.9,pch=19)
})
par(mfrow=c(1,1))Prep data for ML
Using the top 50 DESEQ2 genes as input to train xgboost model.
I will reserve 20 patients for the test. The rest will be used in the training.
head(top)## X10N X10T X11N X11T X12N
## AGER 1494.077763 4.086579 327.76844018 7.4335803 742.54567497
## B4GALNT4 0.119942 23.251226 1.09663871 5.6512599 0.04248459
## CALML3 0.239884 51.058753 0.04217841 0.5216548 0.00000000
## CDC6 2.638724 58.903106 8.64657447 40.1239452 3.61119020
## CLDN18 179.953011 4.227496 30.15756463 0.9563671 134.84609065
## CLIC5 293.458142 2.019803 182.37945378 3.9558819 264.72148419
## X12T X14N X14T X15N X15T
## AGER 263.55858967 1319.5604666 83.824302 549.28380264 203.7734137
## B4GALNT4 0.55468054 0.1314216 4.366323 0.27900636 0.2259129
## CALML3 0.08533547 0.5256864 17.010466 0.03985805 0.1355477
## CDC6 4.39477655 3.2417328 31.064568 2.79006358 6.2803779
## CLDN18 43.05174308 174.3526544 5.685316 76.84632258 67.0509415
## CLIC5 69.42040236 458.1356935 35.840233 244.40956954 88.3319343
## X16N X16T X17N X17T X18N X18T
## AGER 1.010643e+03 209.3178773 1.542161e+03 17.560922 854.9968712 42.55561
## B4GALNT4 7.160571e-02 0.5634398 8.048458e-01 11.734635 0.2055182 22.27270
## CALML3 3.580285e-02 0.2817199 3.832599e-02 1.887389 0.0000000 10.01309
## CDC6 1.754340e+00 21.1289916 2.337886e+00 75.741733 1.6030421 44.80214
## CLDN18 1.180062e+02 20.2838320 1.651850e+02 3.364476 136.9573374 13.47915
## CLIC5 2.660510e+02 64.3126260 3.219000e+02 10.339608 340.6258868 28.62716
## X19N X19T X20N X20T X21N X21T
## AGER 1.221352e+03 5.229616 693.84458933 105.0777079 1440.7305118 5.503928
## B4GALNT4 1.607941e-01 11.825527 1.02405615 1.5807762 0.3653524 6.646253
## CALML3 8.039707e-02 2.308569 0.04654801 0.1756418 0.4018876 9.813608
## CDC6 1.849133e+00 76.748144 3.16526446 8.6942692 3.2151008 59.712429
## CLDN18 2.835605e+02 2.402796 81.92449196 16.3346876 158.7090669 1.090401
## CLIC5 5.148628e+02 2.826819 220.68410023 51.5508688 362.2103336 61.685535
## X22N X22T X23N X23T X24N X24T
## AGER 857.57599073 5.232740 1.296817e+03 16.2089420 751.6361546 3.321511
## B4GALNT4 0.48670601 16.977336 3.557310e-01 2.3398735 0.2081979 25.217683
## CALML3 0.08111767 16.589725 4.446637e-02 0.2127158 0.0000000 1.741375
## CDC6 4.21811876 73.297128 5.069167e+00 4.3394018 2.3318168 70.751402
## CLDN18 104.52011578 2.170618 1.012499e+02 5.7433259 135.6201294 6.610774
## CLIC5 254.70947890 1.783008 3.338535e+02 34.4174123 225.1035982 1.031926
## X25N X25T X26N X26T X27N X27T
## AGER 1470.4860724 3.2154525 574.8708877 37.008149 471.939498 9.847035
## B4GALNT4 1.3673606 17.7445340 0.2325153 9.374743 0.068556 8.692804
## CALML3 0.2646504 0.4763633 0.0000000 19.780218 0.102834 24.743831
## CDC6 4.2785155 71.2560145 7.3009812 41.572815 1.782456 18.287350
## CLDN18 107.7568383 0.1984847 75.7999957 5.300902 185.478256 4.400506
## CLIC5 285.5137185 0.1984847 254.7437893 56.248460 210.775419 9.089571
## X28N X28T X29N X29T X2N X2T
## AGER 847.4133996 12.064514 385.2199594 22.706520 1061.4136417 26.049411
## B4GALNT4 0.8650018 10.690075 1.3822649 13.616109 0.1628872 14.784801
## CALML3 0.9370852 61.429754 0.2159789 20.833818 0.0000000 23.316087
## CDC6 7.6408490 78.304801 2.2461805 31.445800 2.7283604 12.258546
## CLDN18 138.6525752 3.856063 31.8352893 7.724899 151.2814763 3.147465
## CLIC5 292.6949736 5.994078 177.4914570 8.973367 319.6661073 19.216100
## X30N X30T X31N X31T X32N X32T
## AGER 314.84648515 2.2078301 1.066893e+03 6.2318373 1293.7659660 7.238105
## B4GALNT4 0.92291683 7.8909482 3.429299e-01 8.3860526 0.4206001 4.894833
## CALML3 0.06592263 0.2862002 7.620664e-02 36.6985973 0.7570802 82.951805
## CDC6 4.97715863 48.1634043 2.248096e+00 43.8921378 6.2669418 36.086379
## CLDN18 71.13051873 2.9028877 9.769691e+01 0.8462989 123.6564350 1.353890
## CLIC5 137.61349199 1.2265723 3.275361e+02 2.8081736 376.4371065 13.747192
## X33N X33T X34N X34T X35N X35T
## AGER 813.37729515 446.292046 1.072677e+03 3.8025254 677.39948856 599.752392
## B4GALNT4 0.30831644 4.541829 1.795696e-01 0.1810726 0.16889386 3.186625
## CALML3 0.03853955 1.009295 7.182784e-02 0.1086436 0.03377877 1.519205
## CDC6 2.31237326 14.497151 2.837200e+00 2.1366571 5.20193085 5.891550
## CLDN18 83.97768899 35.875861 2.056431e+02 1.1588649 131.60209472 61.842749
## CLIC5 240.02434467 116.252472 3.413259e+02 1.5210101 231.08057750 160.628112
## X36N X36T X37N X37T X38N X38T
## AGER 1079.9607244 9.522532 63.90353350 14.8263543 925.5271908 244.324417
## B4GALNT4 0.2358422 2.426414 0.25841984 0.5702444 0.2283561 5.698123
## CALML3 0.2358422 1.098754 0.07383424 5.0100043 0.0000000 10.979311
## CDC6 2.4370360 11.628476 9.59845102 9.5719595 2.6946017 18.484157
## CLDN18 211.5897572 0.915628 8.78627439 6.9651280 114.8174369 27.274555
## CLIC5 598.3709593 5.219080 95.98451017 24.2761185 231.0963526 89.745445
## X39N X39T X3N X3T X4N X4T
## AGER 417.35375051 150.5498329 667.02379195 9.461049 765.008507 3.5831458
## B4GALNT4 0.37629100 1.0479107 0.45975218 7.065333 1.475510 11.8789074
## CALML3 0.09407275 0.2716805 0.03831268 43.082286 0.000000 0.5452613
## CDC6 6.16176505 9.8193111 2.06888482 27.692855 3.124609 30.9630538
## CLDN18 58.27806795 15.4857910 84.21127482 1.624214 81.283239 0.7399975
## CLIC5 175.53974945 86.6272820 216.27508935 5.725356 268.542809 1.7915729
## X5N X5T X6N X6T X7N X7T
## AGER 663.65237509 2.144060 1217.1591347 1.840647 952.1658869 72.677722
## B4GALNT4 0.08286333 2.728804 0.0000000 7.202533 0.5594900 3.553317
## CALML3 0.00000000 103.499632 0.1248496 100.515353 0.3012638 32.434344
## CDC6 1.53297152 48.728641 1.6230453 32.171315 2.9265630 54.208738
## CLDN18 129.34965133 1.130504 171.7931038 3.161112 82.0728777 11.610256
## CLIC5 223.23379929 0.896607 367.0163233 2.240788 241.1832232 25.740887
## X8N X8T X9N X9T
## AGER 201.1404097 63.399919 1552.9668318 6.903154
## B4GALNT4 3.6730465 5.532435 0.7493524 5.108334
## CALML3 0.8771454 19.898919 0.3085569 55.547377
## CDC6 13.2120029 36.719305 2.7329325 50.669148
## CLDN18 23.4088185 6.558612 133.7814520 1.886862
## CLIC5 184.7487546 29.536063 382.4783060 4.418018my_data <- as.data.frame(t(top))
my_labels <- as.numeric(grepl("T",rownames(my_data)))
Ntest=20
Ntrain=nrow(my_data) - Ntest
test_data <- my_data[1:Ntest,]
test_labels <- my_labels[1:Ntest]
train_data <- my_data[(Ntest+1):nrow(my_data),]
train_labels <- my_labels[(Ntest+1):length(my_labels)]
dim(test_data)## [1] 20 50dim(train_data)## [1] 54 50length(test_labels)## [1] 20length(train_labels)## [1] 54Train model
model <- xgboost(x = train_data, y = train_labels, nrounds = 100, objective = "reg:squarederror")
model## XGBoost model object
## Call:
## xgboost(x = train_data, y = train_labels, objective = "reg:squarederror",
## nrounds = 100)
## Objective: reg:squarederror
## Number of iterations: 100
## Number of features: 50Predict
pred <- predict(model, test_data)
prediction <- as.numeric(pred > 0.5)
table(test_labels == prediction)##
## FALSE TRUE
## 1 19table(test_labels == prediction)[2] / length(prediction)## TRUE
## 0.95Calculate shap values
iml package is used to calculate SHAP values.
# Create a Predictor object
predictor <- iml::Predictor$new(model, data = train_data)
# Example for how to calculate SHAP values for one observations
shapley <- iml::Shapley$new(predictor, x.interest = train_data[1, ])
shapley## Interpretation method: Shapley
## Predicted value: 0.000138, Average prediction: 0.500000 (diff = -0.499862)
##
## Analysed predictor:
## Prediction task: unknown
##
##
## Analysed data:
## Sampling from data.frame with 54 rows and 50 columns.
##
##
## Head of results:
## feature phi phi.var feature.value
## 1 AGER 0 0 AGER=1440.73051182176
## 2 B4GALNT4 0 0 B4GALNT4=0.365352363904692
## 3 CALML3 0 0 CALML3=0.401887600295161
## 4 CDC6 0 0 CDC6=3.21510080236129
## 5 CLDN18 0 0 CLDN18=158.709066880198
## 6 CLIC5 0 0 CLIC5=362.210333575111# Now calculate all SHAP values
# mclapply causes this to freeze, so we need to do it single threaded for now
shapley <- lapply(1:nrow(train_data), function(i) {
iml::Shapley$new(predictor, x.interest = train_data[i, ])
})
shapley[[1]]## Interpretation method: Shapley
## Predicted value: 0.000138, Average prediction: 0.500000 (diff = -0.499862)
##
## Analysed predictor:
## Prediction task: unknown
##
##
## Analysed data:
## Sampling from data.frame with 54 rows and 50 columns.
##
##
## Head of results:
## feature phi phi.var feature.value
## 1 AGER 0 0 AGER=1440.73051182176
## 2 B4GALNT4 0 0 B4GALNT4=0.365352363904692
## 3 CALML3 0 0 CALML3=0.401887600295161
## 4 CDC6 0 0 CDC6=3.21510080236129
## 5 CLDN18 0 0 CLDN18=158.709066880198
## 6 CLIC5 0 0 CLIC5=362.210333575111shap_list <- lapply(1:length(shapley), function(i) { shapley[[i]]$results })
shap_df <- do.call(rbind, shap_list)
meanabs <- function(x) { mean(abs(x)) }
shap_dfa <- aggregate(phi ~ feature,shap_df,FUN=meanabs )
shap_dfa <- shap_dfa[order(-shap_dfa$phi),]
head(shap_dfa)## feature phi
## 38 SERPINE2 0.5113404
## 1 AGER 0.0000000
## 2 B4GALNT4 0.0000000
## 3 CALML3 0.0000000
## 4 CDC6 0.0000000
## 5 CLDN18 0.0000000In this simple model, only the SERPINE2 gene was used to classify the tumor and control samples. So it is not a great example.
Remove SERPINE2 from data
Let’s see if the model is different or performs worse than the original one.
train_data2 <- train_data[,-which(colnames(train_data)=="SERPINE2")]
model2 <- xgboost(x = train_data2, y = train_labels, nrounds = 100, objective = "reg:squarederror")
test_data2 <- test_data[,which(colnames(test_data)!="SERPINE2")]
pred2 <- predict(model2, test_data2)
prediction2 <- as.numeric(pred2 > 0.5)
table(test_labels == prediction2)##
## FALSE TRUE
## 1 19table(test_labels == prediction2)[2] / length(prediction2)## TRUE
## 0.95predictor2 <- iml::Predictor$new(model2, data = train_data2)
shapley2 <- iml::Shapley$new(predictor2, x.interest = train_data2[1, ])
shapley2 <- lapply(1:nrow(train_data2), function(i) {
iml::Shapley$new(predictor2, x.interest = train_data2[i, ])
})
shapley2[[1]]## Interpretation method: Shapley
## Predicted value: 0.000012, Average prediction: 0.500000 (diff = -0.499988)
##
## Analysed predictor:
## Prediction task: unknown
##
##
## Analysed data:
## Sampling from data.frame with 54 rows and 49 columns.
##
##
## Head of results:
## feature phi phi.var feature.value
## 1 AGER 2.534488e-02 1.254827e-03 AGER=1440.73051182176
## 2 B4GALNT4 0.000000e+00 0.000000e+00 B4GALNT4=0.365352363904692
## 3 CALML3 -1.397973e-01 6.718077e-02 CALML3=0.401887600295161
## 4 CDC6 0.000000e+00 0.000000e+00 CDC6=3.21510080236129
## 5 CLDN18 -7.258609e-05 1.130917e-08 CLDN18=158.709066880198
## 6 CLIC5 -4.950514e-02 2.397352e-03 CLIC5=362.210333575111shap_list2 <- lapply(1:length(shapley2), function(i) { shapley2[[i]]$results })
shap_df2 <- do.call(rbind, shap_list2)
shap_dfa2 <- aggregate(phi ~ feature,shap_df2,FUN=meanabs )
shap_dfa2 <- shap_dfa2[order(-shap_dfa2$phi),]
head(shap_dfa2)## feature phi
## 36 RGCC 0.25056379
## 3 CALML3 0.14908066
## 24 LMO7 0.04556933
## 6 CLIC5 0.03475354
## 1 AGER 0.01991069
## 48 ZIC2 0.01815962Barplot.
modelgenes <- head(shap_dfa2,10)$feature
modelweights <- head(shap_dfa2,10)$phi
names(modelweights) <- modelgenes
par(mar=c(5.1,8.1,4.1,2.1))
barplot(rev(modelweights),horiz=TRUE,las=1,xlab="Mean Abs Phi values",main="Model weights")
par(mar=c(5.1,4.1,4.1,2.1))Now look at the direction of the SHAP values as well.
SHAP values for normal and cancer samples are separated for charting.
shap_l <- lapply(topgenes,function(x) { shap_df2[which(shap_df2$feature == x),"phi"] } )
names(shap_l) <- topgenes
which(lapply(shap_l,meanabs)>0)## AGER CLDN18 SLC6A4 ZIC2 CLIC5 CALML3 RGCC COL11A1 ZYG11A LMO7
## 3 7 11 13 14 15 16 21 40 41
## CXCL13
## 42mytopmodelgenes <- names(which(lapply(shap_l,meanabs)>0))
shap_top <- shap_l[which(names(shap_l) %in% mytopmodelgenes)]
shap_top <- shap_top[order(unlist(lapply(shap_top,median)))]
shap_top## $RGCC
## [1] -0.2407864 0.2521787 -0.2649532 0.2180791 -0.2705543 0.4060387
## [7] -0.2729865 0.2165776 -0.1539801 0.4293343 -0.2724211 0.1747100
## [13] -0.2352874 0.1992804 -0.3080797 0.2582473 -0.2442303 0.2228216
## [19] -0.2507993 0.2403487 -0.3237022 0.3379525 -0.2244573 0.2710644
## [25] -0.2970060 0.1922484 -0.2225132 0.2668904 -0.2520664 0.3541977
## [31] -0.1619480 -0.1568924 -0.2427500 0.1776787 -0.2363281 0.2838249
## [37] -0.2720199 0.2485171 -0.2422606 0.3277131 -0.2501515 0.1397674
## [43] -0.2642617 0.4584410 -0.2525709 0.2020575 -0.2654936 0.1682303
## [49] -0.1969997 0.2160078 -0.1888147 0.1720752 -0.3183917 0.2134554
##
## $CALML3
## [1] -0.13979727 0.14813166 -0.14318056 0.15702752 -0.15192546 0.03860351
## [7] -0.22138802 0.15665894 -0.10613308 0.01534872 -0.19672475 0.14661209
## [13] -0.08318746 0.20924870 -0.13083176 0.16276811 -0.08283436 0.16437612
## [19] -0.21280345 0.15543663 -0.23137932 0.02188237 -0.16816999 0.17614626
## [25] -0.09606590 0.17063389 -0.27665558 0.18776720 -0.19711466 0.02494426
## [31] -0.13907879 0.39495030 -0.10589949 0.20229389 -0.22572046 0.16669855
## [37] -0.19115420 0.16196853 -0.16620642 0.01947839 -0.15462876 0.17823009
## [43] -0.11986292 0.02491817 -0.12621228 0.18144870 -0.11307529 0.16480703
## [49] -0.14388585 0.23162008 -0.08089183 0.13169108 -0.07769501 0.17416175
##
## $ZIC2
## [1] -0.0235415364 0.0238905209 -0.0168152150 0.0165970146 -0.0179163526
## [6] 0.0214318594 -0.0211846924 0.0204135452 -0.0170846679 0.0237566776
## [11] -0.0242928201 0.0165688127 -0.0181151489 0.0175399722 -0.0253713594
## [16] 0.0227986930 -0.0156320364 0.0181744383 -0.0186031237 0.0149616891
## [21] -0.0222558071 0.0181490274 -0.0184731735 0.0226752693 -0.0148793043
## [26] 0.0154149177 -0.0250346880 -0.0078513533 -0.0199920456 -0.0008329308
## [31] -0.0135874183 0.0432903239 -0.0165844343 0.0223611574 -0.0082983061
## [36] -0.0007387984 -0.0223033137 0.0287469481 -0.0158134279 -0.0013533914
## [41] -0.0185160026 0.0174311954 -0.0166126830 0.0207517629 -0.0170176879
## [46] 0.0160869518 -0.0213763536 0.0150111702 -0.0192955365 0.0229758219
## [51] -0.0062328353 0.0195625153 -0.0206694438 0.0157530874
##
## $CXCL13
## [1] -1.232445e-03 -3.728884e-04 -1.836282e-03 -4.095501e-04 -3.731544e-03
## [6] -1.640526e-03 -8.732468e-04 -2.979636e-06 -1.591905e-03 -7.724327e-04
## [11] -3.355788e-03 -9.168720e-04 -8.104968e-04 -5.219352e-04 -2.036093e-03
## [16] -6.259084e-06 -8.262901e-04 -1.119995e-04 -2.393708e-03 8.571143e-05
## [21] -1.608880e-03 -4.113537e-04 -8.389076e-04 -8.796871e-04 -6.699682e-04
## [26] 9.708002e-03 -4.414516e-03 -1.052804e-03 -8.641490e-04 8.513331e-06
## [31] -8.217273e-04 2.872030e-02 -1.350893e-03 -5.177170e-04 -3.983798e-04
## [36] 1.239300e-05 -3.598327e-03 -1.459318e-03 -2.405909e-03 -7.703120e-04
## [41] -2.620986e-05 9.982319e-05 -4.805533e-04 -8.234417e-04 -2.029460e-03
## [46] -4.689097e-06 -3.591158e-03 -5.408525e-06 -3.873571e-03 9.741978e-03
## [51] -8.159376e-04 7.289055e-03 -2.097562e-03 -1.119930e-04
##
## $AGER
## [1] 0.0253448777 -0.0243708764 0.0185392065 -0.0165518838 0.0212022832
## [6] -0.0244366680 0.0213154887 -0.0295075795 0.0173345099 -0.0248998409
## [11] 0.0266160755 -0.0192928874 0.0214322793 -0.0218555058 0.0305393407
## [16] -0.0225309176 0.0105198589 -0.0201183398 0.0179940335 -0.0164799356
## [21] 0.0158133236 -0.0183671205 0.0187307695 -0.0323571354 0.0207848965
## [26] -0.0259798391 0.0232347014 -0.0007016947 0.0233892185 -0.0085334017
## [31] 0.0125971276 0.0574680254 0.0207344621 -0.0243920662 -0.0060631073
## [36] -0.0060359219 0.0205438657 -0.0173798058 0.0095364162 -0.0051847643
## [41] 0.0197433776 -0.0171398331 0.0175304343 -0.0267031699 0.0182108082
## [46] -0.0166440486 0.0268649141 -0.0147245719 0.0189148396 -0.0312102191
## [51] -0.0006035946 -0.0259424798 0.0269327875 -0.0153019764
##
## $COL11A1
## [1] 1.440428e-04 -9.132223e-05 -1.545089e-05 -1.598905e-04 7.815566e-05
## [6] -1.149365e-04 1.135849e-04 -1.703415e-04 -2.302444e-06 -1.998468e-04
## [11] 1.008698e-04 -1.026616e-04 1.167659e-04 -9.448805e-05 1.353961e-04
## [16] -1.108280e-04 -3.087872e-05 -1.476298e-04 -4.589507e-05 -1.140402e-04
## [21] -5.542764e-05 -6.585182e-05 -3.407445e-05 -1.299089e-04 1.417487e-04
## [26] -4.632339e-05 1.353996e-04 2.466605e-04 1.340185e-04 3.174861e-04
## [31] 9.269189e-05 1.507362e-04 -2.728331e-05 -1.031213e-04 1.271775e-04
## [36] 3.334005e-04 9.722371e-05 -1.222187e-04 9.810491e-05 3.165768e-04
## [41] 1.149404e-04 -1.122011e-04 -1.545660e-05 -1.621795e-04 -1.909154e-05
## [46] -1.121872e-04 -4.136197e-05 -1.494643e-04 -1.638444e-05 -6.995406e-05
## [51] -7.629901e-05 -7.994741e-05 -3.683405e-05 -1.054065e-04
##
## $CLDN18
## [1] -7.258609e-05 9.073373e-06 -4.537106e-06 6.803274e-06 -2.041618e-05
## [6] 1.134071e-05 -3.175452e-05 2.268109e-05 -6.804466e-06 2.721847e-05
## [11] -2.722189e-05 2.268553e-06 -6.577782e-05 1.134276e-05 -8.165778e-05
## [16] 1.587689e-05 1.814507e-05 2.268073e-05 -1.361012e-05 6.805239e-06
## [21] -2.268553e-06 9.072199e-06 -2.268106e-06 1.587742e-05 -8.165787e-05
## [26] 2.041374e-05 -2.948724e-05 6.351542e-05 -2.722102e-05 6.351527e-05
## [31] -2.041572e-05 5.217328e-05 -1.587871e-05 9.072600e-06 6.804373e-06
## [36] 3.175974e-05 -2.495128e-05 9.074211e-06 1.360882e-05 4.763654e-05
## [41] -2.722018e-05 1.134157e-05 -2.268553e-06 2.041574e-05 -1.134239e-05
## [46] 1.134192e-05 -4.537106e-06 1.134157e-05 -1.134127e-05 1.360934e-05
## [51] 1.587693e-05 9.072600e-06 -1.360851e-05 9.073792e-06
##
## $SLC6A4
## [1] 0.000000e+00 0.000000e+00 1.116492e-04 1.315114e-04 1.105748e-04
## [6] -1.364291e-03 5.369373e-06 4.383624e-05 7.318005e-05 7.693737e-05
## [11] 7.479026e-05 0.000000e+00 0.000000e+00 4.920512e-05 7.693826e-05
## [16] 4.383712e-05 1.078908e-04 3.846824e-05 2.483439e-04 1.154029e-04
## [21] 6.942181e-05 0.000000e+00 1.105757e-04 2.684705e-06 1.528025e-04
## [26] 0.000000e+00 7.747450e-05 8.051991e-06 7.693826e-05 6.262243e-05
## [31] 2.157810e-04 4.007876e-05 7.318005e-05 5.367994e-06 -3.045649e-03
## [36] 1.646131e-04 1.479707e-04 2.334961e-04 7.318006e-05 1.154047e-04
## [41] 1.116492e-04 1.609862e-05 7.747379e-05 1.368799e-04 1.528029e-04
## [46] 5.367994e-06 8.052588e-06 5.367994e-06 7.586418e-05 7.962048e-05
## [51] 1.116495e-04 2.683997e-06 7.854864e-05 7.962137e-05
##
## $ZYG11A
## [1] -0.0038973361 0.0042977760 -0.0046777074 0.0040812272 -0.0038412327
## [6] 0.0033832370 -0.0062739558 0.0038063206 -0.0037634685 0.0043372980
## [11] -0.0048145124 0.0032732971 -0.0044628206 0.0040062026 -0.0058897925
## [16] 0.0039311733 0.0006370734 0.0039623627 -0.0040564099 0.0042462746
## [21] 0.0004723142 0.0031906206 -0.0049782067 0.0043011032 -0.0045392219
## [26] 0.0036190065 -0.0058735161 0.0043771196 -0.0056955136 0.0040224835
## [31] -0.0036429374 0.0056742380 -0.0043214043 0.0033480286 -0.0062479651
## [36] 0.0032902640 -0.0045555084 0.0055171253 -0.0067753144 0.0030684245
## [41] -0.0041460502 0.0032215128 -0.0045823829 0.0049250523 -0.0055968577
## [46] 0.0040284644 -0.0041580134 0.0028781555 -0.0046365143 0.0047175263
## [51] 0.0014320585 0.0032025873 -0.0060740759 0.0046128902
##
## $CLIC5
## [1] -0.049505143 0.033995307 -0.053290648 0.028467598 -0.026539660
## [6] 0.027098790 -0.029376160 0.035440419 -0.035376150 0.038985180
## [11] -0.031236082 0.029164878 -0.041223924 0.039599458 -0.059634618
## [16] 0.039787998 -0.040279767 0.034547790 -0.026188004 0.032327517
## [21] 0.005625716 0.030139203 -0.047153206 0.039743054 -0.054954686
## [26] 0.031330806 -0.031165757 0.036152995 -0.023380448 0.033886301
## [31] -0.017385875 0.031151661 -0.044380180 0.033315184 0.019451776
## [36] 0.043513149 -0.026825693 0.038109284 -0.051914021 0.024942475
## [41] -0.026201767 0.026500517 -0.020065707 0.041188873 -0.022738170
## [46] 0.037758727 -0.049606030 0.028273563 -0.046568883 0.034584229
## [51] -0.036037091 0.025804780 -0.052904038 0.031871988
##
## $LMO7
## [1] -0.06662532 0.04231502 -0.03393435 0.03269179 -0.03688520 0.04072535
## [7] -0.03927780 0.04674606 -0.04947439 0.05401593 -0.04388753 0.04004138
## [13] -0.05839210 0.05274690 -0.07883105 0.05503275 -0.05715079 0.04642510
## [19] -0.04321231 0.03903896 -0.04250526 0.03751265 -0.03472869 0.05947094
## [25] -0.07294827 0.04310747 -0.04764327 0.05382705 -0.03451894 0.04148403
## [31] -0.02647308 0.04352457 -0.06550859 0.04600500 0.01810082 0.04839463
## [37] -0.04026144 0.04576458 0.01626360 0.03128361 -0.03634017 0.04205336
## [43] -0.03169777 0.04724081 -0.02215323 0.04533064 -0.06950156 0.04569482
## [49] -0.06362376 0.05149320 -0.05780098 0.03641192 -0.06910070 0.03552449boxplot(shap_top,cex=0,horizontal=TRUE,las=1,col="white",main="Phi values for training data")
beeswarm(shap_top,horizontal=TRUE,pch=19,cex=0.4,add=TRUE,col="darkgray")
LEN=54
seq(1,LEN,2)## [1] 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
## [26] 51 53NORM <- lapply(shap_top,function(x) { x[seq(1,LEN,2)] })
CANC <- lapply(shap_top,function(x) { x[seq(1,LEN,2)+1] })
boxplot(shap_top,cex=0,horizontal=TRUE,las=1,col="white",main="Phi values for training data")
mtext("blue=normal, red=cancer sample")
beeswarm(NORM,horizontal=TRUE,pch=19,cex=0.4,add=TRUE,col="blue")
beeswarm(CANC,horizontal=TRUE,pch=19,cex=0.4,add=TRUE,col="red")
From this we can see that of all the genes in this model, RGCC is the most inflential. Higher expression of RGCC is associated with cancer classification. RGCC, LMO7, CLIC5 and CALML3 are associated with cancer, while higher AGER is associated with normal tissue.
Conclusion
Here, we trained an XGboost model on a small number of lung cancer datasets (54) using only the top 50 genes from a limma analysis. Testing on a separate 20 samples found it to be approximately 95% accurate. We used SHAP analysis to understand how strongly the genes influence the model. Initially the model exclusively used SERPINE2 to classify. Another model was generated after removing SERPINE2 and it showed RGCC, CALML3, LMO7, CLIC5, AGER and ZIC2 could also be informative in a model. This second model was also 95% accurate. This small example shows how a panel of genes can be selected from an ML model. The SHAP analysis here can be applied to different ML models such as random forest, support vector machine, etc. This is one method to interpret ML models. This type of approach could be used for various omics data such as proteomics, microbiome, or DNA methylation.
Session information
save.image("lung_cancer_gene_expression_shap.Rdata")
sessionInfo()## R version 4.6.0 (2026-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] C
##
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggplot2_4.0.3 iml_0.11.4
## [3] xgboost_3.2.1.1 beeswarm_0.4.0
## [5] gplots_3.3.0 DESeq2_1.52.0
## [7] SummarizedExperiment_1.42.0 Biobase_2.72.0
## [9] MatrixGenerics_1.24.0 matrixStats_1.5.0
## [11] GenomicRanges_1.64.0 Seqinfo_1.2.0
## [13] IRanges_2.46.0 S4Vectors_0.50.1
## [15] BiocGenerics_0.58.1 generics_0.1.4
## [17] kableExtra_1.4.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.57 bslib_0.11.0
## [4] caTools_1.18.3 lattice_0.22-9 vctrs_0.7.3
## [7] tools_4.6.0 bitops_1.0-9 tibble_3.3.1
## [10] pkgconfig_2.0.3 Matrix_1.7-5 KernSmooth_2.23-26
## [13] checkmate_2.3.4 data.table_1.18.4 RColorBrewer_1.1-3
## [16] S7_0.2.2 lifecycle_1.0.5 compiler_4.6.0
## [19] farver_2.1.2 stringr_1.6.0 textshaping_1.0.5
## [22] codetools_0.2-20 htmltools_0.5.9 sass_0.4.10
## [25] yaml_2.3.12 Formula_1.2-5 pillar_1.11.1
## [28] jquerylib_0.1.4 BiocParallel_1.46.0 DelayedArray_0.38.1
## [31] cachem_1.1.0 abind_1.4-8 parallelly_1.47.0
## [34] Metrics_0.1.4 gtools_3.9.5 tidyselect_1.2.1
## [37] locfit_1.5-9.12 digest_0.6.39 future_1.70.0
## [40] stringi_1.8.7 listenv_0.10.1 dplyr_1.2.1
## [43] fastmap_1.2.0 grid_4.6.0 cli_3.6.6
## [46] SparseArray_1.12.2 magrittr_2.0.5 S4Arrays_1.12.0
## [49] dichromat_2.0-0.1 withr_3.0.2 backports_1.5.1
## [52] scales_1.4.0 rmarkdown_2.31 XVector_0.52.0
## [55] globals_0.19.1 otel_0.2.0 evaluate_1.0.5
## [58] knitr_1.51 viridisLite_0.4.3 rlang_1.2.0
## [61] Rcpp_1.1.1-1.1 glue_1.8.1 xml2_1.5.2
## [64] svglite_2.2.2 rstudioapi_0.18.0 jsonlite_2.0.0
## [67] R6_2.6.1 systemfonts_1.3.2