-
Notifications
You must be signed in to change notification settings - Fork 6
/
XCMS_master.Rmd
203 lines (164 loc) · 7.37 KB
/
XCMS_master.Rmd
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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
---
title: "XCMS_master"
author: "Katherine Heal"
date: "November 9, 2016"
output: html_document
---
This code is just the XCMS and CAMERA functions. It should only be run on a high powered computer.
#Load Libraries and source code
```{r Libraries, results='hide', message=FALSE, warning=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(plyr)
library(ggplot2)
library(gridExtra)
library(dplyr)
library(seqinr)
library(lubridate)
library(reshape2)
library(tidyr)
library(Hmisc)
library(gtools)
library(cowplot)
require(RColorBrewer)
library(readr)
library(plotly)
library(stringr)
library(xcms)
library(IPO)
```
#Setup Experiment
This experiment was processed through XCMS on 3/2/17 on personal computer
```{r Set experiment and fraction name. Raw data file -> "rawdat" -> Fraction (i.e. "CyanoDCM") -> Experiment ("Eddy_samples"), results='hide', message=FALSE, warning=FALSE}
Experiment <- "Gradients2016"
Location <- "/Users/katherine2/Google_Drive/00_XCMS_Working/"
FractionList <- c("HILICNeg", "HILICPos", "CyanoAq", "CyanoDCM")
ExpDIR <-paste(Location, Experiment, sep = "", collapse = NULL)
```
If you have not set up your results directories yet, run this. Once you have your xcms run, you don't need to rerun this. As long as the final line is commented out though, its okay.
```{r Set your Directories, grab Params and ISList, results='hide', message=FALSE, warning=FALSE}
FractionDIR <- c()
ResultsDIR <- c()
DataDIR <- c()
for (j in 1:length(FractionList)){
FractionDIR[j] <- paste(ExpDIR, "/", FractionList[j], sep = "", collapse = NULL)
ResultsDIR[j] <- paste(FractionDIR[j], "/Results", Sys.Date(), sep = "", collapse = NULL)
#dir.create(file.path(ResultsDIR[j]), showWarnings = FALSE) #Comment out if you don't need to make the folder i.e. you've already run XCMS and need to redo the next steps
DataDIR <- paste(FractionDIR, "/rawdat", sep = "", collapse = NULL)
}
Dirs <- data.frame(FractionList, FractionDIR, ResultsDIR, DataDIR)
#Modify this to get to your results of choice
Dirs <- Dirs %>% mutate(ResultsDIR = ResultsDIR %>%
str_replace("Results2017-05-11",
"Results2017-05-09"))
rownames(Dirs) <- Dirs$FractionList
#write.csv(Dirs, paste(ExpDIR, "/Dirs.csv", sep="", collapse =NULL))#Write this csv when you are satisfied with it - and then comment it out so you don't accidently rewrite your Dirs CSV.
```
Run this if you have already made your directories (to access older results)
```{r, results='hide', message=FALSE, warning=FALSE}
#Get your directories with preferred results
Dirs <- read.csv(paste(ExpDIR, "/Dirs.csv", sep="", collapse =NULL), row.names=1) #If this doesn't work, make sure last line of chunk B4 isn't commented out
#Get your functions
source(paste(Location, "XCMS_funs.R", sep= "", collapse = NULL))
#Get your parameters
Params <- read.csv(paste(Location, "Params.csv", sep= "", collapse = NULL), row.names=1)
Params <- Params[ , 1:4]
```
#Do XCMS, get peak Lists, Print out TICs
Running with PEAKGROUPS on all fractions on 5/10/2017
```{r perform xcms on all fractions, results='hide', message=FALSE, warning=FALSE}
#for (j in 1:length(FractionList)){
j=4
Fraction <- FractionList[j]
DataDIR <- as.character(Dirs[Fraction, "DataDIR"])
DatFiles <- list.files(DataDIR, recursive = TRUE, full.names = TRUE)
DatFiles <- DatFiles[grepl(".mzXML", DatFiles)]
#Save your parameters used for your xcms-ing
write.csv(Params[,1:4], paste(Dirs[Fraction, "ResultsDIR"], sep = "","/Params_used.csv", collapse = NULL), row.names = TRUE)
#Do the XCMS
setwd(DataDIR)
xset1 <- doPeakPick(DatFiles)
xset1 <- dogrouping(c(xset1))
setwd(as.character(Dirs[Fraction, "ResultsDIR"]))
save(xset1, file="xset1_initial.RData")
#Do RT Correction (of some iterations)
for (k in 1:Params["RTITs", Fraction]){
xset1 <- doRetcorPG(xset1) #Change to doRetcor for obiwarp, doRetcorPG for peak groups (better for cyano)
xset1 <- dogrouping(xset1)
}
save(xset1, file="xset2_RTCorrected.RData")
#RT Correction plot
png("RTCorrplot.png", width = 4, height = 5, units = "in", res = 300)
plotrt(xset1, leg = F, densplit = T)
dev.off()
#Write unfilled dataframe
Peaks.unfilled <- peakTable(xset1)
Peaks.unfilled$MassFeature <- paste("I", round((Peaks.unfilled$mz),digits=4), "R", round( Peaks.unfilled$rt/60, digits=2), sep="")
Peaks.unfilled$groupname <- groupnames(xset1)
Peaks.unfilled$RT <- Peaks.unfilled$rt/60
Peaks.unfilled$Fraction <- Fraction
write.csv(Peaks.unfilled, "xset.unfilled.csv")
#Fill in peaks, write dataframe
xset3 <- fillPeaks(xset1)
save(xset3, file = "xset3_RTCorrected_PeaksFilled.RData")
Peaks.filled <- peakTable(xset3)
Peaks.filled$MassFeature <- paste("I", round((Peaks.filled$mz),digits=4), "R", round( Peaks.filled$rt/60, digits=2), sep="")
Peaks.filled$groupname <- groupnames(xset3)
Peaks.filled$RT <- Peaks.filled$rt/60
Peaks.filled$Fraction <- Fraction
write.csv(Peaks.filled, "xset.allpeaks.csv")
getTICs(xcmsSet=xset3, pdfname="TICs.pdf",rt="corrected")
PPS <- as.data.frame(calcPPS(xset3))
write.csv(PPS, "PPS_results.csv")
#}
```
#Run Camera_Function and run the camera post processing ----
Needs Camera_Fuction.R , CAMERApostprocessing.fuction.R, XCMS results
This was run on on peak group retention time correlated data on 5/10/17
```{r Calculate adducts and isotopes for mass features and calculate neutral mass for ions with two or more adducts, results='hide', message=FALSE, warning=FALSE'}
source(paste(Location, "Camera_Function.R", sep=""))
source(paste(Location, "CAMERApostprocessing.function.R", sep=""))
#for (j in 1:length(FractionList)){
j=4
Fraction <- FractionList[j]
load(paste(as.character(Dirs[Fraction, "ResultsDIR"]), "/xset3_RTCorrected_PeaksFilled.RData", sep=""))
if (Fraction == "HILICNeg") {Pol = "negative"} else {Pol = "positive"}
CAMERAlist <- camera (xset3, Polarity = Pol, PPM = 5)
save(CAMERAlist, file = paste(as.character(Dirs[Fraction, "ResultsDIR"]), "/CAMERAlist.RData", sep=""))
xset.annot <- CAMERAlist[[1]]
xsaFA <- CAMERAlist[[2]]
IonList <- camerapostprocess(xsaFA, xset.annot)
save(IonList, file = paste(as.character(Dirs[Fraction, "ResultsDIR"]), "/IonList.RData", sep=""))
OtherIons <- IonList[[3]]
#}
```
#Do Very basic Filtering on Data
Needs XCMS only
ReRun on 5/31/17, make filter more stark
```{r Making Averages and basic filtering}
for (j in 1:4){
Fraction <- FractionList[j]
setwd(as.character(Dirs[Fraction, "ResultsDIR"]))
Peaks.filled <- read.csv("xset.allpeaks.csv")
#Get data file names
DatNames <- colnames(Peaks.filled)[grepl("X17", colnames(Peaks.filled))]
#Specifiy which runs are Blks, samples, and treatments
Blks <- DatNames[grepl("_Blk_", DatNames)]
Samps <- DatNames[grepl("_Smp_", DatNames)]
Poos <- DatNames[grepl("_Poo_", DatNames)]
#Set rt max and min
rtminparam <- Params["RTMIN", Fraction]
rtmaxparam <- Params["RTMAX", Fraction]
if (length(Blks) > 1) {Peaks.filled$AveBlank<- rowMeans(Peaks.filled[, Blks])} else {Peaks.filled$AveBlank<- Blks[1]}
Peaks.filled$AveSmp<- rowMeans(Peaks.filled[, Samps])
Peaks.filled$AvePoo <- rowMeans(Peaks.filled[, Poos])
Peaks.filtered <- Peaks.filled %>%
filter(rt > rtminparam) %>%
filter(rt < rtmaxparam) %>%
filter(AvePoo > 10*AveBlank) %>%
filter(npeaks > (0.25*length(Samps))) %>%
filter(AvePoo > 5000)
setwd(as.character(Dirs[Fraction, "ResultsDIR"]))
write.csv(Peaks.filtered, "xset.filtered.csv")
}
```
#Now do BMIS on each fraction separately