bp_utils/CompleteCaseNormalize.R at master · kroganlab/bp_utils · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120


# peptideDF is expected to be a table ready to pass to MSstats::dataProcess
# it will be modified in place using data.table functions

CompleteCaseNormalize.inPlace <- function(peptideDF, minIntensityQuantile = 0){
  if (!"data.table" %in% class (peptideDF)){
    setDT(peptideDF)
    data.frame.only <- TRUE
  } else{
    data.frame.only <- FALSE
  }

  minIntensity <- quantile(peptideDF$Intensity, minIntensityQuantile)

  #make a dt of complete cases, by passing through a wide format
  pepDF.wide <- dcast(peptideDF[Intensity > minIntensity & !is.na(Intensity)],
                      ProteinName+PeptideSequence+PrecursorCharge~BioReplicate+Condition+Run, value.var="Intensity")
  completeFeatures <- pepDF.wide[complete.cases(pepDF.wide), .(ProteinName, PeptideSequence, PrecursorCharge)]
  message ("Normalizing based on intensities in ", nrow(completeFeatures), " complete features (having intensity > ", minIntensity, " ")
  setkey (peptideDF, ProteinName, PeptideSequence, PrecursorCharge)
  completeLong <- peptideDF[completeFeatures]

  overallMedian <- median (log2(completeLong$Intensity))
  adjustments <- completeLong[,.(adjustment = overallMedian - median(log2(Intensity))),by = .(Condition,  BioReplicate)]
  print(adjustments)

  setkey(peptideDF, Condition, BioReplicate)

  #normalization:

  peptideDF[adjustments, Intensity := Intensity * 2^adjustment]

  if (data.frame.only)
    setDF(peptideDF)

  invisible(peptideDF)
}


# as above, but for evidence files
CompleteCaseNormalize.MaxQuant.inPlace <- function(peptideDF, minIntensityQuantile = 0){
  if (!"data.table" %in% class (peptideDF)){
    setDT(peptideDF)
    data.frame.only <- TRUE
  } else{
    data.frame.only <- FALSE
  }

  minIntensity <- quantile(peptideDF$Intensity, minIntensityQuantile, na.rm=TRUE)

  #make a dt of complete cases, by passing through a wide format
  pepDF.wide <- dcast(peptideDF[Proteins != "" & !grepl("CON|REV", Proteins) & Intensity > minIntensity & !is.na(Intensity)],
                      Proteins+`Modified sequence`+Charge~`Raw file`, value.var="Intensity", fun.aggregate = sum)


  completeFeatures <- pepDF.wide[,.(anyZero = any (unlist(.SD) == 0)), by = .(Proteins, `Modified sequence`, Charge)
                                 ][anyZero == FALSE, .(Proteins, `Modified sequence`, Charge)]
  message ("Normalizing based on intensities in ", nrow(completeFeatures), " complete features (having intensity > ", minIntensity, " ")
  setkey (peptideDF, Proteins, `Modified sequence`, Charge)
  completeLong <- peptideDF[!is.na(Intensity)][completeFeatures]

  overallMedian <- median (log2(completeLong$Intensity))
  adjustments <- completeLong[,.(adjustment = overallMedian - median(log2(Intensity))),by = .(`Raw file`)]
  print(adjustments)

  setkey(peptideDF, `Raw file`)

  #normalization:

  peptideDF[adjustments, Intensity := Intensity * 2^adjustment]

  if (data.frame.only)
    setDF(peptideDF)

  invisible(peptideDF)
}


NearlyCompleteCaseNormalize.inPlace <- function(peptideDF, require = 0.50, minIntensityQuantile = 0){
  if (!"data.table" %in% class (peptideDF)){
    setDT(peptideDF)
    data.frame.only <- TRUE
  } else{
    data.frame.only <- FALSE
  }

  minIntensity <- quantile(peptideDF$Intensity, minIntensityQuantile)


  # make a dt of complete cases, by passing through a wide format
  totalRuns <- nrow(unique(peptideDF[,.(BioReplicate,Condition,Run)]))
  minRunsReq <- ceiling(totalRuns * require)
  passingFeatures <- peptideDF[!is.na(Intensity) & Intensity > minIntensity,
                               .N,
                               by = .(ProteinName,PeptideSequence,PrecursorCharge)
                               ][N >= minRunsReq,.(ProteinName,PeptideSequence,PrecursorCharge) ]

  message (nrow(passingFeatures), " features present in at least ", minRunsReq, " runs for normalization (having intensity > ", minIntensity, " ")
  completeLong <- peptideDF[passingFeatures, , on = .(ProteinName, PeptideSequence, PrecursorCharge)]

  overallMedian <- median (log2(completeLong$Intensity), na.rm=TRUE)
  adjustments <- completeLong[,.(adjustment = overallMedian - median(log2(Intensity), na.rm=TRUE)),by = .(Condition,  BioReplicate)]
  print(adjustments)

  setkey(peptideDF, Condition, BioReplicate)

  #normalization:

  peptideDF[adjustments, Intensity := Intensity * 2^adjustment]

  if (data.frame.only)
    setDF(peptideDF)

  invisible(peptideDF)


}