1 Rework DIA-Umpire to feed OpenSWATH.

In the following blocks I want to use DIA Umpire to create transition libraries for openswath, then I want to run openswath and score the runs.

cd ~/scratch/proteomics/mycobacterium_tuberculosis_2018
module add openms

type="mzXML"
export VERSION="20190228"
basedir="${HOME}/scratch/proteomics/mycobacterium_tuberculosis_2018"
base_input="${basedir}/results/01${type}/dia/${VERSION}/"
umpire_inputs=$(/usr/bin/find "${base_input}" -name "*.${type}" | sort)
echo "Checking in: ${umpire_inputs}"
for input in ${umpire_inputs};
do
    in_name=$(basename $input ".${type}")
    out_name="${in_name}_Q1.mgf"
    if [[ ! -f "${base_input}/${out_name}" ]]; then
        echo "The output file: ${out_name} already exists."
    else
        java -jar DIA_Umpire_SE.jar ${input} diaumpire_se.params
    fi
done

msconvert --mzXML results/02mgf/*.mgf
mv *.mzXML results/03_dia_umpire_mzxml

comet \
    -Pparameters/comet_dia_umpire_params.txt \
    results/03_dia_umpire_mzxml/*.mzXML

xinteract \
    -dDECOY_ \
    -OARPd \
    -Ninteract.comet.pep.xml \
    results/03_dia_umpire_mzxml/*.pep.xml

mv interact.comet.pep.xml results/04_dia_umpire_xinteract

ProteinProphet \
    results/04_dia_umpire_xinteract/interact.comet.pep.xml \
    results/05_dia_umpire_prophet/combined.prot.xml

InterProphetParser \
    DECOY=DECOY \
    results/04_dia_umpire_xinteract/interact.comet.pep.xml \
    results/05_dia_umpire_prophet/iProphet.pep.xml

Mayu.pl \
    -A results/05_dia_umpire_prophet/iProphet.pep.xml \
    -C reference/mtb_irt.fasta \
    -E DECOY

mayu_output <- "../2019-05-12_12.37.03_main_1.07.csv"
number <- hpgltools::extract_mayu_pps_fdr(mayu_output)
message("The number is: ", number)
## 0.43291

## Rerunning because writing the file failed.
spectrast \
    -cNSpecLib -cICID-QTOF \
    -cf "Protein! ~ DECOY_" \
    -cP0.4237 \
    -c_IRTreference/irt.txt \
    -c_IRR results/05_dia_umpire_prophet/iProphet.pep.xml

spectrast \
    -cNSpecLib_cons \
    -cICID-QTOF \
    -cAC SpecLib.splib

spectrast2tsv.py \
    -l 350,2000 \
    -s b,y \
    -x 1,2 \
    -o 6 \
    -n 6 \
    -p 0.05 \
    -d -e \
    -k openswath \
    -w windows/2018_0817BrikenTrypsinDIA19.txt \
    -a SpecLib_cons_openswath.tsv \
    SpecLib_cons.sptxt

TargetedFileConverter \
    -in SpecLib_cons_openswath.tsv \
    -in_type tsv \
    -out SpecLib_cons_openswath.TraML \
    -out_type TraML

OpenSwathDecoyGenerator \
    -in SpecLib_cons_openswath.TraML \
    -out SpecLib_cons_openswath_decoy.TraML \
    -method shuffle
##    -exclude_similar \
##    -similarity_threshold 0.05 \
##    -identity_threshold 0.7

TargetedFileConverter \
    -in SpecLib_cons_openswath_decoy.TraML \
    -in_type TraML \
    -out SpecLib_cons_openswath_decoy.tsv \
    -out_type tsv

TargetedFileConverter \
    -in SpecLib_cons_openswath_decoy.TraML \
    -in_type TraML \
    -out SpecLib_cons_openswath_decoy.pqp \
    -out_type pqp

export VERSION=${VERSION:-20190327}
echo "Loading environment modules and parameters for version: ${VERSION}."
source "parameters/${VERSION}_settings.sh"

echo "Invoking the OpenSwathWorkflow using local comet-derived transitions."
type="diaumpire"
input_type="mzXML"
export TRANSITION_PREFIX="SpecLib_cons_openswath_decoy"
echo "Checking in, the transition library is: ${TRANSITION_PREFIX}.pqp"
base_mzxmldir="results/01${input_type}/dia/${VERSION}"
swath_inputs=$(/usr/bin/find "${base_mzxmldir}" -name *.${input_type} -print | sort)
echo "Checking in, the inputs are: ${swath_inputs}"
mkdir -p "${SWATH_OUTDIR}_${type}"
pypdir="${PYPROPHET_OUTDIR}_${type}"
mkdir -p "${pypdir}"
for input in ${swath_inputs}
do
    name=$(basename "${input}" ".${input_type}")
    echo "Starting openswath run, library type ${type} for ${name} using ${MZ_WINDOWS} windows at $(date)."
    swath_output_prefix="${SWATH_OUTDIR}_${type}/${name}_${DDA_METHOD}"
    pyprophet_output_prefix="${PYPROPHET_OUTDIR}_${type}/${name}_${DDA_METHOD}"
    echo "Deleting previous swath output file: ${swath_output_prefix}.osw"
    rm -f "${swath_output_prefix}.osw"
    rm -f "${swath_output_prefix}.tsv"
    OpenSwathWorkflow \
        -in "${input}" \
        -force \
        -sort_swath_maps \
        -min_upper_edge_dist 1 \
        -mz_correction_function "quadratic_regression_delta_ppm" \
        -Scoring:TransitionGroupPicker:background_subtraction "original" \
        -Scoring:stop_report_after_feature "5" \
        -swath_windows_file "windows/openswath_${name}.txt" \
        -tr "${TRANSITION_PREFIX}.pqp" \
        -out_tsv "${swath_output_prefix}.tsv"
    OpenSwathWorkflow \
        -in "${input}" \
        -force \
        -sort_swath_maps \
        -min_upper_edge_dist 1 \
        -mz_correction_function "quadratic_regression_delta_ppm" \
        -Scoring:TransitionGroupPicker:background_subtraction "original" \
        -Scoring:stop_report_after_feature "5" \
        -swath_windows_file "windows/openswath_${name}.txt" \
        -tr "${TRANSITION_PREFIX}.pqp" \
        -out_osw "${swath_output_prefix}.osw"
    ##2>"${swath_output_prefix}_osw.log" 1>&2
done
swath_out=$(dirname ${swath_output_prefix})
pyprophet_out="$(dirname "${pyprophet_output_prefix}")/openswath_merged.osw"
echo "Merging osw files to ${pyprophet_out}"
pyprophet merge \
          --template "${TRANSITION_PREFIX}.pqp" \
          --out="${pyprophet_out}" \
          ${swath_out}/*.osw
pyprophet score --in="${pyprophet_out}"
pyprophet export --in="${pyprophet_out}" --out "test.tsv"
## pyprophet always exports to the current working directory.
final_name="$(dirname ${pyprophet_out})/$(basename ${pyprophet_out} ".osw").tsv"
echo $final_name
mv "test.tsv"
ls -ld "${pyprophet_out}"

tric_tb="${TRIC_OUTDIR}_tuberculist"
mkdir -p "${tric_tb}"
feature_alignment.py \
    --force \
    --in "./${pypdir}/"*.tsv \
    --out "${tric_tb}/${SEARCH_METHOD}_${DDA_METHOD}.tsv" \
    --out_matrix "${tric_tb}/${DDA_METHOD}_outmatrix.tsv" \
    --out_meta "${tric_tb}/${DDA_METHOD}_meta.tsv"
2>"${tric_tb}/feature_alignment.err" \
 1>"${tric_tb}/feature_alignment.out"
echo "Wrote final output to ${tric_tb}/${SEARCH_METHOD}_${DDA_METHOD}.tsv"

2 DEP usage

Thanks to Vivek, I now am aware of DEP, which does everything I wish MSstats did. The matrix given to me by tric’s feature_alignment.py I think gives me what DEP requires, along with my annotations and sample sheet.

Let us see if this is true.

mtb_gff <- "reference/mycobacterium_tuberculosis_h37rv_2.gff.gz"
mtb_genome <- "reference/mtuberculosis_h37rv_genbank.fasta"
mtb_cds <- "reference/mtb_cds.fasta"

mtb_annotations <- sm(load_gff_annotations(mtb_gff, type="gene"))
colnames(mtb_annotations) <- gsub(pattern="\\.", replacement="", x=colnames(mtb_annotations))
mtb_annotations[["description"]] <- gsub(pattern="\\+", replacement=" ",
                                         x=mtb_annotations[["description"]])
mtb_annotations[["function"]] <- gsub(pattern="\\+", replacement=" ",
                                      x=mtb_annotations[["function"]])
rownames(mtb_annotations) <- mtb_annotations[["ID"]]

mtb_microbes <- load_microbesonline_annotations(id=83332)

## The species being downloaded is: Mycobacterium tuberculosis H37Rv

tric_data <- read.csv(
  paste0("results/tric/", ver, "/whole_8mz_dia_umpire/comet_HCD.tsv"), sep="\t")
tric_data[["ProteinName"]] <- gsub(pattern="^(.*)_.*$", replacement="\\1",
                                   x=tric_data[["ProteinName"]])
sample_annot <- extract_metadata(paste0("sample_sheets/Mtb_dia_samples_20190521.xlsx"))
rownames(sample_annot)

##  [1] "s2018_0315Briken01"           "s2018_0315Briken02"          
##  [3] "s2018_0315Briken03"           "s2018_0315Briken04"          
##  [5] "s2018_0315Briken05"           "s2018_0315Briken06"          
##  [7] "s2018_0315Briken21"           "s2018_0315Briken22"          
##  [9] "s2018_0315Briken23"           "s2018_0315Briken24"          
## [11] "s2018_0315Briken25"           "s2018_0315Briken26"          
## [13] "s2018_0502BrikenDIA01"        "s2018_0502BrikenDIA02"       
## [15] "s2018_0502BrikenDIA03"        "s2018_0502BrikenDIA04"       
## [17] "s2018_0502BrikenDIA05"        "s2018_0502BrikenDIA06"       
## [19] "s2018_0502BrikenDIA07"        "s2018_0502BrikenDIA08"       
## [21] "s2018_0502BrikenDIA09"        "s2018_0502BrikenDIA10"       
## [23] "s2018_0502BrikenDIA11"        "s2018_0502BrikenDIA12"       
## [25] "s2018_0726Briken01"           "s2018_0726Briken02"          
## [27] "s2018_0726Briken03"           "s2018_0726Briken04"          
## [29] "s2018_0726Briken05"           "s2018_0726Briken06"          
## [31] "s2018_0726Briken07"           "s2018_0726Briken08"          
## [33] "s2018_0726Briken09"           "s2018_0726Briken11"          
## [35] "s2018_0726Briken12"           "s2018_0726Briken13"          
## [37] "s2018_0726Briken14"           "s2018_0726Briken15"          
## [39] "s2018_0726Briken16"           "s2018_0726Briken17"          
## [41] "s2018_0726Briken18"           "s2018_0726Briken19"          
## [43] "s2018_0817BrikenTrypsinDIA01" "s2018_0817BrikenTrypsinDIA02"
## [45] "s2018_0817BrikenTrypsinDIA03" "s2018_0817BrikenTrypsinDIA04"
## [47] "s2018_0817BrikenTrypsinDIA05" "s2018_0817BrikenTrypsinDIA06"
## [49] "s2018_0817BrikenTrypsinDIA07" "s2018_0817BrikenTrypsinDIA08"
## [51] "s2018_0817BrikenTrypsinDIA09" "s2018_0817BrikenTrypsinDIA11"
## [53] "s2018_0817BrikenTrypsinDIA12" "s2018_0817BrikenTrypsinDIA13"
## [55] "s2018_0817BrikenTrypsinDIA14" "s2018_0817BrikenTrypsinDIA15"
## [57] "s2018_0817BrikenTrypsinDIA16" "s2018_0817BrikenTrypsinDIA17"
## [59] "s2018_0817BrikenTrypsinDIA18" "s2018_0817BrikenTrypsinDIA19"

devtools::load_all("~/scratch/git/SWATH2stats_myforked")

## Loading SWATH2stats

s2s_exp <- SWATH2stats::sample_annotation(data=tric_data,
                             sample_annotation=sample_annot, verbose=TRUE,
                             fullpeptidename_column="fullpeptidename")

## Found the same mzXML files in the annotations and data.

## results/01mzXML/dia/20190327/2018_0315Briken01.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken02.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken03.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken04.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken05.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken06.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken21.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken22.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken23.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken24.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken25.mzXML

## results/01mzXML/dia/20190327/2018_0315Briken26.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA01.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA02.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA03.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA04.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA05.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA06.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA07.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA08.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA09.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA10.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA11.mzXML

## results/01mzXML/dia/20190327/2018_0502BrikenDIA12.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken01.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken02.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken03.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken04.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken05.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken06.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken07.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken08.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken09.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken11.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken12.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken13.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken14.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken15.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken16.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken17.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken18.mzXML

## results/01mzXML/dia/20190327/2018_0726Briken19.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA01.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA02.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA03.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA04.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA05.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA06.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA07.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA08.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA09.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA11.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA12.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA13.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA14.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA15.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA16.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA17.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA18.mzXML

## results/01mzXML/dia/20190327/2018_0817BrikenTrypsinDIA19.mzXML

## 60 samples were read from the annotations.

## 224796 transitions were read from the data and merged with the annotations.

decoy_lists <- assess_decoy_rate(s2s_exp)

## Number of non-decoy peptides: 17081
## Number of decoy peptides: 1801
## Decoy rate: 0.1054

## This seems a bit high to me, yesno?
fdr_overall <- assess_fdr_overall(s2s_exp, output="Rconsole", plot=TRUE)

byrun_fdr <- assess_fdr_byrun(s2s_exp, FFT=0.7, plot=TRUE, output="Rconsole")

## The average FDR by run on assay level is 0.015

## The average FDR by run on peptide level is 0.016

## The average FDR by run on protein level is 0.001

chosen_mscore <- mscore4assayfdr(s2s_exp, FFT=0.7, fdr_target=0.02)

## Target assay FDR: 0.02

## Required overall m-score cutoff: 0.0031623
## achieving assay FDR: 0.0194

prot_score <- mscore4protfdr(s2s_exp, FFT=0.7, fdr_target=0.02)

## Target protein FDR: 0.02

## Required overall m-score cutoff: 0.01
## achieving protein FDR: 0.00115

filtered_ms <- filter_mscore(s2s_exp, chosen_mscore)

## Original dimension: 221952, new dimension: 211415, difference: 10537.

filtered_fq <- filter_mscore_freqobs(s2s_exp, 0.01, 0.8, rm.decoy=FALSE)

## Peptides need to have been quantified in more conditions than: 48 in order to pass this percentage-based threshold.

## Fraction of peptides selected: 0.00058

## Original dimension: 224796, new dimension: 680, difference: 224116.

filtered_ms_fdr <- filter_mscore_fdr(filtered_ms, FFT=0.7,
                                     overall_protein_fdr_target=prot_score,
                                     upper_overall_peptide_fdr_limit=0.05)

## Target protein FDR: 0.01

## Required overall m-score cutoff: 0.01
## achieving protein FDR: 0

## filter_mscore_fdr is filtering the data...

## finding m-score cutoff to achieve desired protein FDR in protein master list..

## finding m-score cutoff to achieve desired global peptide FDR..

## Target peptide FDR: 0.05

## Required overall m-score cutoff: 0.01

## Achieving peptide FDR: 0

## Proteins selected: 
## Total proteins selected: 2412
## Final target proteins: 2412
## Final decoy proteins: 0

## Peptides mapping to these protein entries selected:
## Total mapping peptides: 16868
## Final target peptides: 16868
## Final decoy peptides: 0

## Total peptides selected from:
## Total peptides: 16868
## Final target peptides: 16868
## Final decoy peptides: 0

## Individual run FDR quality of the peptides was not calculated
## as not every run contains a decoy.

## The decoys have been removed from the returned data.

filtered_ms_fdr_pr <- filter_proteotypic_peptides(filtered_ms_fdr)

## Number of proteins detected: 2363
## Protein identifiers: Rv0577, Rv0242c, Rv3012c, Rv2467, Rv3715c, Rv2220

## Number of proteins detected that are supported by a proteotypic peptide: 2337
## Number of proteotypic peptides detected: 16728

filtered_ms_fdr_pr_all <- filter_all_peptides(filtered_ms_fdr_pr)

## Number of proteins detected: 2337

## First 6 protein identifiers: Rv0577, Rv0242c, Rv3012c, Rv2467, Rv3715c, Rv2220

filtered_ms_fdr_pr_all_str <- filter_on_max_peptides(data=filtered_ms_fdr_pr_all, n_peptides=10)

## Before filtering:
##   Number of proteins: 2337
##   Number of peptides: 16728
## 
## Percentage of peptides removed: 25.94%
## 
## After filtering:
##   Number of proteins: 2331
##   Number of peptides: 12388

filtered_all_filters <- filter_on_min_peptides(data=filtered_ms_fdr_pr_all_str, n_peptides=3)

## Before filtering:
##   Number of proteins: 2331
##   Number of peptides: 12388
## 
## Percentage of peptides removed: 0%
## 
## After filtering:
##   Number of proteins: 2284
##   Number of peptides: 12388

matrix_prefix <- file.path("results", "swath2stats", ver)
if (!file.exists(matrix_prefix)) {
  dir.create(matrix_prefix)
}
protein_matrix_all <- write_matrix_proteins(
  s2s_exp, write.csv=TRUE,
  filename=file.path(matrix_prefix, "protein_all.csv"))

## Protein overview matrix results/swath2stats/20190327/protein_all.csv written to working folder.

dim(protein_matrix_all)

## [1] 2434   61

protein_matrix_mscore <- write_matrix_proteins(
  filtered_ms, write.csv=TRUE,
  filename=file.path(matrix_prefix, "protein_matrix_mscore.csv"))

## Protein overview matrix results/swath2stats/20190327/protein_matrix_mscore.csv written to working folder.

dim(protein_matrix_mscore)

## [1] 2412   61

peptide_matrix_mscore <- write_matrix_peptides(
  filtered_ms, write.csv=TRUE,
  filename=file.path(matrix_prefix, "peptide_matrix_mscore.csv"))

## Peptide overview matrix results/swath2stats/20190327/peptide_matrix_mscore.csv written to working folder.

dim(peptide_matrix_mscore)

## [1] 16868    61

protein_matrix_filtered <- write_matrix_proteins(
  filtered_all_filters, write.csv=TRUE,
  filename=file.path(matrix_prefix, "protein_matrix_filtered.csv"))

## Protein overview matrix results/swath2stats/20190327/protein_matrix_filtered.csv written to working folder.

dim(protein_matrix_filtered)

## [1] 2284   61

peptide_matrix_filtered <- write_matrix_peptides(
  filtered_all_filters, write.csv=TRUE,
  filename=file.path(matrix_prefix, "peptide_matrix_filtered.csv"))

## Peptide overview matrix results/swath2stats/20190327/peptide_matrix_filtered.csv written to working folder.

dim(peptide_matrix_filtered)

## [1] 144819     61

intensities <- protein_matrix_filtered
cols <- gsub(x=colnames(intensities), pattern="^.*(2018.*$)", replacement="s\\1")
cols[[1]] <- "Protein"
colnames(intensities) <- cols
rownames(intensities) <- intensities[["Protein"]]
intensities[["Protein"]] <- NULL

## Check the sample names of the intensity matrix and sample annotations.
colnames(intensities) %in% rownames(sample_annot)

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE

rownames(sample_annot) %in% colnames(intensities)

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE

##rownames(sample_annot)
##colnames(intensities)

reordered <- colnames(intensities)
metadata <- sample_annot[reordered, ]
protein_expt <- create_expt(sample_annot,
                            count_dataframe=intensities,
                            gene_info=mtb_annotations)

## Reading the sample metadata.

## Warning in `[<-.factor`(`*tmp*`, iseq, value = c("undefined",
## "undefined", : invalid factor level, NA generated

## Warning in `[<-.factor`(`*tmp*`, iseq, value = c("undefined",
## "undefined", : invalid factor level, NA generated

## The sample definitions comprises: 60 rows(samples) and 28 columns(metadata fields).

## Matched 2268 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## The final expressionset has 2284 rows and 60 columns.

protein_norm <- normalize_expt(protein_expt, transform="log2", convert="cpm",
                               filter=TRUE, norm="quant")

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 8 low-count genes (2276 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 11255 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

hpgltools::plot_pca(protein_norm)$plot

protein_sub <- subset_expt(protein_norm, subset="batch!='march2018'")

## There were 60, now there are 48 samples.

protein_sub <- subset_expt(protein_sub, subset="batch!='may2018'")

## There were 48, now there are 36 samples.

protein_norm <- normalize_expt(protein_sub, transform="log2", convert="cpm",
                               filter=TRUE, norm="quant", batch="svaseq")

## This function will replace the expt$expressionset slot with:

## log2(svaseq(cpm(quant(cbcb(data)))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Warning in normalize_expt(protein_sub, transform = "log2", convert =
## "cpm", : Quantile normalization and sva do not always play well together.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (2276 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## Step 5: doing batch correction with svaseq.

## Note to self:  If you get an error like 'x contains missing values' The data has too many 0's and needs a stronger low-count filter applied.

## Passing off to all_adjusters.

## batch_counts: Before batch/surrogate estimation, 81936 entries are x>1: 100%.

## batch_counts: Before batch/surrogate estimation, 0 entries are x==0: 0%.

## The be method chose 4 surrogate variable(s).

## Attempting svaseq estimation with 4 surrogates.

hpgltools::plot_pca(protein_norm)$plot

3 Pass the data to DEP and see what happens.

devtools::load_all("~/scratch/git/DEP")

## Loading DEP

wtf <- function (proteins_unique, columns, expdesign) {
  assertthat::assert_that(is.data.frame(proteins_unique), is.integer(columns),
                          is.data.frame(expdesign))
  if (any(!c("name", "ID") %in% colnames(proteins_unique))) {
    stop("'name' and/or 'ID' columns are not present in '",
         deparse(substitute(proteins_unique)), "'.\nRun make_unique() to obtain the required columns",
         call. = FALSE)
  }
  if (any(!c("label", "condition", "replicate") %in% colnames(expdesign))) {
    stop("'label', 'condition' and/or 'replicate' columns",
         "are not present in the experimental design", call. = FALSE)
  }
  if (any(!apply(proteins_unique[, columns], 2, is.numeric))) {
    stop("specified 'columns' should be numeric", "\nRun make_se_parse() with the appropriate columns as argument",
         call. = FALSE)
  }
  if (tibble::is.tibble(proteins_unique))
    proteins_unique <- as.data.frame(proteins_unique)
  if (tibble::is.tibble(expdesign))
    expdesign <- as.data.frame(expdesign)
  rownames(proteins_unique) <- proteins_unique$name
  raw <- proteins_unique[, columns]
  raw[raw == 0] <- NA
  raw <- log2(raw)
  expdesign <- mutate(expdesign, condition = make.names(condition))
  ## I changed the following because it didn't make sense to me.
  if (is.null(expdesign[["ID"]])) {
    expdesign <- expdesign %>%
      tidyr::unite(condition, replicate, remove=FALSE)
  }

  rownames(expdesign) <- expdesign$ID
  matched <- match(make.names(delete_prefix(expdesign$label)),
                   make.names(delete_prefix(colnames(raw))))
  if (any(is.na(matched))) {
    stop("None of the labels in the experimental design match ",
         "with column names in 'proteins_unique'", "\nRun make_se() with the correct labels in the experimental design",
         "and/or correct columns specification")
  }
  colnames(raw)[matched] <- expdesign$ID
  raw <- raw[, !is.na(colnames(raw))][rownames(expdesign)]
  row_data <- proteins_unique[, -columns]
  rownames(row_data) <- row_data$name
  se <- SummarizedExperiment(assays = as.matrix(raw), colData = expdesign,
                             rowData = row_data)
  return(se)
}

library(DEP)
library(SummarizedExperiment)

## Loading required package: GenomicRanges

## Loading required package: stats4

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:DEP':
## 
##     first, rename

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:DEP':
## 
##     collapse, desc, slice

## Loading required package: GenomeInfoDb

## Loading required package: DelayedArray

## Loading required package: matrixStats

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:DEP':
## 
##     count, rowRanges

## The following objects are masked from 'package:hpgltools':
## 
##     anyMissing, rowMedians

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## Loading required package: BiocParallel

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following objects are masked from 'package:base':
## 
##     aperm, apply

design <- pData(protein_expt)
design[["label"]] <- design[["sampleid"]]
design[["replicate"]] <- design[["sampleid"]]

my_se <- SummarizedExperiment(
  assays=exprs(protein_expt),
  colData=design,
  rowData=fData(protein_expt))

mtb_unique <- as.data.frame(exprs(protein_expt))
mtb_unique[["name"]] <- rownames(mtb_unique)
mtb_unique[["ID"]] <- rownames(mtb_unique)
intensity_columns <- 1:60
## mtb_se <- wtf(mtb_unique, intensity_columns, design)
mtb_se <- make_se(mtb_unique, intensity_columns, design)

DEP::plot_frequency(mtb_se)

dim(assay(mtb_se))

## [1] 2284   60

mtb_filt <- DEP::filter_missval(mtb_se, thr=2)
dim(assay(mtb_filt))

## [1] 848  60

DEP::plot_numbers(mtb_se)

DEP::plot_coverage(mtb_se)

mtb_norm <- DEP::normalize_vsn(mtb_se)

## vsn2: 2284 x 60 matrix (1 stratum).

## Please use 'meanSdPlot' to verify the fit.

DEP::plot_normalization(mtb_se, mtb_norm)

DEP::plot_missval(mtb_se)

DEP::plot_detect(mtb_se)

mtb_imp <- DEP::impute(mtb_norm, fun="MinProb", q=0.01)

## Loading required package: imputeLCMD

## Loading required package: tmvtnorm

## Loading required package: mvtnorm

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:S4Vectors':
## 
##     expand

## Loading required package: gmm

## Loading required package: sandwich

## Loading required package: norm

## Loading required package: pcaMethods

## 
## Attaching package: 'pcaMethods'

## The following object is masked from 'package:stats':
## 
##     loadings

## Loading required package: impute

## [1] 1.708

mtb_imp_man <- DEP::impute(mtb_norm, fun="man", shift=1.8, scale=0.3)
DEP::plot_imputation(mtb_norm, mtb_imp)

mtb_diff <- DEP::test_diff(mtb_imp, type="all")

## Tested contrasts: wt_filtrate_vs_wt_whole, wt_filtrate_vs_delta_filtrate, wt_filtrate_vs_comp_filtrate, wt_filtrate_vs_delta_whole, wt_filtrate_vs_comp_whole, wt_whole_vs_delta_filtrate, wt_whole_vs_comp_filtrate, wt_whole_vs_delta_whole, wt_whole_vs_comp_whole, delta_filtrate_vs_comp_filtrate, delta_filtrate_vs_delta_whole, delta_filtrate_vs_comp_whole, comp_filtrate_vs_delta_whole, comp_filtrate_vs_comp_whole, delta_whole_vs_comp_whole

mtb_diff <- DEP::test_diff(mtb_imp, type="manual",
                           test=c("wt_filtrate_vs_wt_whole",
                                  "delta_filtrate_vs_wt_filtrate",
                                  "comp_filtrate_vs_wt_filtrate",
                                  "wt_filtrate_vs_delta_filtrate",
                                  "wt_filtrate_vs_comp_filtrate",
                                  "wt_whole_vs_delta_whole",
                                  "wt_whole_vs_comp_whole"))

## Tested contrasts: wt_filtrate_vs_wt_whole, delta_filtrate_vs_wt_filtrate, comp_filtrate_vs_wt_filtrate, wt_filtrate_vs_delta_filtrate, wt_filtrate_vs_comp_filtrate, wt_whole_vs_delta_whole, wt_whole_vs_comp_whole

mtb_dep <- DEP::add_rejections(mtb_diff, alpha=0.05, lfc=0.6)
## mtb_pca <- DEP::plot_pca(mtb_dep)
## The PCA plotter provided by DEP has some problems.

DEP::plot_cor(mtb_dep)

DEP::plot_heatmap(mtb_dep, type="centered", kmeans=TRUE)

DEP::plot_heatmap(mtb_dep, type="contrast", kmeans=TRUE)

DEP::plot_volcano(mtb_dep, contrast="wt_whole_vs_delta_whole")

DEP::plot_volcano(mtb_dep, contrast="wt_filtrate_vs_delta_filtrate")

DEP::plot_cond(mtb_dep)

mtb_result <- DEP::get_results(mtb_dep)
plot_single(mtb_dep, proteins = c("Rv0287", "Rv0288"))

written_dep <- hpgltools::write_xls(data=mtb_result, excel="excel/dep_result.xlsx")

## Saving to: excel/dep_result.xlsx

## Note: zip::zip() is deprecated, please use zip::zipr() instead

4 Take the imputed matrix into hpgltools

imputed_mtrx <- assay(mtb_imp)
colnames(imputed_mtrx) <- gsub(x=colnames(imputed_mtrx),
                               pattern="^.*(2018.*$)",
                               replacement="s\\1")
imputed_mtrx <- 2 ^ imputed_mtrx
dim(imputed_mtrx)

## [1] 2284   60

exprset <- protein_expt[["expressionset"]]
exprset_rows <- rownames(exprs(exprset))
imputed_mtrx <- imputed_mtrx[exprset_rows, ]
exprs(exprset) <- imputed_mtrx
protein_expt[["expressionset"]] <- exprset

my_norm <- normalize_expt(protein_expt, filter=TRUE, convert="cpm", norm="quant", transform="log2")

## This function will replace the expt$expressionset slot with:

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: cbcb

## Removing 0 low-count genes (2284 remaining).

## Step 2: normalizing the data with quant.

## Using normalize.quantiles.robust due to a thread error in preprocessCore.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## Step 5: not doing batch correction.

my_pca <- hpgltools::plot_pca(my_norm)
my_pca$plot

5 Minor change to plot_missval

DEP has a neat function to plot missing values. Sadly, it does not return the actual matrix, only the plot. This is nice and all, but I need the matrix, ergo this minor change.

count_defined <- function(se) {
  df <- as.data.frame(assay(se))
  na_idx <- is.na(df)
  defined_mtrx <- ifelse(is.na(df), 0, 1)
  return(defined_mtrx)
}

def_mtrx <- count_defined(mtb_se)
summary(def_mtrx)

##  wt_filtrate_2018_0315Briken01 wt_filtrate_2018_0315Briken02
##  Min.   :0.000                 Min.   :0.000                
##  1st Qu.:0.000                 1st Qu.:0.000                
##  Median :0.000                 Median :0.000                
##  Mean   :0.355                 Mean   :0.374                
##  3rd Qu.:1.000                 3rd Qu.:1.000                
##  Max.   :1.000                 Max.   :1.000                
##  wt_filtrate_2018_0315Briken03 wt_whole_2018_0315Briken04
##  Min.   :0.000                 Min.   :0.000             
##  1st Qu.:0.000                 1st Qu.:0.000             
##  Median :0.000                 Median :1.000             
##  Mean   :0.403                 Mean   :0.731             
##  3rd Qu.:1.000                 3rd Qu.:1.000             
##  Max.   :1.000                 Max.   :1.000             
##  wt_whole_2018_0315Briken05 wt_whole_2018_0315Briken06
##  Min.   :0.000              Min.   :0.000             
##  1st Qu.:1.000              1st Qu.:1.000             
##  Median :1.000              Median :1.000             
##  Mean   :0.783              Mean   :0.757             
##  3rd Qu.:1.000              3rd Qu.:1.000             
##  Max.   :1.000              Max.   :1.000             
##  wt_whole_2018_0315Briken21 wt_whole_2018_0315Briken22
##  Min.   :0.000              Min.   :0.000             
##  1st Qu.:1.000              1st Qu.:1.000             
##  Median :1.000              Median :1.000             
##  Mean   :0.795              Mean   :0.817             
##  3rd Qu.:1.000              3rd Qu.:1.000             
##  Max.   :1.000              Max.   :1.000             
##  wt_filtrate_2018_0315Briken23 wt_filtrate_2018_0315Briken24
##  Min.   :0.00                  Min.   :0.000                
##  1st Qu.:0.00                  1st Qu.:0.000                
##  Median :0.00                  Median :0.000                
##  Mean   :0.39                  Mean   :0.412                
##  3rd Qu.:1.00                  3rd Qu.:1.000                
##  Max.   :1.00                  Max.   :1.000                
##  wt_filtrate_2018_0315Briken25 wt_whole_2018_0315Briken26
##  Min.   :0.000                 Min.   :0.000             
##  1st Qu.:0.000                 1st Qu.:1.000             
##  Median :0.000                 Median :1.000             
##  Mean   :0.362                 Mean   :0.799             
##  3rd Qu.:1.000                 3rd Qu.:1.000             
##  Max.   :1.000                 Max.   :1.000             
##  delta_filtrate_2018_0502BrikenDIA01 delta_filtrate_2018_0502BrikenDIA02
##  Min.   :0.000                       Min.   :0.00                       
##  1st Qu.:0.000                       1st Qu.:0.00                       
##  Median :0.000                       Median :0.00                       
##  Mean   :0.305                       Mean   :0.35                       
##  3rd Qu.:1.000                       3rd Qu.:1.00                       
##  Max.   :1.000                       Max.   :1.00                       
##  delta_filtrate_2018_0502BrikenDIA03 comp_filtrate_2018_0502BrikenDIA04
##  Min.   :0.000                       Min.   :0.000                     
##  1st Qu.:0.000                       1st Qu.:0.000                     
##  Median :0.000                       Median :0.000                     
##  Mean   :0.296                       Mean   :0.366                     
##  3rd Qu.:1.000                       3rd Qu.:1.000                     
##  Max.   :1.000                       Max.   :1.000                     
##  comp_filtrate_2018_0502BrikenDIA05 comp_filtrate_2018_0502BrikenDIA06
##  Min.   :0.000                      Min.   :0.000                     
##  1st Qu.:0.000                      1st Qu.:0.000                     
##  Median :0.000                      Median :0.000                     
##  Mean   :0.401                      Mean   :0.359                     
##  3rd Qu.:1.000                      3rd Qu.:1.000                     
##  Max.   :1.000                      Max.   :1.000                     
##  delta_whole_2018_0502BrikenDIA07 delta_whole_2018_0502BrikenDIA08
##  Min.   :0.000                    Min.   :0.000                   
##  1st Qu.:0.000                    1st Qu.:0.000                   
##  Median :1.000                    Median :1.000                   
##  Mean   :0.747                    Mean   :0.738                   
##  3rd Qu.:1.000                    3rd Qu.:1.000                   
##  Max.   :1.000                    Max.   :1.000                   
##  delta_whole_2018_0502BrikenDIA09 comp_whole_2018_0502BrikenDIA10
##  Min.   :0.000                    Min.   :0.0000                 
##  1st Qu.:0.000                    1st Qu.:0.0000                 
##  Median :1.000                    Median :0.0000                 
##  Mean   :0.535                    Mean   :0.0814                 
##  3rd Qu.:1.000                    3rd Qu.:0.0000                 
##  Max.   :1.000                    Max.   :1.0000                 
##  comp_whole_2018_0502BrikenDIA11 comp_whole_2018_0502BrikenDIA12
##  Min.   :0.00                    Min.   :0.000                  
##  1st Qu.:1.00                    1st Qu.:0.000                  
##  Median :1.00                    Median :1.000                  
##  Mean   :0.75                    Mean   :0.679                  
##  3rd Qu.:1.00                    3rd Qu.:1.000                  
##  Max.   :1.00                    Max.   :1.000                  
##  delta_filtrate_2018_0726Briken01 delta_filtrate_2018_0726Briken02
##  Min.   :0.000                    Min.   :0.000                   
##  1st Qu.:0.000                    1st Qu.:0.000                   
##  Median :0.000                    Median :0.000                   
##  Mean   :0.063                    Mean   :0.102                   
##  3rd Qu.:0.000                    3rd Qu.:0.000                   
##  Max.   :1.000                    Max.   :1.000                   
##  delta_filtrate_2018_0726Briken03 comp_filtrate_2018_0726Briken04
##  Min.   :0.0000                   Min.   :0.0000                 
##  1st Qu.:0.0000                   1st Qu.:0.0000                 
##  Median :0.0000                   Median :0.0000                 
##  Mean   :0.0517                   Mean   :0.0517                 
##  3rd Qu.:0.0000                   3rd Qu.:0.0000                 
##  Max.   :1.0000                   Max.   :1.0000                 
##  comp_filtrate_2018_0726Briken05 comp_filtrate_2018_0726Briken06
##  Min.   :0.0000                  Min.   :0.000                  
##  1st Qu.:0.0000                  1st Qu.:0.000                  
##  Median :0.0000                  Median :0.000                  
##  Mean   :0.0823                  Mean   :0.113                  
##  3rd Qu.:0.0000                  3rd Qu.:0.000                  
##  Max.   :1.0000                  Max.   :1.000                  
##  wt_filtrate_2018_0726Briken07 wt_filtrate_2018_0726Briken08
##  Min.   :0.000                 Min.   :0.000                
##  1st Qu.:0.000                 1st Qu.:0.000                
##  Median :0.000                 Median :0.000                
##  Mean   :0.115                 Mean   :0.109                
##  3rd Qu.:0.000                 3rd Qu.:0.000                
##  Max.   :1.000                 Max.   :1.000                
##  wt_filtrate_2018_0726Briken09 delta_whole_2018_0726Briken11
##  Min.   :0.00                  Min.   :0.000                
##  1st Qu.:0.00                  1st Qu.:0.000                
##  Median :0.00                  Median :0.000                
##  Mean   :0.14                  Mean   :0.255                
##  3rd Qu.:0.00                  3rd Qu.:1.000                
##  Max.   :1.00                  Max.   :1.000                
##  delta_whole_2018_0726Briken12 delta_whole_2018_0726Briken13
##  Min.   :0.000                 Min.   :0.00                 
##  1st Qu.:0.000                 1st Qu.:0.00                 
##  Median :0.000                 Median :0.00                 
##  Mean   :0.304                 Mean   :0.23                 
##  3rd Qu.:1.000                 3rd Qu.:0.00                 
##  Max.   :1.000                 Max.   :1.00                 
##  comp_whole_2018_0726Briken14 comp_whole_2018_0726Briken15
##  Min.   :0.000                Min.   :0.000               
##  1st Qu.:0.000                1st Qu.:0.000               
##  Median :0.000                Median :0.000               
##  Mean   :0.226                Mean   :0.255               
##  3rd Qu.:0.000                3rd Qu.:1.000               
##  Max.   :1.000                Max.   :1.000               
##  comp_whole_2018_0726Briken16 wt_whole_2018_0726Briken17
##  Min.   :0.000                Min.   :0.000             
##  1st Qu.:0.000                1st Qu.:0.000             
##  Median :0.000                Median :0.000             
##  Mean   :0.288                Mean   :0.262             
##  3rd Qu.:1.000                3rd Qu.:1.000             
##  Max.   :1.000                Max.   :1.000             
##  wt_whole_2018_0726Briken18 wt_whole_2018_0726Briken19
##  Min.   :0.00               Min.   :0.000             
##  1st Qu.:0.00               1st Qu.:0.000             
##  Median :0.00               Median :0.000             
##  Mean   :0.35               Mean   :0.186             
##  3rd Qu.:1.00               3rd Qu.:0.000             
##  Max.   :1.00               Max.   :1.000             
##  delta_filtrate_2018_0817BrikenTrypsinDIA01
##  Min.   :0.000                             
##  1st Qu.:0.000                             
##  Median :0.000                             
##  Mean   :0.246                             
##  3rd Qu.:0.000                             
##  Max.   :1.000                             
##  delta_filtrate_2018_0817BrikenTrypsinDIA02
##  Min.   :0.000                             
##  1st Qu.:0.000                             
##  Median :0.000                             
##  Mean   :0.254                             
##  3rd Qu.:1.000                             
##  Max.   :1.000                             
##  delta_filtrate_2018_0817BrikenTrypsinDIA03
##  Min.   :0.000                             
##  1st Qu.:0.000                             
##  Median :0.000                             
##  Mean   :0.251                             
##  3rd Qu.:1.000                             
##  Max.   :1.000                             
##  comp_filtrate_2018_0817BrikenTrypsinDIA04
##  Min.   :0.000                            
##  1st Qu.:0.000                            
##  Median :0.000                            
##  Mean   :0.311                            
##  3rd Qu.:1.000                            
##  Max.   :1.000                            
##  comp_filtrate_2018_0817BrikenTrypsinDIA05
##  Min.   :0.000                            
##  1st Qu.:0.000                            
##  Median :0.000                            
##  Mean   :0.283                            
##  3rd Qu.:1.000                            
##  Max.   :1.000                            
##  comp_filtrate_2018_0817BrikenTrypsinDIA06
##  Min.   :0.00                             
##  1st Qu.:0.00                             
##  Median :0.00                             
##  Mean   :0.29                             
##  3rd Qu.:1.00                             
##  Max.   :1.00                             
##  wt_filtrate_2018_0817BrikenTrypsinDIA07
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :0.000                          
##  Mean   :0.201                          
##  3rd Qu.:0.000                          
##  Max.   :1.000                          
##  wt_filtrate_2018_0817BrikenTrypsinDIA08
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :0.000                          
##  Mean   :0.222                          
##  3rd Qu.:0.000                          
##  Max.   :1.000                          
##  wt_filtrate_2018_0817BrikenTrypsinDIA09
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :0.000                          
##  Mean   :0.236                          
##  3rd Qu.:0.000                          
##  Max.   :1.000                          
##  delta_whole_2018_0817BrikenTrypsinDIA11
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :1.000                          
##  Mean   :0.656                          
##  3rd Qu.:1.000                          
##  Max.   :1.000                          
##  delta_whole_2018_0817BrikenTrypsinDIA12
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :1.000                          
##  Mean   :0.655                          
##  3rd Qu.:1.000                          
##  Max.   :1.000                          
##  delta_whole_2018_0817BrikenTrypsinDIA13
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :1.000                          
##  Mean   :0.697                          
##  3rd Qu.:1.000                          
##  Max.   :1.000                          
##  comp_whole_2018_0817BrikenTrypsinDIA14
##  Min.   :0.000                         
##  1st Qu.:0.000                         
##  Median :1.000                         
##  Mean   :0.671                         
##  3rd Qu.:1.000                         
##  Max.   :1.000                         
##  comp_whole_2018_0817BrikenTrypsinDIA15
##  Min.   :0.000                         
##  1st Qu.:0.000                         
##  Median :1.000                         
##  Mean   :0.723                         
##  3rd Qu.:1.000                         
##  Max.   :1.000                         
##  comp_whole_2018_0817BrikenTrypsinDIA16
##  Min.   :0.000                         
##  1st Qu.:0.000                         
##  Median :1.000                         
##  Mean   :0.719                         
##  3rd Qu.:1.000                         
##  Max.   :1.000                         
##  wt_filtrate_2018_0817BrikenTrypsinDIA17
##  Min.   :0.000                          
##  1st Qu.:0.000                          
##  Median :0.000                          
##  Mean   :0.499                          
##  3rd Qu.:1.000                          
##  Max.   :1.000                          
##  wt_whole_2018_0817BrikenTrypsinDIA18 wt_whole_2018_0817BrikenTrypsinDIA19
##  Min.   :0.00                         Min.   :0.000                       
##  1st Qu.:0.00                         1st Qu.:0.000                       
##  Median :1.00                         Median :1.000                       
##  Mean   :0.65                         Mean   :0.654                       
##  3rd Qu.:1.00                         3rd Qu.:1.000                       
##  Max.   :1.00                         Max.   :1.000

defined_by_protein <- rowSums(def_mtrx)
head(defined_by_protein)

## iRT_Kit  Rv0001  Rv0002  Rv0003  Rv0005  Rv0006 
##      35      22      57      16      50      57

defined_by_sample <- colSums(def_mtrx)
defined_by_sample

##              wt_filtrate_2018_0315Briken01 
##                                        810 
##              wt_filtrate_2018_0315Briken02 
##                                        855 
##              wt_filtrate_2018_0315Briken03 
##                                        920 
##                 wt_whole_2018_0315Briken04 
##                                       1669 
##                 wt_whole_2018_0315Briken05 
##                                       1789 
##                 wt_whole_2018_0315Briken06 
##                                       1728 
##                 wt_whole_2018_0315Briken21 
##                                       1817 
##                 wt_whole_2018_0315Briken22 
##                                       1867 
##              wt_filtrate_2018_0315Briken23 
##                                        891 
##              wt_filtrate_2018_0315Briken24 
##                                        941 
##              wt_filtrate_2018_0315Briken25 
##                                        828 
##                 wt_whole_2018_0315Briken26 
##                                       1825 
##        delta_filtrate_2018_0502BrikenDIA01 
##                                        697 
##        delta_filtrate_2018_0502BrikenDIA02 
##                                        799 
##        delta_filtrate_2018_0502BrikenDIA03 
##                                        677 
##         comp_filtrate_2018_0502BrikenDIA04 
##                                        835 
##         comp_filtrate_2018_0502BrikenDIA05 
##                                        915 
##         comp_filtrate_2018_0502BrikenDIA06 
##                                        819 
##           delta_whole_2018_0502BrikenDIA07 
##                                       1706 
##           delta_whole_2018_0502BrikenDIA08 
##                                       1685 
##           delta_whole_2018_0502BrikenDIA09 
##                                       1223 
##            comp_whole_2018_0502BrikenDIA10 
##                                        186 
##            comp_whole_2018_0502BrikenDIA11 
##                                       1714 
##            comp_whole_2018_0502BrikenDIA12 
##                                       1552 
##           delta_filtrate_2018_0726Briken01 
##                                        144 
##           delta_filtrate_2018_0726Briken02 
##                                        233 
##           delta_filtrate_2018_0726Briken03 
##                                        118 
##            comp_filtrate_2018_0726Briken04 
##                                        118 
##            comp_filtrate_2018_0726Briken05 
##                                        188 
##            comp_filtrate_2018_0726Briken06 
##                                        259 
##              wt_filtrate_2018_0726Briken07 
##                                        263 
##              wt_filtrate_2018_0726Briken08 
##                                        248 
##              wt_filtrate_2018_0726Briken09 
##                                        319 
##              delta_whole_2018_0726Briken11 
##                                        583 
##              delta_whole_2018_0726Briken12 
##                                        694 
##              delta_whole_2018_0726Briken13 
##                                        526 
##               comp_whole_2018_0726Briken14 
##                                        516 
##               comp_whole_2018_0726Briken15 
##                                        582 
##               comp_whole_2018_0726Briken16 
##                                        658 
##                 wt_whole_2018_0726Briken17 
##                                        598 
##                 wt_whole_2018_0726Briken18 
##                                        800 
##                 wt_whole_2018_0726Briken19 
##                                        426 
## delta_filtrate_2018_0817BrikenTrypsinDIA01 
##                                        561 
## delta_filtrate_2018_0817BrikenTrypsinDIA02 
##                                        580 
## delta_filtrate_2018_0817BrikenTrypsinDIA03 
##                                        574 
##  comp_filtrate_2018_0817BrikenTrypsinDIA04 
##                                        711 
##  comp_filtrate_2018_0817BrikenTrypsinDIA05 
##                                        646 
##  comp_filtrate_2018_0817BrikenTrypsinDIA06 
##                                        662 
##    wt_filtrate_2018_0817BrikenTrypsinDIA07 
##                                        459 
##    wt_filtrate_2018_0817BrikenTrypsinDIA08 
##                                        508 
##    wt_filtrate_2018_0817BrikenTrypsinDIA09 
##                                        538 
##    delta_whole_2018_0817BrikenTrypsinDIA11 
##                                       1498 
##    delta_whole_2018_0817BrikenTrypsinDIA12 
##                                       1495 
##    delta_whole_2018_0817BrikenTrypsinDIA13 
##                                       1592 
##     comp_whole_2018_0817BrikenTrypsinDIA14 
##                                       1532 
##     comp_whole_2018_0817BrikenTrypsinDIA15 
##                                       1652 
##     comp_whole_2018_0817BrikenTrypsinDIA16 
##                                       1642 
##    wt_filtrate_2018_0817BrikenTrypsinDIA17 
##                                       1139 
##       wt_whole_2018_0817BrikenTrypsinDIA18 
##                                       1485 
##       wt_whole_2018_0817BrikenTrypsinDIA19 
##                                       1494

if (!isTRUE(get0("skip_load"))) {
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
  message(paste0("Saving to ", this_save))
  tmp <- sm(saveme(filename=this_save))
  pander::pander(sessionInfo())
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 0abc58e173be7300595d30d407b7efd4e4a512d6

## This is hpgltools commit: Thu May 9 14:56:34 2019 -0400: 0abc58e173be7300595d30d407b7efd4e4a512d6

## Saving to dia_umpire_20190308-v20190327.rda.xz

R version 3.5.3 (2019-03-11)

Platform: x86_64-pc-linux-gnu (64-bit)

locale: LC_CTYPE=en_US.UTF-8, LC_NUMERIC=C, LC_TIME=en_US.UTF-8, LC_COLLATE=en_US.UTF-8, LC_MONETARY=en_US.UTF-8, LC_MESSAGES=en_US.UTF-8, LC_PAPER=en_US.UTF-8, LC_NAME=C, LC_ADDRESS=C, LC_TELEPHONE=C, LC_MEASUREMENT=en_US.UTF-8 and LC_IDENTIFICATION=C

attached base packages: stats4, parallel, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: imputeLCMD(v.2.0), impute(v.1.56.0), pcaMethods(v.1.74.0), norm(v.1.0-9.5), tmvtnorm(v.1.4-10), gmm(v.1.6-2), sandwich(v.2.5-1), Matrix(v.1.2-17), mvtnorm(v.1.0-10), SummarizedExperiment(v.1.12.0), DelayedArray(v.0.8.0), BiocParallel(v.1.16.6), matrixStats(v.0.54.0), GenomicRanges(v.1.34.0), GenomeInfoDb(v.1.18.2), IRanges(v.2.16.0), S4Vectors(v.0.20.1), DEP(v.1.5.3), ruv(v.0.9.7), SWATH2stats(v.1.13.5), testthat(v.2.1.1), hpgltools(v.1.0), Biobase(v.2.42.0) and BiocGenerics(v.0.28.0)

loaded via a namespace (and not attached): shinydashboard(v.0.7.1), tidyselect(v.0.2.5), lme4(v.1.1-21), htmlwidgets(v.1.3), RSQLite(v.2.1.1), AnnotationDbi(v.1.44.0), grid(v.3.5.3), devtools(v.2.0.2), munsell(v.0.5.0), codetools(v.0.2-16), preprocessCore(v.1.44.0), DT(v.0.6), withr(v.2.1.2), colorspace(v.1.4-1), GOSemSim(v.2.8.0), knitr(v.1.23), rstudioapi(v.0.10), DOSE(v.3.8.2), mzID(v.1.20.1), labeling(v.0.3), urltools(v.1.7.3), GenomeInfoDbData(v.1.2.0), polyclip(v.1.10-0), bit64(v.0.9-7), farver(v.1.1.0), rprojroot(v.1.3-2), xfun(v.0.7), R6(v.2.4.0), doParallel(v.1.0.14), locfit(v.1.5-9.1), bitops(v.1.0-6), fgsea(v.1.8.0), gridGraphics(v.0.4-1), assertthat(v.0.2.1), promises(v.1.0.1), scales(v.1.0.0), ggraph(v.1.0.2), enrichplot(v.1.2.0), gtable(v.0.3.0), affy(v.1.60.0), sva(v.3.30.1), processx(v.3.3.1), rlang(v.0.3.4), genefilter(v.1.64.0), mzR(v.2.16.2), GlobalOptions(v.0.1.0), splines(v.3.5.3), rtracklayer(v.1.42.2), lazyeval(v.0.2.2), selectr(v.0.4-1), europepmc(v.0.3), BiocManager(v.1.30.4), yaml(v.2.2.0), reshape2(v.1.4.3), GenomicFeatures(v.1.34.8), backports(v.1.1.4), httpuv(v.1.5.1), qvalue(v.2.14.1), clusterProfiler(v.3.10.1), tools(v.3.5.3), usethis(v.1.5.0), ggplotify(v.0.0.3), ggplot2(v.3.1.1), affyio(v.1.52.0), gplots(v.3.0.1.1), RColorBrewer(v.1.1-2), sessioninfo(v.1.1.1), MSnbase(v.2.8.3), ggridges(v.0.5.1), Rcpp(v.1.0.1), plyr(v.1.8.4), base64enc(v.0.1-3), progress(v.1.2.2), zlibbioc(v.1.28.0), purrr(v.0.3.2), RCurl(v.1.95-4.12), ps(v.1.3.0), prettyunits(v.1.0.2), GetoptLong(v.0.1.7), viridis(v.0.5.1), cowplot(v.0.9.4), zoo(v.1.8-5), cluster(v.2.0.9), ggrepel(v.0.8.1), colorRamps(v.2.3), fs(v.1.3.1), variancePartition(v.1.12.3), magrittr(v.1.5), data.table(v.1.12.2), DO.db(v.2.9), openxlsx(v.4.1.0), circlize(v.0.4.6), triebeard(v.0.3.0), packrat(v.0.5.0), ProtGenerics(v.1.14.0), pkgload(v.1.0.2), mime(v.0.6), hms(v.0.4.2), evaluate(v.0.13), xtable(v.1.8-4), pbkrtest(v.0.4-7), XML(v.3.98-1.19), shape(v.1.4.4), gridExtra(v.2.3), compiler(v.3.5.3), biomaRt(v.2.38.0), tibble(v.2.1.1), KernSmooth(v.2.23-15), ncdf4(v.1.16.1), crayon(v.1.3.4), minqa(v.1.2.4), htmltools(v.0.3.6), later(v.0.8.0), mgcv(v.1.8-28), corpcor(v.1.6.9), tidyr(v.0.8.3), DBI(v.1.0.0), tweenr(v.1.0.1), ComplexHeatmap(v.1.20.0), MASS(v.7.3-51.4), boot(v.1.3-22), readr(v.1.3.1), cli(v.1.1.0), vsn(v.3.50.0), gdata(v.2.18.0), igraph(v.1.2.4.1), pkgconfig(v.2.0.2), rvcheck(v.0.1.3), GenomicAlignments(v.1.18.1), MALDIquant(v.1.19.3), xml2(v.1.2.0), foreach(v.1.4.4), annotate(v.1.60.1), XVector(v.0.22.0), rvest(v.0.3.4), stringr(v.1.4.0), callr(v.3.2.0), digest(v.0.6.19), Biostrings(v.2.50.2), rmarkdown(v.1.12), fastmatch(v.1.1-0), edgeR(v.3.24.3), curl(v.3.3), shiny(v.1.3.2), Rsamtools(v.1.34.1), gtools(v.3.8.1), rjson(v.0.2.20), nloptr(v.1.2.1), nlme(v.3.1-140), jsonlite(v.1.6), desc(v.1.2.0), viridisLite(v.0.3.0), limma(v.3.38.3), pillar(v.1.4.0), lattice(v.0.20-38), httr(v.1.4.0), pkgbuild(v.1.0.3), survival(v.2.44-1.1), GO.db(v.3.7.0), glue(v.1.3.1), remotes(v.2.0.4), fdrtool(v.1.2.15), zip(v.2.0.2), UpSetR(v.1.3.3), iterators(v.1.0.10), pander(v.0.6.3), bit(v.1.1-14), ggforce(v.0.2.2), stringi(v.1.4.3), blob(v.1.1.1), caTools(v.1.17.1.2), memoise(v.1.1.0) and dplyr(v.0.8.1)

tmp <- loadme(filename=this_save)