cutadapt bugfixes + vserch tutorial

Author: VascoElbrecht

.gitignore

@@ -1,16 +1,3 @@
-.Rhistory
-.Rapp.history
-
-*.history
-JAMP/R/.Rapp.history
-JAMP/R/.Rapp.history
-*.history
-*.history
-JAMP/R/.Rapp.history
-JAMP/R/log.txt
-JAMP/R/.Rapp.history
-*.history
-.DS_Store
-*.history
-*.history
-*.history
+*.Rhistory
+*.Rapp.history
+*.DS_Store

JAMP/R/Cutadapt.R

@@ -1,18 +1,13 @@
 # Cutadapt v0.1
 
-Cutadapt <- function(files="latest", forward=NA, reverse=NA, bothsides=F, anchoring=T, fastq=T, LDist=F, cores=1, delete_data=T, ambsequ=F, error=0.1, exe="cutadapt", exeU="usearch"){
+Cutadapt <- function(files="latest", forward=NA, reverse=NA, bothsides=F, anchoring=T, LDist=F, cores=1, delete_data=T, ambsequ=F, error=0.1, exe="cutadapt", exe2="vsearch", fastq=T){
 
 # do checks
 FW_only <- F
 if (is.na(reverse[1])){FW_only <- T}
 
 
 
-
-
-
-
-
 folder <- Core(module="Cutadapt", delete_data=delete_data)
 cat(file="log.txt", c("Module Version: v0.2", "\n"), append=T, sep="\n")
 
@@ -39,6 +34,7 @@ if(!fastq){ #rename to fasta if fastq=F
 new_names <- sub(".fastq", ".fasta", new_names)
 }
 
+
 # get primer sequences / names
 primers <- read.csv(paste(system.file(package="JAMP"), "/primers.csv", sep=""), stringsAsFactors=F)
 
@@ -86,15 +82,15 @@ rw[i] <- paste(rev(comp(unlist(strsplit(rw[i], "")), forceToLower=F, ambiguous=T
 
 # add: write down used primers in log!
 if (FW_only){
-cmd1 <- paste("-g ", if(anchoring){"^"}, fw, " -o \"", new_names, "\" \"", files, "\"", " -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") # forward adapter
+cmd1 <- paste("-g ", if(anchoring){"^"}, fw, " -o \"", new_names, "\" \"", files, "\"", " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") # forward adapter
 files_to_delete <- c(files_to_delete, new_names)
 
 temp <- paste("Starting to remove adapters (primers) on forward direction only ", length(cmd1), " files:", sep="")
 cat(file="log.txt", temp , append=T, sep="\n")
 message(temp)
 }else{
-cmd1 <- paste("-g ", if(anchoring){"^"}, fw, " -o ", folder, "/_data/temp.txt \"", files, "\"", " -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") # forward adapter
-cmd2 <- paste("-a ", rw, if(anchoring){"$"}," -o \"", new_names, "\" ", folder, "/_data/temp.txt -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") #rverse adapter
+cmd1 <- paste("-g ", if(anchoring){"^"}, fw, " -o ", folder, "/_data/temp.txt \"", files, "\"", " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") # forward adapter
+cmd2 <- paste("-a ", rw, if(anchoring){"$"}," -o \"", new_names, "\" ", folder, "/_data/temp.txt ", " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error, sep="") #rverse adapter
 files_to_delete <- c(files_to_delete, new_names)
 
 temp <- paste("Starting to remove adapters (primers) on both ends in ", length(cmd1), " files:", sep="")
@@ -129,28 +125,28 @@ cat(file=log_names[i], A, append=T, sep="\n")
 #cat(file=log_names[i], A, append=T, sep="\n")
 
 # trimm RW
-rev_primer <- paste("-a ", rw[i], if(anchoring){"$"}," -o ", folder, "/_data/temp_A.txt ", folder, "/_data/temp.txt -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, sep="", if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error) #rverse adapter
+rev_primer <- paste("-a ", rw[i], if(anchoring){"$"}," -o ", folder, "/_data/temp_A.txt ", folder, "/_data/temp.txt ", " --discard-untrimmed --cores=", cores, sep="", if(ambsequ){" --match-read-wildcards"}, " --error-rate ", error) #rverse adapter
 
 A <- system2(exe, rev_primer, stdout=T, stderr=T)
 cat(file=log_names[i], A, append=T, sep="\n")
 
 # build revcomp
 
 revcom_cmd <- paste("-fastx_revcomp ", files[i],  if(fastq){" -fastqout "} else {" -fastaout "}, folder, "/_data/temp2_RC.txt  -label_suffix _RC",  sep="")
-A <- system2(exeU, revcom_cmd, stdout=T, stderr=T)
+A <- system2(exe2, revcom_cmd, stdout=T, stderr=T)
 
-cat(file=log_names[i], paste("usearch ", revcom_cmd, sep=""), append=T, sep="\n")
+cat(file=log_names[i], paste(exe2, revcom_cmd, sep=""), append=T, sep="\n")
 cat(file=log_names[i], A, append=T, sep="\n")
 
 # cut primers again! simmillar commands
 
-fw_primer_cmd <- paste("-g ", if(anchoring){"^"}, fw[i], " -o ", folder, "/_data/temp.txt ", folder, "/_data/temp2_RC.txt", " -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, sep="", " --error-rate ", error) # forward adapter
+fw_primer_cmd <- paste("-g ", if(anchoring){"^"}, fw[i], " -o ", folder, "/_data/temp.txt ", folder, "/_data/temp2_RC.txt", " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, sep="", " --error-rate ", error) # forward adapter
 
 A <- system2(exe, fw_primer_cmd, stdout=T, stderr=T)
 cat(file=log_names[i], A, append=T, sep="\n")
 
 
-rev_primer <- paste("-a ", rw[i], if(anchoring){"$"}," -o ", folder, "/_data/temp_B.txt ", folder, "/_data/temp.txt -f ", if(fastq){"fastq"}else{"fasta"}, " --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, sep="", " --error-rate ", error) #rverse adapter
+rev_primer <- paste("-a ", rw[i], if(anchoring){"$"}," -o ", folder, "/_data/temp_B.txt ", folder, "/_data/temp.txt --discard-untrimmed --cores=", cores, if(ambsequ){" --match-read-wildcards"}, sep="", " --error-rate ", error) #rverse adapter
 
 A <- system2(exe, rev_primer, stdout=T, stderr=T)
 cat(file=log_names[i], A, append=T, sep="\n")

JAMP/R/Demultiplexing_shifted.R

@@ -90,7 +90,7 @@ write.csv(combos2, file= combinations, row.names=F)
 
 combos <- combos2
 
-temp <- paste("FileName detected! \"R1.txt\" and \"R2.txt\" added automatically to filenames. Please check the table used for primer combinations / filenames in \"", folder, "/combinations_auto_edited.csv\" for details.", sep="", "\n")
+temp <- paste("FileName detected! \"R1.fastq\" and \"R2.fastq\" added automatically to filenames. Please check the table used for primer combinations / filenames in \"", folder, "/combinations_auto_edited.csv\" for details.", sep="", "\n")
 message(temp)
 cat(file="log.txt", temp, append=T, sep="\n")
 

JAMP/R/Minmax.R

@@ -39,7 +39,7 @@ max <- plusminus[1] + plusminus[2]
 
 
 # make cmd
-cmd <- paste("\"", files, "\" -o \"", folder, "/", new_names, "\" -f ", if(fastq){"fastq"}else{"fasta"}, if(!is.na(min[1])){" -m "}, if(!is.na(min[1])){min}, if(!is.na(min[1])){" -M "}, if(!is.na(min[1])){max}, sep="")
+cmd <- paste("\"", files, "\" -o \"", folder, "/", new_names, "\"", if(!is.na(min[1])){" -m "}, if(!is.na(min[1])){min}, if(!is.na(min[1])){" -M "}, if(!is.na(min[1])){max}, sep="")
 
 files_to_delete <- c(files_to_delete, paste(folder, "/", new_names, sep=""))
 

Tutorial/package_tutorial.R

@@ -1,10 +1,6 @@
-# 170201 jamp pipeline, short tutorial
+# 12052021 jamp pipeline, short tutorial
 #
-# download the raw data from and place it into the tutorial folder (Study we are looking at https://peerj.com/articles/3006/) 
-# MiSeq Run1, R1 direction:
-# https://dx.doi.org/10.6084/m9.figshare.4039821.v1
-# MiSeq Run1, R2 direction:
-# https://dx.doi.org/10.6084/m9.figshare.4039860.v1
+
 # 
 # Installing dependencies needed fro JAMP
 install.packages(c("bold", "XML", "seqinr", "devtools", "fastqcr"), dependencies=T)
@@ -13,8 +9,21 @@ library("devtools")
 install_github("VascoElbrecht/JAMP", subdir="JAMP")
 
 
+# set workin directory
+setwd("~/Desktop/JAMP_pipeline/Tutorial")
+# download the raw data from and place it into the tutorial folder (Study we are looking at https://peerj.com/articles/3006/) 
+# MiSeq Run1, R1 direction:
+# https://dx.doi.org/10.6084/m9.figshare.4039821.v1
+# MiSeq Run1, R2 direction:
+# https://dx.doi.org/10.6084/m9.figshare.4039860.v1
+
+download.file("https://ndownloader.figshare.com/files/6503952", "16_S10_L001_R1_001_run1.fastq.gz")
+download.file("https://ndownloader.figshare.com/files/6503991", "16_S10_L001_R2_001_run1.fastq.gz")
+system2("gunzip", "16_S10_L001_R1_001_run1.fastq.gz")
+system2("gunzip", "16_S10_L001_R2_001_run1.fastq.gz")
+
+
 
-setwd("~/Desktop/JAMP_pipeline")
 list.files()
 
 library("JAMP")
@@ -30,33 +39,41 @@ Remove_last_folder()
 # In this example we are dealing with sequence raw data that is not jet demultiplexed. To demultiplex run:
 Demultiplexing_shifted(file1="16_S10_L001_R1_001_run1.fastq", file2="16_S10_L001_R2_001_run1.fastq", tags="_converter/indexe_1.csv", combinations="_converter/combos_1.csv")
 
-# Paired end merging
-U_merge_PE(fastq_pctid=75)
+
 
 # check for PhiX
-system2("usearch", "-usearch_global A_Demultiplexing_shifted/_data/N_debris_r1.txt -db PhiX.fasta -id 0.9 -strand both -blast6out PhiX_table.txt")
+# only subsample 10000 reads
+if(T){
+system2("vsearch", "-fastx_subsample A_Demultiplexing_shifted/_data/N_debris_R1.fastq -fastaout A_Demultiplexing_shifted/_data/N_debris_R1.fasta -sample_size 10000")
+} else { # check all reads in read 1 for PhiX
+system2("paste", " - - - - < A_Demultiplexing_shifted/_data/N_debris_R1.fastq | cut -f 1,2 | sed 's/^@/>/' | tr \"\t\" \"\n\" > A_Demultiplexing_shifted/_data/N_debris_R1.fasta")
+}
+system2("vsearch", "-usearch_global A_Demultiplexing_shifted/_data/N_debris_R1.fasta -db PhiX.fasta -id 0.9 -strand both -blast6out PhiX_table.txt -maxrejects 1 -maxaccepts 1")
 
 
+# Paired end merging
+Merge_PE()
+
 
-# trimm primers
-Cutadapt(forward="GGWACWGGWTGAACWGTWTAYCCYCC", # mlCOIintF
-reverse="TANACYTCNGGRTGNCCRAARAAYCA", bothsides=T) # jgHCO, I (inosin) replaced with N
+# trimm primers (mlCOIintF and jgHCO)
+Cutadapt(forward="GGWACWGGWTGAACWGTWTAYCCYCC", reverse="TAIACYTCIGGRTGICCRAARAAYCA", bothsides=T)
 #by using "bothsides=T", forward or reverse primers are detected on both ends. This is not nessesary for fusion primers.
 
 # discard with non target length
 Minmax(min=(313-10), max=(313+10))
 
 # discard reads above 1 expected error
-U_max_ee(max_ee=1)
+Max_ee(max_ee=1)
 
 # subsample to lowest sample size, should be done if samples are widely different in sequencing depth (as one starts with)
+# needs to be umdated, right now still on userch! Vserch not supported right now
 U_subset(sample_size=60000)
 
 
 
 
 #cluster OTUs
-U_cluster_otus(filter=0.01)
+Cluster_otus(filter=0.01)
 file.rename("G_U_cluster_otus", "G_U_cluster_otus - 60k")
 
 #cluster OTUs (without subsetting)