Differential splicing analysis for proteins regulated by SF3B1 mutations
Junyan Lu
2020-03-17
Last updated: 2020-03-18
Checks: 5 2
Knit directory: Proteomics/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.6.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown is untracked by Git. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
The global environment had objects present when the code in the R Markdown file was run. These objects can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment. Use wflow_publish or wflow_build to ensure that the code is always run in an empty environment.
The following objects were defined in the global environment when these results were created:
| Name | Class | Size |
|---|---|---|
| myPaths | character | 608 bytes |
The command set.seed(20200227) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 8ba49d7. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: analysis/.Rhistory
Ignored: analysis/compareProteomicsRNAseq_cache/
Ignored: analysis/correlateCLLPD_cache/
Ignored: analysis/correlateGenomic_cache/
Ignored: analysis/correlateGenomic_noBlock_MCLL_cache/
Ignored: analysis/correlateGenomic_noBlock_UCLL_cache/
Ignored: analysis/correlateGenomic_noBlock_cache/
Ignored: analysis/correlateGenomic_removePC_cache/
Ignored: analysis/correlateMIR_cache/
Ignored: analysis/correlateMethylationCluster_cache/
Ignored: analysis/predictOutcome_cache/
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: ._.DS_Store
Untracked: analysis/._.DS_Store
Untracked: analysis/analysisSF3B1.Rmd
Untracked: analysis/correlateGenomic_noBlock.Rmd
Untracked: analysis/correlateGenomic_noBlock_MCLL.Rmd
Untracked: analysis/correlateGenomic_noBlock_UCLL.Rmd
Untracked: code/utils.R
Untracked: data/._.DS_Store
Untracked: data/190909_CLL_prot_abund_med_norm.tsv
Untracked: data/190909_CLL_prot_abund_no_norm.tsv
Untracked: data/20190423_Proteom_submitted_samples_bereinigt.xlsx
Untracked: data/20191025_Proteom_submitted_samples_final.xlsx
Untracked: data/LUMOS/
Untracked: data/LUMOS_peptides/
Untracked: data/LUMOS_protAnnotation.csv
Untracked: data/SampleAnnotation_cleaned.xlsx
Untracked: data/facTab_IC50atLeast3New.RData
Untracked: data/gmts/
Untracked: data/mapEnsemble.txt
Untracked: data/mapSymbol.txt
Untracked: data/pyprophet_export_aligned.csv
Untracked: data/timsTOF_protAnnotation.csv
Untracked: output/._.DS_Store
Untracked: output/LUMOS_processed.RData
Untracked: output/dxdCLL.RData
Untracked: output/proteomic_LUMOS_20200227.RData
Untracked: output/proteomic_timsTOF_20200227.RData
Untracked: output/timsTOF_processed.RData
Unstaged changes:
Modified: analysis/correlateGenomic.Rmd
Modified: analysis/correlateMIR.Rmd
Modified: analysis/correlateMethylationCluster.Rmd
Modified: analysis/index.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
There are no past versions. Publish this analysis with wflow_publish() to start tracking its development.
Identify differentially expressed proteins related to SF3B1
Process data
protCLL$SF3B1 <- patMeta[match(colnames(protCLL), patMeta$Patient.ID),]$SF3B1
protMat <- assay(protCLL)DE test without blocking for IGHV&trisomy12
SF3B1 <- protCLL$SF3B1
fit <- proDA(protMat, design = ~ SF3B1)
resTab.noBlock <- test_diff(fit, contrast = "SF3B11")
resTab.noBlock$test <- "noBlocking"
hist(resTab.noBlock$pval, breaks =50, main = "P value histogram")
#how may below 10% FDR?
nrow(filter(resTab.noBlock, adj_pval < 0.1))[1] 20DE test with blocking for IGHV&trisomy12
designMat <- data.frame(row.names = colnames(protCLL),
SF3B1 = protCLL$SF3B1,trisomy12 = protCLL$trisomy12, IGHV = protCLL$IGHV.status)
fit <- proDA(protMat, design = ~ SF3B1 + trisomy12 + IGHV, col_data = designMat)
resTab.block <- test_diff(fit, contrast = "SF3B11")
resTab.block$test <- "withBlocking"
hist(resTab.block$pval, breaks =50, main = "P value histogram")
#how may below 10% FDR?
nrow(filter(resTab.block, adj_pval < 0.1))[1] 18Compare two results
resTab.all <- bind_rows(resTab.noBlock, resTab.block) %>%
filter(adj_pval < 0.1) %>%
mutate(direct = ifelse(diff > 0, "up","down")) %>%
mutate(group = paste0(test,"_",direct),
symbol = rowData(protCLL[name,])$hgnc_symbol) %>%
group_by(group) %>% summarise(l = list(symbol))
diffList <- resTab.all$l
names(diffList) <- resTab.all$group
upset(fromList(diffList))
Most cases are consistent
The test result with blocking for IGHV/trisomy12 is used for the analyses below
List of differentially expressed proteins (10% FDR)
resTab <- filter(resTab.block, adj_pval <= 0.1) %>%
mutate(symbol = rowData(protCLL[name,])$hgnc_symbol,
id = rowData(protCLL[name,])$ensembl_gene_id) %>%
dplyr::select(symbol, name,id, pval, adj_pval, diff, n_obs) %>%
mutate(n_obs = as.character(n_obs)) %>%
arrange(pval)
#add ensemble id to SUGP1
resTab[resTab$symbol == "SUGP1", ]$id <- "ENSG00000105705"
resTab %>% mutate_if(is.numeric, formatC, digits =2, format ="e") %>%
DT::datatable()Test if SF3B1 leads to differential splicing of the proteins listed above
Processing splicing dataset
dxdCLL <- dxdCLL[,dxdCLL$diag %in% "CLL"]
dxdCLL$SF3B1 <- factor(patMeta[match(dxdCLL$patID, patMeta$Patient.ID),]$SF3B1)
dxdCLL$trisomy12 <- factor(patMeta[match(dxdCLL$patID, patMeta$Patient.ID),]$trisomy12)
dxdCLL$IGHV <- factor(patMeta[match(dxdCLL$patID, patMeta$Patient.ID),]$IGHV.status)
dxdCLL.sub <- dxdCLL[rowData(dxdCLL)$groupID %in% resTab$id,
!is.na(dxdCLL$SF3B1) & !is.na(dxdCLL$trisomy12) & !is.na(dxdCLL$IGHV)]
#add gene symbol to SUGP1
rowData(dxdCLL.sub)[rowData(dxdCLL.sub)$groupID == "ENSG00000105705",]$symbol <- "SUGP1"Are all the proteins present in the splicing dataset?
all(resTab$id %in% rowData(dxdCLL.sub)$groupID)[1] TRUEYes
How many samples in SF3B1 mutated and unmutated group?
sumTab <- colData(dxdCLL.sub) %>% data.frame() %>%
distinct(sample,.keep_all = TRUE)
table(sumTab$SF3B1)
0 1
169 31 Differential exon usage test using DEXseq
dxdCLL.sub$sample <- droplevels(dxdCLL.sub$sample)
dxdCLL.sub$condition <- dxdCLL.sub$SF3B1
formulaFullModel <- ~ sample + exon + condition:exon + IGHV:exon + trisomy12:exon
formulaReducedModel <- ~ sample + exon + IGHV:exon + trisomy12:exon
dxdCLL.sub <- estimateDispersions(dxdCLL.sub, formula = formulaFullModel)
dxdCLL.sub <- testForDEU(dxdCLL.sub, reducedModel = formulaReducedModel,
fullModel = formulaFullModel)
save(dxdCLL.sub, file = "../output/dxdCLL.RData")#load results
load("../output/dxdCLL.RData")Any significant associations?
resDxd <- DEXSeqResults(dxdCLL.sub)
resTab <- resDxd %>% data.frame() %>%
rownames_to_column("id") %>%
filter(pvalue < 0.05) %>%
mutate(symbol = rowData(dxdCLL.sub[id,])$symbol) %>%
select(symbol, featureID, groupID, pvalue, padj)
resTab symbol featureID groupID pvalue padj
1 SUGP1 E013 ENSG00000105705 9.356828e-14 1.387150e-11
2 SUGP1 E016 ENSG00000105705 4.555527e-22 1.350714e-19
3 SUGP1 E017 ENSG00000105705 4.670237e-27 2.769451e-24
4 SUGP1 E025 ENSG00000105705 3.031991e-21 5.993236e-19
5 MSH6 E001 ENSG00000116062 2.159328e-05 2.560963e-03
6 TPP2 E005 ENSG00000134900 1.300002e-02 7.709012e-01
7 TPP2 E015 ENSG00000134900 1.220560e-02 7.709012e-01
8 MICAL1 E044 ENSG00000135596 4.161688e-02 9.999451e-01
9 ATM E038 ENSG00000149311 3.217624e-02 9.999451e-01
10 SEPT2 E027 ENSG00000168385 9.304578e-03 7.248072e-01
11 SEPT2 E028 ENSG00000168385 2.991883e-03 2.956977e-01
12 SEPT2 E029 ENSG00000168385 9.778174e-03 7.248072e-01
13 NT5DC1 E016 ENSG00000178425 4.308154e-02 9.999451e-01
14 UBA7 E049 ENSG00000182179 2.047708e-02 9.999451e-01Two genes pass 10% FDR, SUGP1 and MSH6
Plot exon usage
SUGP1 (ENSG00000105705)
plotDEXSeq(resDxd, "ENSG00000105705", displayTranscripts = TRUE, legend = TRUE, norCounts = TRUE)
MSH6 (ENSG00000116062)
plotDEXSeq(resDxd, "ENSG00000116062", displayTranscripts = TRUE, legend = TRUE, norCounts = TRUE)
sessionInfo()R version 3.6.2 (2019-12-12)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /g/easybuild/x86_64/CentOS/7/haswell/software/OpenBLAS/0.3.7-GCC-8.3.0/lib/libopenblas_haswellp-r0.3.7.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] forcats_0.4.0 stringr_1.4.0
[3] dplyr_0.8.3 purrr_0.3.3
[5] readr_1.3.1 tidyr_1.0.0
[7] tibble_2.1.3 ggplot2_3.2.1
[9] tidyverse_1.3.0 DEXSeq_1.32.0
[11] RColorBrewer_1.1-2 AnnotationDbi_1.48.0
[13] DESeq2_1.26.0 SummarizedExperiment_1.16.1
[15] DelayedArray_0.12.2 matrixStats_0.55.0
[17] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0
[19] IRanges_2.20.2 S4Vectors_0.24.3
[21] Biobase_2.46.0 BiocGenerics_0.32.0
[23] BiocParallel_1.20.1 UpSetR_1.4.0
[25] proDA_1.0.0 jyluMisc_0.1.5
[27] pheatmap_1.0.12 cowplot_1.0.0
loaded via a namespace (and not attached):
[1] shinydashboard_0.7.1 tidyselect_0.2.5 RSQLite_2.1.4
[4] htmlwidgets_1.5.1 grid_3.6.2 maxstat_0.7-25
[7] munsell_0.5.0 codetools_0.2-16 statmod_1.4.32
[10] DT_0.10 withr_2.1.2 colorspace_1.4-1
[13] highr_0.8 knitr_1.26 rstudioapi_0.10
[16] ggsignif_0.6.0 labeling_0.3 git2r_0.26.1
[19] slam_0.1-46 GenomeInfoDbData_1.2.2 hwriter_1.3.2
[22] KMsurv_0.1-5 farver_2.0.1 bit64_0.9-7
[25] rprojroot_1.3-2 vctrs_0.2.0 generics_0.0.2
[28] TH.data_1.0-10 xfun_0.11 BiocFileCache_1.10.2
[31] sets_1.0-18 R6_2.4.1 locfit_1.5-9.1
[34] bitops_1.0-6 fgsea_1.12.0 assertthat_0.2.1
[37] promises_1.1.0 scales_1.1.0 multcomp_1.4-11
[40] nnet_7.3-12 gtable_0.3.0 extraDistr_1.8.11
[43] sandwich_2.5-1 workflowr_1.6.1 rlang_0.4.2
[46] zeallot_0.1.0 genefilter_1.68.0 splines_3.6.2
[49] lazyeval_0.2.2 acepack_1.4.1 broom_0.5.3
[52] checkmate_1.9.4 modelr_0.1.5 yaml_2.2.0
[55] abind_1.4-5 crosstalk_1.0.0 backports_1.1.5
[58] httpuv_1.5.2 Hmisc_4.3-0 tools_3.6.2
[61] relations_0.6-9 gplots_3.0.1.1 Rcpp_1.0.3
[64] plyr_1.8.5 base64enc_0.1-3 visNetwork_2.0.9
[67] progress_1.2.2 zlibbioc_1.32.0 RCurl_1.95-4.12
[70] prettyunits_1.0.2 ggpubr_0.2.4 rpart_4.1-15
[73] openssl_1.4.1 zoo_1.8-6 haven_2.2.0
[76] cluster_2.1.0 exactRankTests_0.8-31 fs_1.3.1
[79] magrittr_1.5 data.table_1.12.8 openxlsx_4.1.4
[82] reprex_0.3.0 survminer_0.4.6 mvtnorm_1.0-11
[85] hms_0.5.2 shinyjs_1.1 mime_0.7
[88] evaluate_0.14 xtable_1.8-4 XML_3.98-1.20
[91] rio_0.5.16 readxl_1.3.1 gridExtra_2.3
[94] compiler_3.6.2 biomaRt_2.42.0 KernSmooth_2.23-16
[97] crayon_1.3.4 htmltools_0.4.0 later_1.0.0
[100] Formula_1.2-3 geneplotter_1.64.0 lubridate_1.7.4
[103] DBI_1.1.0 dbplyr_1.4.2 rappdirs_0.3.1
[106] MASS_7.3-51.4 Matrix_1.2-18 car_3.0-5
[109] cli_2.0.0 marray_1.64.0 gdata_2.18.0
[112] igraph_1.2.4.2 pkgconfig_2.0.3 km.ci_0.5-2
[115] foreign_0.8-72 piano_2.2.0 xml2_1.2.2
[118] annotate_1.64.0 XVector_0.26.0 drc_3.0-1
[121] rvest_0.3.5 digest_0.6.23 Biostrings_2.54.0
[124] rmarkdown_2.0 cellranger_1.1.0 fastmatch_1.1-0
[127] survMisc_0.5.5 htmlTable_1.13.3 curl_4.3
[130] Rsamtools_2.2.1 shiny_1.4.0 gtools_3.8.1
[133] lifecycle_0.1.0 nlme_3.1-143 jsonlite_1.6
[136] carData_3.0-3 fansi_0.4.0 askpass_1.1
[139] limma_3.42.1 pillar_1.4.2 lattice_0.20-38
[142] fastmap_1.0.1 httr_1.4.1 plotrix_3.7-7
[145] survival_3.1-8 glue_1.3.1 zip_2.0.4
[148] bit_1.1-14 stringi_1.4.3 blob_1.2.0
[151] latticeExtra_0.6-28 caTools_1.17.1.3 memoise_1.1.0