Mercurial > repos > iuc > dada2_mergepairs

comparison test-data/gentest.R @ 3:0b884b080bb6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/dada2 commit 8533fe71d1d50f09348da2dc34941724407a1ffe"
author iuc
date Tue, 14 Jul 2020 07:41:00 -0400
parents 69900ffd3b8e
children 84743da21318
comparison
equal deleted inserted replaced
2:69900ffd3b8e 3:0b884b080bb6
1 library(dada2, quietly=T) 1 library(dada2, quietly = T)
2 library(ggplot2, quietly=T) 2 library(ggplot2, quietly = T)
3 3
4 sample.names <- c('F3D0_S188_L001', 'F3D141_S207_L001') 4 sample_names <- c("F3D0_S188_L001", "F3D141_S207_L001")
5 fwd <- c('F3D0_S188_L001_R1_001.fastq.gz', 'F3D141_S207_L001_R1_001.fastq.gz') 5 fwd <- c("F3D0_S188_L001_R1_001.fastq.gz", "F3D141_S207_L001_R1_001.fastq.gz")
6 rev <- c('F3D0_S188_L001_R2_001.fastq.gz', 'F3D141_S207_L001_R2_001.fastq.gz') 6 rev <- c("F3D0_S188_L001_R2_001.fastq.gz", "F3D141_S207_L001_R2_001.fastq.gz")
7 7
8 filt.fwd <- c('filterAndTrim_F3D0_R1.fq.gz', 'filterAndTrim_F3D141_R1.fq.gz') 8 filt_fwd <- c("filterAndTrim_F3D0_R1.fq.gz", "filterAndTrim_F3D141_R1.fq.gz")
9 filt.rev <- c('filterAndTrim_F3D0_R2.fq.gz', 'filterAndTrim_F3D141_R2.fq.gz') 9 filt_rev <- c("filterAndTrim_F3D0_R2.fq.gz", "filterAndTrim_F3D141_R2.fq.gz")
10 10
11 print("filterAndTrim") 11 print("filterAndTrim")
12 12
13 for(i in 1:length(fwd)){ 13 for (i in seq_len(fwd)) {
14 ftout <- filterAndTrim(fwd[i], filt.fwd[i], rev[i], filt.rev[i]) 14 ftout <- dada2::filterAndTrim(fwd[i], filt_fwd[i], rev[i], filt_rev[i])
15 b <- paste(strsplit(fwd[i], ".", fixed=T)[[1]][1], "tab", sep=".") 15 b <- paste(strsplit(fwd[i], ".", fixed = T)[[1]][1], "tab", sep = ".")
16 write.table(ftout, b, quote=F, sep="\t", col.names=NA) 16 write.table(ftout, b, quote = F, sep = "\t", col.names = NA)
17 } 17 }
18 18
19 # In the test only the 1st data set is used 19 # In the test only the 1st data set is used
20 t <- data.frame() 20 t <- data.frame()
21 t <- rbind(t, ftout[1,]) 21 t <- rbind(t, ftout[1, ])
22 colnames(t) <- colnames(ftout) 22 colnames(t) <- colnames(ftout)
23 rownames(t) <- rownames(ftout)[1] 23 rownames(t) <- rownames(ftout)[1]
24 write.table(t, "filterAndTrim.tab", quote=F, sep="\t", col.names=NA) 24 write.table(t, "filterAndTrim.tab", quote = F, sep = "\t", col.names = NA)
25 25
26 names(fwd) <- sample.names 26 names(fwd) <- sample_names
27 names(rev) <- sample.names 27 names(rev) <- sample_names
28 names(filt.fwd) <- sample.names 28 names(filt_fwd) <- sample_names
29 names(filt.rev) <- sample.names 29 names(filt_rev) <- sample_names
30 30
31 # Plot quality profile (just for one file, Galaxy compares with sim_size) 31 # Plot quality profile (just for one file, Galaxy compares with sim_size)
32 print("plots") 32 print("plots")
33 qp <- plotQualityProfile(fwd) 33 qp <- dada2::plotQualityProfile(fwd)
34 ggsave('qualityProfile_fwd.pdf', qp, width = 20,height = 15,units = c("cm")) 34 ggsave("qualityProfile_fwd.pdf", qp, width = 20, height = 15, units = c("cm"))
35 qp <- plotQualityProfile(rev) 35 qp <- dada2::plotQualityProfile(rev)
36 ggsave('qualityProfile_rev.pdf', qp, width = 20,height = 15,units = c("cm")) 36 ggsave("qualityProfile_rev.pdf", qp, width = 20, height = 15, units = c("cm"))
37 qp <- plotQualityProfile(fwd[1]) 37 qp <- dada2::plotQualityProfile(fwd[1])
38 ggsave('qualityProfile.pdf', qp, width = 20,height = 15,units = c("cm")) 38 ggsave("qualityProfile.pdf", qp, width = 20, height = 15, units = c("cm"))
39 39
40 # Plot complexity (just for one file, Galaxy compares with sim_size) 40 # Plot complexity (just for one file, Galaxy compares with sim_size)
41 41
42 cp <- plotComplexity(fwd) 42 cp <- dada2::plotComplexity(fwd)
43 ggsave('complexity_fwd.pdf', cp, width = 20,height = 15,units = c("cm")) 43 ggsave("complexity_fwd.pdf", cp, width = 20, height = 15, units = c("cm"))
44 cp <- plotComplexity(rev) 44 cp <- dada2::plotComplexity(rev)
45 ggsave('complexity_rev.pdf', cp, width = 20,height = 15,units = c("cm")) 45 ggsave("complexity_rev.pdf", cp, width = 20, height = 15, units = c("cm"))
46 cp <- plotComplexity(fwd[1]) 46 cp <- dada2::plotComplexity(fwd[1])
47 ggsave('complexity.pdf', cp, width = 20,height = 15,units = c("cm")) 47 ggsave("complexity.pdf", cp, width = 20, height = 15, units = c("cm"))
48 48
49 49
50 # learn Errors 50 # learn Errors
51 print("learnErrors") 51 print("learnErrors")
52 err.fwd <- learnErrors(filt.fwd) 52 err_fwd <- dada2::learnErrors(filt_fwd)
53 saveRDS(err.fwd, file='learnErrors_R1.Rdata') 53 saveRDS(err_fwd, file = "learnErrors_R1.Rdata")
54 plot <- plotErrors(err.fwd) 54 plot <- dada2::plotErrors(err_fwd)
55 ggsave('learnErrors_R1.pdf', plot, width = 20,height = 15,units = c("cm")) 55 ggsave("learnErrors_R1.pdf", plot, width = 20, height = 15, units = c("cm"))
56 56
57 err.rev <- learnErrors(filt.rev) 57 err_rev <- dada2::learnErrors(filt_rev)
58 saveRDS(err.rev, file='learnErrors_R2.Rdata') 58 saveRDS(err_rev, file = "learnErrors_R2.Rdata")
59 plot <- plotErrors(err.rev) 59 plot <- dada2::plotErrors(err_rev)
60 ggsave('learnErrors.pdf', plot, width = 20,height = 15,units = c("cm")) 60 ggsave("learnErrors.pdf", plot, width = 20, height = 15, units = c("cm"))
61 61
62 # dada 62 # dada
63 print("dada") 63 print("dada")
64 dada.fwd <- dada(filt.fwd, err.fwd) 64 dada_fwd <- dada2::dada(filt_fwd, err_fwd)
65 dada.rev <- dada(filt.rev, err.rev) 65 dada_rev <- dada2::dada(filt_rev, err_rev)
66 for( id in sample.names ){ 66 for (id in sample_names) {
67 saveRDS(dada.fwd[[id]], file=paste("dada_", id,"_R1.Rdata", sep="")) 67 saveRDS(dada_fwd[[id]], file = paste("dada_", id, "_R1.Rdata", sep = ""))
68 saveRDS(dada.rev[[id]], file=paste("dada_", id,"_R2.Rdata", sep="")) 68 saveRDS(dada_rev[[id]], file = paste("dada_", id, "_R2.Rdata", sep = ""))
69 } 69 }
70 70
71 # merge pairs 71 # merge pairs
72 print("mergePairs") 72 print("mergePairs")
73 merged <- mergePairs(dada.fwd, filt.fwd, dada.rev, filt.rev) 73 merged <- dada2::mergePairs(dada_fwd, filt_fwd, dada_rev, filt_rev)
74 for( id in sample.names ){ 74 for (id in sample_names) {
75 saveRDS(merged[[id]], file=paste("mergePairs_", id,".Rdata", sep="")) 75 saveRDS(merged[[id]], file = paste("mergePairs_", id, ".Rdata", sep = ""))
76 } 76 }
77 77
78 78
79 # make sequence table 79 # make sequence table
80 print("makeSequenceTable") 80 print("makeSequenceTable")
81 seqtab <- makeSequenceTable(merged) 81 seqtab <- makeSequenceTable(merged)
82 write.table(t(seqtab), file="makeSequenceTable.tab", quote=F, sep="\t", row.names = T, col.names = NA) 82 write.table(t(seqtab), file = "makeSequenceTable.tab", quote = F, sep = "\t", row.names = T, col.names = NA)
83 83
84 reads.per.seqlen <- tapply(colSums(seqtab), factor(nchar(getSequences(seqtab))), sum) 84 reads_per_seqlen <- tapply(colSums(seqtab), factor(nchar(getSequences(seqtab))), sum)
85 df <- data.frame(length=as.numeric(names(reads.per.seqlen)), count=reads.per.seqlen) 85 df <- data.frame(length = as.numeric(names(reads_per_seqlen)), count = reads_per_seqlen)
86 pdf( 'makeSequenceTable.pdf' ) 86 pdf("makeSequenceTable.pdf")
87 ggplot(data=df, aes(x=length, y=count)) + 87 ggplot(data = df, aes(x = length, y = count)) +
88 geom_col() + 88 geom_col() +
89 theme_bw() 89 theme_bw()
90 bequiet <- dev.off() 90 bequiet <- dev.off()
91 91
92 # remove bimera 92 # remove bimera
93 print("removeBimera") 93 print("removeBimera")
94 seqtab.nochim <- removeBimeraDenovo(seqtab) 94 seqtab_nochim <- dada2::removeBimeraDenovo(seqtab)
95 write.table(t(seqtab), file="removeBimeraDenovo.tab", quote=F, sep="\t", row.names = T, col.names = NA) 95 write.table(t(seqtab), file = "removeBimeraDenovo.tab", quote = F, sep = "\t", row.names = T, col.names = NA)
96 96
97 # assign taxonomy/species 97 # assign taxonomy/species
98 tl <- 'Level1,Level2,Level3,Level4,Level5' 98 tl <- "Level1,Level2,Level3,Level4,Level5"
99 tl <- strsplit(tl, ",")[[1]] 99 tl <- strsplit(tl, ",")[[1]]
100 100
101 set.seed(42) 101 set.seed(42)
102 print("assignTaxonomyAndSpecies") 102 print("assignTaxonomyAndSpecies")
103 taxa <- assignTaxonomy(seqtab.nochim, 'reference.fa.gz', outputBootstraps = T, taxLevels=tl, multithread = 1) 103 taxa <- dada2::assignTaxonomy(seqtab_nochim, "reference.fa.gz", outputBootstraps = T, taxLevels = tl, multithread = 1)
104 104
105 taxa$tax <- addSpecies(taxa$tax, 'reference_species.fa.gz') 105 taxa$tax <- dada2::addSpecies(taxa$tax, "reference_species.fa.gz")
106 write.table(taxa$tax, file = 'assignTaxonomyAddspecies.tab', quote = F, sep = "\t", row.names = T, col.names = NA) 106 write.table(taxa$tax, file = "assignTaxonomyAddspecies.tab", quote = F, sep = "\t", row.names = T, col.names = NA)
107 107
108 write.table(taxa$boot, file = 'assignTaxonomyAddspecies_boot.tab', quote = F, sep = "\t", row.names = T, col.names = NA) 108 write.table(taxa$boot, file = "assignTaxonomyAddspecies_boot.tab", quote = F, sep = "\t", row.names = T, col.names = NA)
109 109
110 110
111 ## Generate extra test data for parameter testing 111 ## Generate extra test data for parameter testing
112 print("alternatives") 112 print("alternatives")
113 filterAndTrim(fwd, c('filterAndTrim_single_F3D0_R1.fq.gz', 'filterAndTrim_single_F3D141_R1.fq.gz'), rm.phix = T, orient.fwd = 'TACGG') 113 dada2::filterAndTrim(fwd, c("filterAndTrim_single_F3D0_R1.fq.gz", "filterAndTrim_single_F3D141_R1.fq.gz"), rm.phix = T, orient.fwd = "TACGG")
114 114
115 filterAndTrim(fwd, c('filterAndTrim_single_trimmers_F3D0_R1.fq.gz', 'filterAndTrim_single_trimmers_F3D141_R1.fq.gz'), truncQ = 30, truncLen = 2, trimLeft = 150, trimRight = 2) 115 dada2::filterAndTrim(fwd, c("filterAndTrim_single_trimmers_F3D0_R1.fq.gz", "filterAndTrim_single_trimmers_F3D141_R1.fq.gz"), truncQ = 30, truncLen = 2, trimLeft = 150, trimRight = 2)
116 116
117 filterAndTrim(fwd, c('filterAndTrim_single_filters_F3D0_R1.fq.gz', 'filterAndTrim_single_filters_F3D141_R1.fq.gz'), maxLen = 255, minLen = 60, maxN = 100, minQ = 13, maxEE = 1) 117 dada2::filterAndTrim(fwd, c("filterAndTrim_single_filters_F3D0_R1.fq.gz", "filterAndTrim_single_filters_F3D141_R1.fq.gz"), maxLen = 255, minLen = 60, maxN = 100, minQ = 13, maxEE = 1)
118 118
119 119
120 merged_nondef <- mergePairs(dada.fwd, filt.fwd, dada.rev, filt.rev, minOverlap = 8, maxMismatch = 1, justConcatenate = TRUE, trimOverhang = TRUE) 120 merged_nondef <- dada2::mergePairs(dada_fwd, filt_fwd, dada_rev, filt_rev, minOverlap = 8, maxMismatch = 1, justConcatenate = TRUE, trimOverhang = TRUE)
121 for( id in sample.names ){ 121 for (id in sample_names) {
122 saveRDS(merged_nondef[[id]], file=paste("mergePairs_", id,"_nondefault.Rdata", sep="")) 122 saveRDS(merged_nondef[[id]], file = paste("mergePairs_", id, "_nondefault.Rdata", sep = ""))
123 } 123 }
124 rb.dada.fwd <- removeBimeraDenovo(dada.fwd[["F3D0_S188_L001"]]) 124 rb_dada_fwd <- dada2::removeBimeraDenovo(dada_fwd[["F3D0_S188_L001"]])
125 write.table(rb.dada.fwd, file = 'removeBimeraDenovo_F3D0_dada_uniques.tab', quote = F, sep = "\t", row.names = T, col.names = F) 125 write.table(rb_dada_fwd, file = "removeBimeraDenovo_F3D0_dada_uniques.tab", quote = F, sep = "\t", row.names = T, col.names = F)
126 126
127 rb.merged <- removeBimeraDenovo(merged, method="pooled") 127 rb_merged <- dada2::removeBimeraDenovo(merged, method = "pooled")
128 saveRDS(rb.merged, file='removeBimeraDenovo_F3D0_mergepairs.Rdata') 128 saveRDS(rb_merged, file = "removeBimeraDenovo_F3D0_mergepairs.Rdata")
129 129
130 # SeqCounts 130 # SeqCounts
131 getN <- function(x){ sum(getUniques(x)) } 131 get_n <- function(x) {
132 132 sum(dada2::getUniques(x))
133 read.uniques <- function ( fname ) {
134 p <- read.table(fname, header=F, sep="\t")
135 n <-x[,2]
136 names(n)<-x[,1]
137 } 133 }
138 134
139
140 print("seqCounts ft") 135 print("seqCounts ft")
141 samples = list() 136 samples <- list()
142 samples[["F3D0_S188_L001_R1_001.tab"]] <- read.table("F3D0_S188_L001_R1_001.tab", header=T, sep="\t", row.names=1) 137 samples[["F3D0_S188_L001_R1_001.tab"]] <- read.table("F3D0_S188_L001_R1_001.tab", header = T, sep = "\t", row.names = 1)
143 dname <- "filter" 138 dname <- "filter"
144 tdf <- samples[["F3D0_S188_L001_R1_001.tab"]] 139 tdf <- samples[["F3D0_S188_L001_R1_001.tab"]]
145 names(tdf) <- paste( dname, names(tdf) ) 140 names(tdf) <- paste(dname, names(tdf))
146 tdf <- cbind( data.frame(samples=names( samples )), tdf) 141 tdf <- cbind(data.frame(samples = names(samples)), tdf)
147 write.table(tdf, "seqCounts_filter.tab", quote=F, sep="\t", row.names = F, col.names = T) 142 write.table(tdf, "seqCounts_filter.tab", quote = F, sep = "\t", row.names = F, col.names = T)
148 143
149 samples = list() 144 samples <- list()
150 samples[["F3D0_S188_L001_R1_001.tab"]] <- read.table("F3D0_S188_L001_R1_001.tab", header=T, sep="\t", row.names=1) 145 samples[["F3D0_S188_L001_R1_001.tab"]] <- read.table("F3D0_S188_L001_R1_001.tab", header = T, sep = "\t", row.names = 1)
151 samples[["F3D141_S207_L001_R1_001.tab"]] <- read.table("F3D141_S207_L001_R1_001.tab", header=T, sep="\t", row.names=1) 146 samples[["F3D141_S207_L001_R1_001.tab"]] <- read.table("F3D141_S207_L001_R1_001.tab", header = T, sep = "\t", row.names = 1)
152 dname <- "filter" 147 dname <- "filter"
153 tdf <- samples[["F3D0_S188_L001_R1_001.tab"]] 148 tdf <- samples[["F3D0_S188_L001_R1_001.tab"]]
154 tdf <- rbind(tdf, samples[["F3D141_S207_L001_R1_001.tab"]]) 149 tdf <- rbind(tdf, samples[["F3D141_S207_L001_R1_001.tab"]])
155 names(tdf) <- paste( dname, names(tdf) ) 150 names(tdf) <- paste(dname, names(tdf))
156 tdf <- cbind( data.frame(samples=names( samples )), tdf) 151 tdf <- cbind(data.frame(samples = names(samples)), tdf)
157 write.table(tdf, "seqCounts_filter_both.tab", quote=F, sep="\t", row.names = F, col.names = T) 152 write.table(tdf, "seqCounts_filter_both.tab", quote = F, sep = "\t", row.names = F, col.names = T)
158 153
159 print("seqCounts dada") 154 print("seqCounts dada")
160 samples = list() 155 samples <- list()
161 samples[["dada_F3D0_S188_L001_R1.Rdata"]] <- readRDS('dada_F3D0_S188_L001_R1.Rdata') 156 samples[["dada_F3D0_S188_L001_R1.Rdata"]] <- readRDS("dada_F3D0_S188_L001_R1.Rdata")
162 samples[["dada_F3D141_S207_L001_R1.Rdata"]] <- readRDS('dada_F3D141_S207_L001_R1.Rdata') 157 samples[["dada_F3D141_S207_L001_R1.Rdata"]] <- readRDS("dada_F3D141_S207_L001_R1.Rdata")
163 dname <- "dadaF" 158 dname <- "dadaF"
164 tdf <- data.frame( samples = names(samples) ) 159 tdf <- data.frame(samples = names(samples))
165 tdf[[ dname ]] <- sapply(samples, getN) 160 tdf[[dname]] <- sapply(samples, get_n)
166 write.table(tdf, "seqCounts_dadaF.tab", quote=F, sep="\t", row.names = F, col.names = T) 161 write.table(tdf, "seqCounts_dadaF.tab", quote = F, sep = "\t", row.names = F, col.names = T)
167 162
168 print("seqCounts mp") 163 print("seqCounts mp")
169 samples = list() 164 samples <- list()
170 samples[["mergePairs_F3D0_S188_L001.Rdata"]] <- readRDS('mergePairs_F3D0_S188_L001.Rdata') 165 samples[["mergePairs_F3D0_S188_L001.Rdata"]] <- readRDS("mergePairs_F3D0_S188_L001.Rdata")
171 samples[["mergePairs_F3D141_S207_L001.Rdata"]] <- readRDS('mergePairs_F3D141_S207_L001.Rdata') 166 samples[["mergePairs_F3D141_S207_L001.Rdata"]] <- readRDS("mergePairs_F3D141_S207_L001.Rdata")
172 dname <- "merge" 167 dname <- "merge"
173 tdf <- data.frame( samples = names(samples) ) 168 tdf <- data.frame(samples = names(samples))
174 tdf[[ dname ]] <- sapply(samples, getN) 169 tdf[[dname]] <- sapply(samples, get_n)
175 write.table(tdf, "seqCounts_merge.tab", quote=F, sep="\t", row.names = F, col.names = T) 170 write.table(tdf, "seqCounts_merge.tab", quote = F, sep = "\t", row.names = F, col.names = T)
176 171
177 print("seqCounts st") 172 print("seqCounts st")
178 samples = list() 173 samples <- list()
179 samples <- t(as.matrix( read.table("makeSequenceTable.tab", header=T, sep="\t", row.names=1) )) 174 samples <- t(as.matrix(read.table("makeSequenceTable.tab", header = T, sep = "\t", row.names = 1)))
180 dname <- "seqtab" 175 dname <- "seqtab"
181 tdf <- data.frame( samples = row.names(samples) ) 176 tdf <- data.frame(samples = row.names(samples))
182 tdf[[ dname ]] <- rowSums(samples) 177 tdf[[dname]] <- rowSums(samples)
183 write.table(tdf, "seqCounts_seqtab.tab", quote=F, sep="\t", row.names = F, col.names = T) 178 write.table(tdf, "seqCounts_seqtab.tab", quote = F, sep = "\t", row.names = F, col.names = T)
184 179
185 print("seqCounts rb") 180 print("seqCounts rb")
186 samples = list() 181 samples <- list()
187 samples <- t(as.matrix( read.table("removeBimeraDenovo.tab", header=T, sep="\t", row.names=1) )) 182 samples <- t(as.matrix(read.table("removeBimeraDenovo.tab", header = T, sep = "\t", row.names = 1)))
188 dname <- "nochim" 183 dname <- "nochim"
189 tdf <- data.frame( samples = row.names(samples) ) 184 tdf <- data.frame(samples = row.names(samples))
190 tdf[[ dname ]] <- rowSums(samples) 185 tdf[[dname]] <- rowSums(samples)
191 write.table(tdf, "seqCounts_nochim.tab", quote=F, sep="\t", row.names = F, col.names = T) 186 write.table(tdf, "seqCounts_nochim.tab", quote = F, sep = "\t", row.names = F, col.names = T)
192