Workflow Steps and Code Snippets

267 tagged steps and code snippets that match keyword Cutadapt

Snakemake based analysis pipeline to identify m6As from eCLIP data 

107
108
109
110
111
112
113
114
shell:
    """
    rna_seq="{params.rna_seq}"; \
    rna_seq_rc=`echo -e ">seq\n${{rna_seq}}" | seqkit seq -p -r -s -t RNA --rna2dna --quiet`; \
    randomer_size="{params.randomer_size}"; \
    randomer_size=`echo "{{"${{randomer_size}}"}}"`; \
    cutadapt {params.args} -j {threads} -a "ssDNA=${{rna_seq_rc}};min_overlap=5...N${{randomer_size}}{params.truseq_read1};min_overlap=15" -A ssRNA={params.rna_seq} -A TruSeq={params.truseq_read2} -o {output.read1} -p {output.read2} {input} > {log}
    """
SnakeMake Cutadapt seqkit From line 107 of workflow/Snakefile 
253
254
255
256
257
258
259
260
shell:
    """
    rna_seq="{params.rna_seq}"; \
    rna_seq_rc=`echo -e ">seq\n${{rna_seq}}" | seqkit seq -p -r -s -t RNA --rna2dna --quiet`; \
    randomer_size="{params.randomer_size}"; \
    randomer_size=`echo "{{"${{randomer_size}}"}}"`; \
    cutadapt {params.args} -j {threads} -a "ssDNA=${{rna_seq_rc}};min_overlap=5...N${{randomer_size}}{params.truseq_read1};min_overlap=15" -A ssRNA={params.rna_seq} -A TruSeq={params.truseq_read2} -o {output.read1} -p {output.read2} {input} > {log}
    """
SnakeMake Cutadapt seqkit From line 253 of workflow/Snakefile

Bioinformatics pipeline for processing 16s marker gene metagenomic sequence data. 

30
31
32
33
34
shell:
        "mkdir -p output/temp/cutadapt/logs; cutadapt -e 0 -O 10 -m 50 -n 2 -g {config[fwd_primer]} "
        "-G {config[rev_primer]} -a {config[rev_primer_rc]} "
        "-A {config[fwd_primer_rc]} -o {output.r1} -p {output.r2} "
        "{input.r1} {input.r2} > output/temp/cutadapt/logs/{params.pre}.cutadapt.log"
SnakeMake Cutadapt From line 30 of main/Snakefile 
  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
library(dada2)
library(tidyverse)
library(colorRamps)

## Set variables
list_of_filenames <- snakemake@input$listfiles

## Set filtering parameters from config file
trimleft <- c(snakemake@config$trimleft_forward, snakemake@config$trimleft_reverse)
expected_errors <- c(snakemake@config$expected_errors_forward, snakemake@config$expected_errors_reverse)
truncate <- c(snakemake@config$truncate_forward, snakemake@config$truncate_reverse)
readlength <- snakemake@config$readlength

## Cutadapt setting
trimmed <- snakemake@config$run_cutadapt
path_to_raw <- snakemake@config$path

## Exploring parameter space options
trunc_param <- snakemake@config$explore_truncation_params
ee_param <- snakemake@config$explore_quality_params

## Make directory for quality plots
dir.create("output/quality_plots/")

#########################################################


## Make a vector of sample names
samples <- scan(list_of_filenames, what = "character")

## file names of forward and reverse reads, before quality filtering
if (trimmed == T){
	forward_reads <- file.path("output/temp/cutadapt", paste0(samples, "_r1_cutadapt.fastq.gz"))
	reverse_reads <- file.path("output/temp/cutadapt", paste0(samples, "_r2_cutadapt.fastq.gz"))
}

if (trimmed == F){
	forward_reads <- file.path(path_to_raw, paste0(samples, "_R1_001.fastq.gz"))
	reverse_reads <- file.path(path_to_raw, paste0(samples, "_R2_001.fastq.gz"))
}

## file names of forward and reverse reads, after quality filtering
filtered_forward_reads <- file.path("output/temp/filtered", paste0(samples, "_r1_filtered.fastq.gz"))
filtered_reverse_reads <- file.path("output/temp/filtered", paste0(samples, "_r2_filtered.fastq.gz"))

#########################################################
######## Step 0: Exploring filtering parameters

if (trunc_param == T){
## Explore truncation parameters
results <- NULL
## Select a set of samples at random to inspect
test <- sample(c(1:length(samples)), 12)

for (i in seq(from = readlength-45, to = readlength, by = 5)){
  for (j in seq(from = readlength-90, to = readlength, by = 10)){
  	truncparam <- c()
    out <- filterAndTrim(forward_reads[test],
                         filtered_forward_reads[test],
                         reverse_reads[test],
                         filtered_reverse_reads[test], 
						 truncLen=c(i,j),
                         maxEE=expected_errors, rm.phix=TRUE,
                         compress=TRUE, multithread=TRUE, trimLeft = trimleft)
    res <- data.frame(Sample = rownames(out), perc = out[,2]/out[,1], for_trunc = i, rev_trunc = j)
    results <- rbind(results, res)
  }
}

results <- results %>% separate(Sample, c("Name", "Sample"), sep = "_S")

gg <- ggplot(results, aes(x = for_trunc, y = perc, colour = as.factor(rev_trunc))) +
  geom_point(size = 2) +
  geom_line(size = 1) +
  scale_colour_manual(values = sample(primary.colors(20), 10), name = "Truncate Reverse") +
  xlab("Truncate Forward") +
  ylab("Percentage reads passed filtering") +
  ggtitle("Truncation parameters") +
  theme_minimal() +
  facet_wrap(~Name)


pdf(file.path("output", "truncation_parameters.pdf"))
print(gg)
dev.off()
}

if (ee_param == T){
## Explore expected error values
results <- NULL

for (i in 1:5){
	for (j in 1:5){
		out <- filterAndTrim(forward_reads[test], 
							filtered_forward_reads[test], 
							reverse_reads[test],
							filtered_reverse_reads[test], 
							truncLen=truncate,
							maxEE=c(i,j), rm.phix=TRUE,
							compress=TRUE, multithread=TRUE, trimLeft = trimleft)
		res <- data.frame(Sample = rownames(out), perc = out[,2]/out[,1], for_error = i, rev_error = j)
 		results <- rbind(results, res)

	}
}

results <- results %>% separate(Sample, c("Name", "Sample"), sep = "_S")

gg <- ggplot(results, aes(x = for_error, y = perc, colour = as.factor(rev_error))) +
	geom_point(size = 2) +
	geom_line(size = 1) +
	scale_colour_manual(values = rainbow(5, v = 0.8), name = "Error Rate Reverse") +
	xlab("Error Rate Forward") +
	ylab("Percentage reads passed filtering") +
	ggtitle("Expected error parameters") +
	theme_minimal() +
	facet_wrap(~Name)


pdf(file.path("output", "expected_error_parameters.pdf"))
print(gg)
dev.off()
}

#########################################################
####### Step 1: Quality filtering

out <- filterAndTrim(forward_reads, filtered_forward_reads, reverse_reads, filtered_reverse_reads, maxEE = expected_errors, multithread = TRUE, rm.phix=TRUE, trimLeft = trimleft, compress = TRUE, truncLen = truncate)

saveRDS(out, file.path("output/temp", "filt_out.rds"))

####### Step 2: Plot quality profiles
## Select 10 samples at random to inspect/print
toplot <- sample(c(1:length(samples)), 10)

## Plotting forward versus reverse quality
for (i in 1:10){
  pdf(paste("output/quality_plots/quality_", samples[toplot[i]], ".pdf", sep = ""))
  print(plotQualityProfile(c(filtered_forward_reads[toplot[i]], forward_reads[toplot[i]], filtered_reverse_reads[toplot[i]], reverse_reads[toplot[i]])))
  dev.off()
}
R ggplot2 tidyverse Cutadapt dada2 colorRamps From line 8 of scripts/filter.R

Automated pipeline for amplicon sequence analysis 

  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
import os
import subprocess
from sys import stdin
#import benchmark_utils
from benchmark_utils import countFasta

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'Y':'R', 'R':'Y','S':'S','W':'W','K':'M','M':'K','N':'N','B':'V','V':'B','D':'H','H':'D'} 
    bases = list(seq) 
    bases = [complement[base] for base in bases] 
    return ''.join(bases)


def reverse_complement(s):
    return complement(s[::-1])

#from Bio.Seq import Seq


primer_by_sample={}
uniq_primers={}
idx_fw_primer=-1   # default for qiime (col 3)
idx_rv_primer=-1  # new field 
idx_rv_revcomp_primer=-1
isRC = False
foundSample=False  
primer=""
if snakemake.config["primers"]["remove"].lower() == "metadata":
    with open(snakemake.input[1]) as mappingFile:
        l=0
        for line in mappingFile:
            l=l+1;
            columns = line.split('\t')
            #the header is always at row 1 and must contain these first 3 fields (qiime specs):
            #SampleID BarcodeSequence LinkerPrimerSequence Description
            if l==1 :
                c=0
                #Find target headers
                for col in columns:
                    if col == "ReversePrimer" or col == "LinkerPrimerSequenceReverse"  or col == "ReverseLinkerPrimerSequence"  or col == "RvLinkerPrimerSequence" or col == "ReversePrimerSequence" :
                        idx_rv_primer=c
                    elif col == "LinkerPrimerSequence":
                        idx_fw_primer=c
                    elif col == "ReverseLinkerPrimerSequenceRevCom"  or col == "ReversePrimerRevCom":
                        idx_rv_revcomp_primer=c
                        isRC=True
                    c=c+1
                if isRC:
                    idx_rv_primer=idx_rv_revcomp_primer 
            elif line.startswith(snakemake.params[4]):
                foundSample=True
                if idx_rv_primer != -1:
                    if isRC:
                        #fw_primer=columns[idx_fw_primer]
                        #rv_primer=columns[idx_rv_primer]
                        primer="-g "+columns[idx_fw_primer]+"..."+columns[idx_rv_primer]
                    else:
                        #fw_primer=columns[idx_fw_primer]
                        #rv_primer=reverse_complement(columns[idx_rv_primer])
                        primer="-g "+columns[idx_fw_primer]+"..."+reverse_complement(columns[idx_rv_primer])
                else:
                    #fw_primer=columns[idx_fw_primer]
                    primer="-g "+columns[idx_fw_primer]


    if not foundSample:
        print("\033[91m" +"No primers found for sample:"+ snakemake.params[4]+ " \033[0m")
        print("\033[91mPlease make sure to have the sample included in the mapping file: "+snakemake.input[1]+"  \033[0m")
        print("\033[91m Aborting the pipeline \033[0m")
        exit(1)

elif snakemake.config["primers"]["remove"].lower() == "cfg":
    primer="-g " + snakemake.config["primers"]["fw_primer"]
    if snakemake.config["primers"]["rv_primer"].len() > 2 :
        primer=primer+"..."+reverse_complement(snakemake.config["primers"]["rv_primer"]) 


discard = True
if "--discard-untrimmed" in snakemake.params[0]:
    extra=snakemake.params[0].replace("--discard-untrimmed","")
else: 
    extra=snakemake.params[0]
    discard = False

#This file will contain the untrimmed reads for the first pass
no_primer=" --untrimmed-output " + snakemake.params[2]+".tmp"

if snakemake.config["primers"]["remove"].lower() == "metadata":
    subprocess.run( ["cutadapt  "+ primer +" "+ extra+" -o "+snakemake.output[0] + ".1 "+ no_primer +" " + snakemake.input[0]+ ">"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
else:
    subprocess.run( ["cutadapt  "+ primer +" "+extra+" -o "+snakemake.output[0] + ".1 "+ no_primer +" " + snakemake.input[0]+ ">"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
#    primer=snakemake.config["cutadapt"]["adapters"]
#comment above line because we just add the primer generation in the elif above....


initialReads=countFasta(snakemake.input[0],False)
disscardedReads=countFasta(snakemake.params[2]+".tmp",False)

#The "extra" var returns to the original values in the sense that if the user wants to disscard reads
# this option will be present on the final cutadapt command 
extra=snakemake.params[0]
#if we disscarded reads
if disscardedReads>0:
    #reverse complement disscardedReads
    subprocess.run( ["vsearch --fastx_revcomp "+ snakemake.params[2]+".tmp  --fastaout "+ snakemake.params[2]+".tmp2"],stdout=subprocess.PIPE, shell=True)
    if snakemake.config["primers"]["remove"].lower() == "metadata":
        if discard:
        #Run cutadapt on disscarded reads
            subprocess.run( ["cutadapt  "+ primer +" "+ extra+" -o "+snakemake.output[0] + ".2 " + snakemake.params[2]+".tmp2"+ ">>"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
        else:
            print("Running second cutadapt")
            print("cutadapt  "+ primer +" "+ extra+" -o "+snakemake.output[0] + ".2 --untrimmed-output "+ snakemake.output[0] + ".3 " + snakemake.params[2]+".tmp2"+ ">>"+ snakemake.params[5])
            subprocess.run( ["cutadapt  "+ primer +" "+ extra+" -o "+snakemake.output[0] + ".2 --untrimmed-output "+ snakemake.output[0] + ".3 " + snakemake.params[2]+".tmp2"+ ">>"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
            #reverse complement untrimmed disscardedReads
            subprocess.run( ["vsearch --fastx_revcomp "+ snakemake.output[0]+".3  --fastaout "+ snakemake.params[2]+".tmp3"],stdout=subprocess.PIPE, shell=True)
    else:
        if discard:
            subprocess.run( ["cutadapt "+ primer  +" "+extra+" -o "+snakemake.output[0] + ".2 " + snakemake.params[2]+".tmp2 >>"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
        else:
            subprocess.run( ["cutadapt "+ primer  +" "+extra+" -o "+snakemake.output[0] + ".2 --untrimmed-output "+ snakemake.output[0] + ".3 "  + snakemake.params[2]+".tmp2 >>"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
            subprocess.run( ["vsearch --fastx_revcomp "+ snakemake.output[0]+".3  --fastaout "+ snakemake.params[2]+".tmp3"],stdout=subprocess.PIPE, shell=True)

    if discard:        
        #Concatenate results
        subprocess.run( ["cat "+snakemake.output[0] + ".1 "+ snakemake.output[0] + ".2 > "+ snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
        #remove intermediate files: disscarded reads first round, disscarded reads RC, accepted reads first round, accepted reads second round
        subprocess.run( ["rm -f "+ snakemake.params[2]+".tmp "+ snakemake.params[2]+".tmp2 "+snakemake.output[0] + ".1 "+ snakemake.output[0] + ".2"],stdout=subprocess.PIPE, shell=True)
    else:
        #Concatenate results
        subprocess.run( ["cat "+snakemake.output[0] + ".1 "+ snakemake.output[0] + ".2 " + snakemake.params[2] + ".tmp3  > " + snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
        #remove intermediate files: disscarded reads first round, disscarded reads RC, accepted reads first round, accepted reads second round
        #subprocess.run( ["rm -f "+ snakemake.params[2]+".tmp "+ snakemake.params[2]+".tmp2 "+snakemake.output[0] + ".1 "+ snakemake.output[0] + ".2 "+ snakemake.params[2]+".tmp3"],stdout=subprocess.PIPE, shell=True)
else: #no reads to evaluate just rename file
    print("No untrimmed output!!!!")
    subprocess.run( ["mv "+snakemake.output[0] + ".1  > "+ snakemake.output[0]],stdout=subprocess.PIPE, shell=True)

survivingReads=countFasta(snakemake.output[0],False)
prc = float((survivingReads/initialReads)*100)
prc_str = "{:.2f}".format(float((survivingReads/initialReads)*100))

with open(snakemake.params[1], "w") as primers:
        primers.write(primer)
        primers.close()

print("\033[91m This step removes primers \033[0m")
print("\033[93m Total number of initial reads: " + str(initialReads) + " \033[0m")
print("\033[93m Total number of surviving reads: " + str(survivingReads) + " = "+ prc_str + "% \033[0m")
print("\033[93m You can find cutadapt's log file at: " + snakemake.params[5] +"\n \033[0m")
if snakemake.config["interactive"] != "F" or prc < snakemake.config["primers"]["min_prc"]:
    print("\033[92m Do you want to continue?(y/n): \033[0m")
    user_input = stdin.readline() #READS A LINE
    user_input = user_input[:-1]
    if user_input.upper() == "N" or user_input.upper() == "NO":
        subprocess.run( ["rm -f "+ snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
        exit(1)
else:
    print("\033[93m" +" Interactive mode off \033[0m")
    print("\033[93m" +" Removing primers...\033[0m")


if not os.path.exists(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files"):
    os.makedirs(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files")
subprocess.run( ["cat "+ snakemake.output[0]+"| grep '^>' | cut -f1 -d' ' | sed 's/>// ; s/_[0-9]*$//' |  sort | uniq -c | awk '{print $2\"\\t\"$1}' > " + snakemake.params[3]+".tmp1"],stdout=subprocess.PIPE, shell=True)
subprocess.run( ["cat "+ snakemake.input[0]+"| grep '^>' | cut -f1 -d' ' | sed 's/>// ; s/_[0-9]*$//' |  sort | uniq -c | awk '{print $2\"\\t\"$1}'| awk -F'\t' 'NR==FNR{h[$1]=$2;next} BEGIN{print \"Sample\\tReads_before_cutadapt\\tSurviving_reads\\tPrc_surviving_reads\"}{if(h[$1]){print $1\"\\t\"h[$1]\"\\t\"$2\"\\t\"($2/h[$1])*100\"%\"}else{print $1\"\\t\"$2\"\\t0\\t0%\"}}' - "+snakemake.params[3]+".tmp1 > "+ snakemake.params[3]],stdout=subprocess.PIPE, shell=True)
os.remove(snakemake.params[3]+".tmp1")
exit(0)
Python Cutadapt QIIME2.0 VSEARCH benchmark-utils From line 7 of Scripts/remove_adapters_by_sample.py 
  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
import os
import subprocess
from benchmark_utils import countFasta
from sys import stdin

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'Y':'R', 'R':'Y','S':'S','W':'W','K':'M','M':'K','N':'N','B':'V','V':'B','D':'H','H':'D'} 
    bases = list(seq) 
    bases = [complement[base] for base in bases] 
    return ''.join(bases)


def reverse_complement(s):
    return complement(s[::-1])

primer_by_sample={}
uniq_primers={}
idx_fw_primer=-1   # default for qiime (col 3)
idx_rv_primer=-1   # new field 
idx_rv_revcomp_primer=-1
isRC = False  
primer_set = ""
no_primer = ""
extra=snakemake.params[0]
log_by_sample="Sample\tInitial reads\tSurviving reads\n"
if "--discard-untrimmed" in snakemake.params[0]:
    no_primer=" --untrimmed-output " + snakemake.params[2]
    extra=snakemake.params[0].replace("--discard-untrimmed","")

if not os.path.exists(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files"):
    os.makedirs(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files")

if snakemake.config["primers"]["remove"].lower() == "metadata":
    with open(snakemake.input[1]) as mappingFile:
        l=0
        for line in mappingFile:
            l=l+1;
            columns = line.split('\t')
            #the header is always at row 1 and must contain these first 3 fields (qiime specs):
            #SampleID BarcodeSequence LinkerPrimerSequence Description
            if l==1 :
                c=0
                #Find target headers
                for col in columns:
                    if col == "ReversePrimer" or col == "LinkerPrimerSequenceReverse"  or col == "ReverseLinkerPrimerSequence"  or col == "RvLinkerPrimerSequence" or col == "ReversePrimerSequence" :
                        idx_rv_primer=c
                    elif col == "LinkerPrimerSequence":
                        idx_fw_primer=c
                    elif col == "ReverseLinkerPrimerSequenceRevCom"  or col == "ReversePrimerRevCom":
                        idx_rv_revcomp_primer=c
                        isRC=True
                    c=c+1
                if isRC:
                    idx_rv_primer=idx_rv_revcomp_primer 
            elif not line.startswith("#"):
                if idx_rv_primer != -1:
                    #here, we are creating a dic with sample:primer
                    if isRC:  
                        primer_by_sample[columns[0]]=[columns[idx_fw_primer],columns[idx_rv_primer]]
                    else:
                        primer_by_sample[columns[0]]=[columns[idx_fw_primer],reverse_complement(columns[idx_rv_primer])]
                    #for primer in uniq_primers:
                    if columns[idx_fw_primer]+columns[idx_rv_primer] not in uniq_primers:
                        if isRC:
                            uniq_primers[columns[idx_fw_primer]+columns[idx_rv_primer]]=[columns[idx_fw_primer],columns[idx_rv_primer]]
                        else:
                            uniq_primers[columns[idx_fw_primer]+columns[idx_rv_primer]]=[columns[idx_fw_primer],reverse_complement(columns[idx_rv_primer])]    
                else:
                    primer_by_sample[columns[0]]=[columns[idx_fw_primer]]
                    if columns[idx_fw_primer] not in uniq_primers:
                        uniq_primers[columns[idx_fw_primer]]=[columns[idx_fw_primer]]
        mappingFile.close()

    #If we have more than one different pair of primers, we run cutadapt by sample
    #otherwise we run only one instance
    if len(uniq_primers) >1:
        #create tmp dir
        if not os.path.exists(snakemake.params[4]):
            os.makedirs(snakemake.params[4])
        else: #it exists and most lickly we want to delete all its content.
            subprocess.run( ["rm -fr " + snakemake.params[4]+"*"],stdout=subprocess.PIPE, shell=True)
        #split the reads
        #If we are running this, it comes from our demultiplexing, and thus we have fasta headers like this:
        #><sample>_###  so we remove the _###
        subprocess.run(["cat "+ snakemake.input[0]+ " |  awk '{if($0 ~ \"^>\"){sample=$1; header=$0; gsub(\">\",\"\",sample);gsub(\"_[0-9].*\",\"\",sample);}else{print header\"\\n\"$0 >> \""+snakemake.params[4]+"\"sample\".fasta\"} }'"],stdout=subprocess.PIPE, shell=True)
        all_primers=""
        for file in os.listdir(snakemake.params[4]):
            #file only has the name of the file, the path is already discarded
            #the function os.path.splitext strip the extension
            sample=os.path.splitext(file)[0]
            no_primer=""
            extra=""
            if "--discard-untrimmed" in snakemake.params[0]:
                no_primer=" --untrimmed-output " + snakemake.params[4]+sample+"_untrimmed.fasta"
                extra=snakemake.params[0].replace("--discard-untrimmed","")
            tmp_out = snakemake.params[4]+sample+"_trimmed.fasta"
            tmp_log = snakemake.params[4]+sample+".log"
            if sample in primer_by_sample:
                if len(primer_by_sample[sample])>1:
                    primer_set=" -g "+primer_by_sample[sample][0]+"..."+primer_by_sample[sample][1]+" "
                else:
                    primer_set=" -g "+primer_by_sample[sample][0]
                #run cutadapt by sample
                subprocess.run(["echo \"Processing sample\" " + sample + "\n >> "+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
                subprocess.run( ["cutadapt  "+ primer_set +" "+extra+" -o "+tmp_out + " "+ no_primer +" " + snakemake.params[4]+file+ ">>"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
                #stats by sample
                initialReads=countFasta(snakemake.params[4]+file,False)
                survivingReads=countFasta(tmp_out,False)
                prc = "{:.2f}".format(float((survivingReads/initialReads)*100))
                log_by_sample=log_by_sample+sample+"\t"+str(initialReads)+"\t"+str(survivingReads)+" ("+prc+"%)\n"
                all_primers=all_primers+sample+"\t"+primer_set+"\n"
            else:
                print("\033[91mNo primers found for sample:"+ sample+ " \033[0m")
                print("\033[91mPlease make sure to have the sample included in the mapping file: "+snakemake.input[1]+"  \033[0m")
                print("\033[91mAborting the pipeline \033[0m")
                exit(1)

        #merge results
        subprocess.run( ["cat  "+ snakemake.params[4]+"*_trimmed.fasta > "+ snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
        with open(snakemake.params[1], "a") as primers:
            primers.write(all_primers)
            primers.close()
        #subprocess.run( ["cat  "+ snakemake.params[4]+"*_untrimmed.fasta > " snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
    else: #only run one cutadapt instance
        new_key = list(uniq_primers)
        if len(uniq_primers[new_key[0]])>1: #is PE?
            primer_set=" -g "+uniq_primers[new_key[0]][0]+"..."+uniq_primers[new_key[0]][1]+" "
        else: #is SE
            primer_set=" -g "+uniq_primers[new_key[0]][0]
        subprocess.run( ["cutadapt  "+ primer_set +" "+extra+" -o "+snakemake.output[0] + " "+ no_primer +" " + snakemake.input[0]+ ">"+  snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
        with open(snakemake.params[1], "a") as primers:
            primers.write(primer_set)
            primers.close()
else: #values come at the CFG, run only once
    primer_set="-g " + snakemake.config["primers"]["fw_primer"]
    if snakemake.config["primers"]["rv_primer"].len() > 2 :
        primer_set=primer_set+"..."+reverse_complement(snakemake.config["primers"]["rv_primer"])

    subprocess.run( ["cutadapt  "+ primer_set  +" "+extra+" -o "+snakemake.output[0] + " "+ no_primer +" " + snakemake.input[0]+ ">"+ snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
  #  primer_set = snakemake.config["cutadapt"]["adapters"]
    with open(snakemake.params[1], "w") as primers:
        primers.write(primer_set)
        primers.close()

initialReads=countFasta(snakemake.input[0],False)
survivingReads=countFasta(snakemake.output[0],False)
prc=float((survivingReads/initialReads)*100)
prc_str = "{:.2f}".format(float((survivingReads/initialReads)*100))

user_input="0"
while (user_input != "1" and user_input !=  "2"):
    print("\033[91m This step removes primers \033[0m")
    print("\033[93m Total number of initial reads: " + str(initialReads) + " \033[0m")
    print("\033[93m Total number of surviving reads: " + str(survivingReads) + " = "+ prc_str + "% \033[0m")
    print("\033[93m You can find cutadapt's log file at: " + snakemake.params[5] +"\n \033[0m")
    if snakemake.config["interactive"] != "F" or prc < snakemake.config["primers"]["min_prc"]:
        print("\033[92m What would you like to do? \033[0m")
        print("\033[92m 1. Continue with the workflow. \033[0m")
        print("\033[92m 2. Interrupt the workflow. \033[0m")
        if snakemake.config["primers"]["remove"].lower() == "metadata" and  len(uniq_primers)>1:
            print("\033[92m 3. Print results by sample. \033[0m")
        user_input = stdin.readline() #READS A LINE
        user_input = user_input[:-1]
        if user_input == "2":
            print("\033[91m Aborting workflow... \033[0m")
            #delete target outpu (snakemake also does it)
            subprocess.run( ["rm -f "+ snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
            #delete cutadapt mp directory
            subprocess.run( ["rm -fr " + snakemake.params[4]],stdout=subprocess.PIPE, shell=True)
            #delete all the concatenated log files
            subprocess.run( ["rm -f " + snakemake.params[5]],stdout=subprocess.PIPE, shell=True)
            #delete primers file
            subprocess.run( ["rm -f " + snakemake.params[1]],stdout=subprocess.PIPE, shell=True)
            exit(1)
        if user_input == "3":
            print(log_by_sample)
    else:
        print("\033[93m" +" Interactive mode off \033[0m")
        print("\033[93m" +" Removing primers...\033[0m")
        user_input="1"
# if we ran multiple cutadap tasks, now delete tmp files and logs. 
if snakemake.config["primers"]["remove"].lower() == "metadata" and  len(uniq_primers)>1: 
    print("\033[96mCleaning intermediate files...\033[0m")
    subprocess.run( ["rm -fr " + snakemake.params[4]],stdout=subprocess.PIPE, shell=True)

#Summarize results    
subprocess.run( ["cat "+ snakemake.output[0]+"| grep '^>' | cut -f1 -d' ' | sed 's/>// ; s/_[0-9]*$//' |  sort | uniq -c | awk '{print $2\"\\t\"$1}' > " + snakemake.params[3]+".tmp1"],stdout=subprocess.PIPE, shell=True)
subprocess.run( ["cat "+ snakemake.input[0]+"| grep '^>' | cut -f1 -d' ' | sed 's/>// ; s/_[0-9]*$//' |  sort | uniq -c | awk '{print $2\"\\t\"$1}'| awk -F'\t' 'NR==FNR{h[$1]=$2;next} BEGIN{print \"Sample\\tReads_before_cutadapt\\tSurviving_reads\\tPrc_surviving_reads\"}{if(h[$1]){print $1\"\\t\"h[$1]\"\\t\"$2\"\\t\"($2/h[$1])*100\"%\"}else{print $1\"\\t\"$2\"\\t0\\t0%\"}}' - "+snakemake.params[3]+".tmp1 > "+ snakemake.params[3]],stdout=subprocess.PIPE, shell=True)
os.remove(snakemake.params[3]+".tmp1")
exit(0)
Python Cutadapt QIIME2.0 benchmark-utils From line 7 of Scripts/remove_adapters_v2.py 
  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
import os
import subprocess
from sys import stdin
from benchmark_utils import countFasta
from benchmark_utils import countFastaGZ
import shutil

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'Y':'R', 'R':'Y','S':'S','W':'W','K':'M','M':'K','N':'N','B':'V','V':'B','D':'H','H':'D'} 
    bases = list(seq) 
    bases = [complement[base] for base in bases] 
    return ''.join(bases)


def reverse_complement(s):
    return complement(s[::-1])

# List files
fq_files = [f for f in os.listdir(snakemake.params[0]) if f.endswith("_1."+snakemake.params[2])]
if not os.path.exists(snakemake.params[0]+"/reads_discarded_primer/") and "--discard-untrimmed" in snakemake.params[1]:
    os.makedirs(snakemake.params[0]+"/reads_discarded_primer/")
if not os.path.exists(snakemake.params[0]+"/primer_removed/"):
    os.makedirs(snakemake.params[0]+"/primer_removed/")
if not os.path.exists(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files"):
    os.makedirs(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files")

summ_file = open(snakemake.output[0],"w") # this iss a log for the wf
summ_file2 = open(snakemake.params[4],"w") # this is for the report
summ_file.write("Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n");
summ_file2.write("Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n");
log_str = "Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n"
log_zero = "Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n"
has_zero_length_reads = False
zero_samples = 0;
to_remove = []

for fw in fq_files:
    sample=fw.replace("_1."+snakemake.params[2],"")
    fw_fq= snakemake.params[0]+"/"+fw
    rv=fw.replace("_1."+snakemake.params[2],"_2."+snakemake.params[2])
    rv_fq= snakemake.params[0]+"/"+rv
    discard_untrimmed=""
    extra_params=snakemake.params[1]
#Count reads before trimming
    if snakemake.params[2].endswith("gz"):
        reads_ori=countFastaGZ(fw_fq,True)
    else:
        reads_ori=countFasta(fw_fq,True)
#no cutadapt if no reads 
    if reads_ori > 0:
        if snakemake.params[3] == "PE":
            if "--discard-untrimmed" in snakemake.params[1]:
                discard_untrimmed=" --untrimmed-output "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_1.fastq.gz --untrimmed-paired-output  "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_2.fastq.gz"
                extra_params=snakemake.params[1].replace("--discard-untrimmed","")
            subprocess.run(["cutadapt -g "+ snakemake.config["primers"]["fw_primer"]  + " -G " + snakemake.config["primers"]["rv_primer"]  + " " +extra_params+" -O "+ snakemake.config["primers"]["min_overlap"]+" -m "+ snakemake.config["primers"]["min_length"] +" -o "+snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz -p "+snakemake.params[0]+"/primer_removed/"+sample+"_2.fastq.gz "+discard_untrimmed +" "+ fw_fq + " " +  rv_fq + " >> "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log"],stdout=subprocess.PIPE, shell=True)
        elif snakemake.params[3] == "SE":
            if "--discard-untrimmed" in snakemake.params[0]:
                discard_untrimmed=" --untrimmed-output "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_1.fastq.gz"
                extra_params=snakemake.params[1].replace("--discard-untrimmed","") 
            subprocess.run(["cutadapt -g "+ snakemake.config["primers"]["fw_primer"] +" " +extra_params+" -O "+ snakemake.config["primers"]["min_overlap"]+" -m "+ snakemake.config["primers"]["min_length"] +" -o "+snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz "+ discard_untrimmed + " " + fw_fq + " >> "+ snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log"],stdout=subprocess.PIPE, shell=True)  

        if snakemake.params[2].endswith("gz"):
            reads_after=countFastaGZ(snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz",True)
        else:
            reads_after=countFasta(snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq",True)

        prcOK="{:.2f}".format(float((reads_after/reads_ori)*100))

    else:
        reads_after = 0
        prcOK="{:.2f}".format(float((reads_after/1)*100))
        to_copy=snakemake.params[0]+"/primer_removed/"+sample+"_1."+snakemake.params[2]
        os.symlink(fw_fq,to_copy)
        if snakemake.params[3] == "PE":
            to_copy_rv=snakemake.params[0]+"/primer_removed/"+sample+"_2."+snakemake.params[2]
            os.symlink(rv_fq,to_copy_rv)

    if reads_after < 1:
        has_zero_length_reads = True
        log_zero = log_zero + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
        to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_1."+snakemake.params[2])
        if snakemake.params[3] == "PE":
            to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_2."+snakemake.params[2])
        zero_samples = zero_samples + 1

    log_str = log_str + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
    summ_file.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");
    summ_file2.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");

summ_file.close()
summ_file2.close()

user_input="0"
show_menu = True
if zero_samples > 0:
    while show_menu:
        print("\033[91m\n###########  Primer removal validation    ###########\033[0m")
        print("\033[91m You have " + str(zero_samples) + " samples without reads surviving filters. \033[0m")
        print("\033[92m LIBRARY: "+snakemake.wildcards.sample+" \033[0m")
        print("\033[92m cutadapt_log: "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log \033[0m")
        print("\033[93m Please select one of the following options: \033[0m")
        print("\033[93m   1. Print samples with 0 reads \033[0m")
        print("\033[93m   2. Print summary (all the samples) \033[0m")
        print("\033[93m   3. Remove from this analysis samples with 0 reads\033[0m")
        print("\033[93m      and continue with the workflow. \033[0m")
        print("\033[93m   4. Interrupt the workflow and re-do primer removal step. \033[0m")
        print("\033[93m      Adjust primer values in your configuration and/or mapping file \033[0m")
        print("\033[93m      and restart the pipeline. \033[0m")
        print("\033[93m      This action will remove:"+snakemake.params[0]+"/primer_removed \033[0m")
        print("\033[93m   5. Interrupt the workflow \033[0m")
        print("\033[93m   6. Continue with the workflow\n      (an error will be raised during dada2)\n      Pointless option... \033[0m")
        print("\033[93m Select an option: \033[0m")
        user_input = stdin.readline() #READS A LINE
        user_input = user_input[:-1]
        if user_input == "1":
            print(log_zero)
        elif user_input == "2":
            print(log_str)
        elif user_input == "3":
            for file in to_remove:
                newn = file+"_NOK"
                os.rename(file, newn)
                show_menu = False
        elif user_input == "4":
            shutil.rmtree(snakemake.params[0]+"/primer_removed")
            exit(1) 
        elif user_input == "5":
            exit(1)

exit(0)
Python Cutadapt benchmark-utils From line 7 of Scripts/removePrimersDemultiplex_cfg.py 
  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
203
204
205
206
207
208
209
210
211
212
213
import os
import subprocess
from benchmark_utils import countFasta
from benchmark_utils import countFastaGZ
from sys import stdin
import shutil

def complement(seq):
    complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'Y':'R', 'R':'Y','S':'S','W':'W','K':'M','M':'K','N':'N','B':'V','V':'B','D':'H','H':'D'} 
    bases = list(seq) 
    bases = [complement[base] for base in bases] 
    return ''.join(bases)


def reverse_complement(s):
    return complement(s[::-1])

#from Bio.Seq import Seq


primer_by_sample={}
uniq_primers={}
idx_fw_primer=-1   # default for qiime (col 3)
idx_rv_primer=-1  # new field 
idx_rv_revcomp_primer=-1
isRC = False  
with open(snakemake.input[0]) as mappingFile:
    l=0
    for line in mappingFile:
        l=l+1;
        columns = line.split('\t')
        #the header is always at row 1 and must contain these first 3 fields (qiime specs):
        #SampleID BarcodeSequence LinkerPrimerSequence Description
        if l==1 :
            c=0
            #Find target headers
            for col in columns:
                if col == "ReverseLinkerPrimerSequence"  or col == "RvLinkerPrimerSequence" or col == "ReversePrimer" or col == "ReversePrimerSequence"  :
                    idx_rv_primer=c
                elif col == "LinkerPrimerSequence":
                    idx_fw_primer=c
                elif col == "ReverseLinkerPrimerSequenceRevCom"  or col == "RvLinkerPrimerSequenceRevCom" or col == "ReversePrimerRevCom":
                    idx_rv_revcomp_primer=c                   
                c=c+1
            #if there is no "ReverseLinkerPrimerSequence" we look for the ReverseLinkerPrimerSequenceRevCom
            if idx_rv_primer == -1 and idx_rv_revcomp_primer !=1:
                idx_rv_primer=idx_rv_revcomp_primer
                isRC = True 
        elif not line.startswith("#"):
            if idx_rv_primer != -1:
                #if the valuee is the reverse complemented now we want the 5' to 3' orientation so rev com again.
                if isRC:
                    primer_by_sample[columns[0]]=[columns[idx_fw_primer],reverse_complement(columns[idx_rv_primer])]
                else:
                    primer_by_sample[columns[0]]=[columns[idx_fw_primer],columns[idx_rv_primer]]
            elif idx_fw_primer != -1:
                primer_by_sample[columns[0]]=[columns[idx_fw_primer]]
            else:
                print("\033[91m ERROR: LinkerPrimerSequence not found on mapping file: "+ snakemake.input[0] +" \033[0m")
                exit(1)


# List files
fq_files = [f for f in os.listdir(snakemake.params[0]) if f.endswith("_1."+snakemake.params[2])]
if not os.path.exists(snakemake.params[0]+"/reads_discarded_primer/"):
    os.makedirs(snakemake.params[0]+"/reads_discarded_primer/")
if not os.path.exists(snakemake.params[0]+"/primer_removed/"):
    os.makedirs(snakemake.params[0]+"/primer_removed/")
if not os.path.exists(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files"):
    os.makedirs(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files")
summ_file = open(snakemake.output[0],"w")
summ_file2 = open(snakemake.params[4],"w")
summ_file.write("Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n")
summ_file2.write("Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n")
log_str = "Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n"
log_zero = "Sample\tReads_before_cutadapt\tSurviving_reads\tPrc_surviving_reads\n"
has_zero_length_reads = False
zero_samples = 0;
to_remove = []

for fw in fq_files:
    sample=fw.replace("_1."+snakemake.params[2],"")
    fw_fq= snakemake.params[0]+"/"+fw
    rv=fw.replace("_1."+snakemake.params[2],"_2."+snakemake.params[2])
    rv_fq= snakemake.params[0]+"/"+rv
    #Count reads before trimming
    if snakemake.params[2].endswith("gz"):
        reads_ori=countFastaGZ(fw_fq,True)
    else:
        reads_ori=countFasta(fw_fq,True)
    #no cutadapt if no reads
    #if reads_ori > 0:

    if sample in primer_by_sample:
        runCutAdapt=False
        discard_untrimmed=""
        extra_params=snakemake.params[1] 
        if len(primer_by_sample[sample])>1 and reads_ori < 1:
            reads_after = 0
            prcOK="{:.2f}".format(float((reads_after/1)*100))
            to_copy=snakemake.params[0]+"/primer_removed/"+sample+"_1."+snakemake.params[2]
            os.symlink(fw_fq,to_copy)
            if snakemake.params[3] == "PE":
                to_copy_rv=snakemake.params[0]+"/primer_removed/"+sample+"_2."+snakemake.params[2]
                os.symlink(rv_fq,to_copy_rv)
            runCutAdapt=True

        elif len(primer_by_sample[sample])>1 and snakemake.params[3] == "PE" :
            if "--discard-untrimmed" in snakemake.params[1]:
                discard_untrimmed=" --untrimmed-output "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_1.fastq.gz --untrimmed-paired-output  "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_2.fastq.gz"
                extra_params=snakemake.params[1].replace("--discard-untrimmed","")
            #print("cutadapt -g "+ primer_by_sample[sample][0] + " -G " + primer_by_sample[sample][1] + " " +extra_params+" -O "+ snakemake.config["primers"]["min_overlap"] +" -o "+snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz -p "+snakemake.params[0]+"/primer_removed/"+sample+"_2.fastq.gz "+discard_untrimmed +" "+ fw_fq + " " +  rv_fq + " >> "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log")
            subprocess.run(["cutadapt -g "+ primer_by_sample[sample][0] + " -G " + primer_by_sample[sample][1]+" -m "+ snakemake.config["primers"]["min_length"]  + " " +extra_params+" -O "+ snakemake.config["primers"]["min_overlap"] +" -o "+snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz -p "+snakemake.params[0]+"/primer_removed/"+sample+"_2.fastq.gz "+discard_untrimmed +" "+ fw_fq + " " +  rv_fq + " >> "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log"],stdout=subprocess.PIPE, shell=True)
            runCutAdapt=True
            #subprocess.run(["grep \"(passing filters)\" "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log | awk '{print \""+sample+"\t\"$5\"\t\"$6}' >> "+snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
            #subprocess.run( ["cutadapt "+ primer_set +" "+snakemake.params[0]+" -o "+snakemake.output[0] + " " + snakemake.input[0]+ ">"+ snakemake.output[1]],stdout=subprocess.PIPE, shell=True)
        elif len(primer_by_sample[sample])>=1 and snakemake.params[3] == "SE":
            if "--discard-untrimmed" in snakemake.params[0]:
                discard_untrimmed=" --untrimmed-output "+snakemake.params[0]+"/reads_discarded_primer/"+sample+"_1.fastq.gz"
                extra_params=snakemake.params[1].replace("--discard-untrimmed","") 
            subprocess.run(["cutadapt -g "+ primer_by_sample[sample][0] +" -m "+ snakemake.config["primers"]["min_length"] +" " +extra_params+" -O "+ snakemake.config["primers"]["min_overlap"] +" -o "+snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz "+ discard_untrimmed + " " + fw_fq + " >> "+ snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log"],stdout=subprocess.PIPE, shell=True)  
            #subprocess.run(["grep \"(passing filters)\" "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log | awk '{print \""+sample+"\t\"$5\"\t\"$6}' >> "+snakemake.output[0]],stdout=subprocess.PIPE, shell=True)
            runCutAdapt=True
        elif len(primer_by_sample[sample])==1 and snakemake.params[3] == "PE":
            print("\033[91m ERROR: Found forward and reverse reads, but only one primer was supplied \033[0m")
            print("sample: "+sample + " primer " + primer_by_sample[sample][0])
            summ_file.close()
            summ_file2.close()
            exit(1)

        if runCutAdapt and reads_ori > 0:
            if snakemake.params[2].endswith("gz"):
                reads_ori=countFastaGZ(fw_fq,True)
                reads_after=countFastaGZ(snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq.gz",True)
            else:
                reads_ori=countFasta(fw_fq,True)
                reads_after=countFasta(snakemake.params[0]+"/primer_removed/"+sample+"_1.fastq",True)
            prcOK="{:.2f}".format(float((reads_after/reads_ori)*100))
            summ_file.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");
            summ_file2.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");
            log_str = log_str + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
            if reads_after < 1:
                has_zero_length_reads = True
                log_zero = log_zero + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
                to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_1."+snakemake.params[2])
                if snakemake.params[3] == "PE":
                    to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_2."+snakemake.params[2])
                zero_samples = zero_samples + 1


        elif runCutAdapt and reads_ori < 1:
            summ_file.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");
            summ_file2.write(sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n");
            log_str = log_str + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
            log_zero = log_zero + sample+"\t"+str(reads_ori)+"\t"+str(reads_after)+"\t"+prcOK+"\n"
            zero_samples = zero_samples + 1
            has_zero_length_reads = True
            to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_1."+snakemake.params[2])
            if snakemake.params[3] == "PE":
                to_remove.append(snakemake.params[0]+"/primer_removed/"+sample+"_2."+snakemake.params[2])


    else:
        print("\033[93m WARNING: No primers found for sample: "+sample +" \033[0m")
        summ_file.close()
        summ_file2.close()
        exit(1)
summ_file.close() 
summ_file2.close()

user_input="0"
show_menu = True
if zero_samples > 0:
    while show_menu:
        print("\033[91m\n###########  Primer removal validation    ###########\033[0m")
        print("\033[91m You have " + str(zero_samples) + " samples without reads surviving filters. \033[0m")
        print("\033[92m LIBRARY: "+snakemake.wildcards.sample+" \033[0m")
        print("\033[92m cutadapt_log: "+snakemake.params[0]+"/primer_removed/"+sample+".cutadapt.log \033[0m")
        print("\033[93m Please select one of the following options: \033[0m")
        print("\033[93m   1. Print samples with 0 reads \033[0m")
        print("\033[93m   2. Print summary (all the samples) \033[0m")
        print("\033[93m   3. Remove from this analysis samples with 0 reads\033[0m")
        print("\033[93m      and continue with the workflow. \033[0m")
        print("\033[93m   4. Interrupt the workflow and re-do primer removal step. \033[0m")
        print("\033[93m      Adjust primer values in your configuration and/or mapping file \033[0m")
        print("\033[93m      and restart the pipeline. \033[0m")
        print("\033[93m      This action will remove:"+snakemake.params[0]+"/primer_removed \033[0m")
        print("\033[93m   5. Interrupt the workflow \033[0m")
        print("\033[93m Select an option: \033[0m")
        user_input = stdin.readline() #READS A LINE
        user_input = user_input[:-1]
        if user_input == "1":
            print(log_zero)
        elif user_input == "2":
            print(log_str)
        elif user_input == "3":
            for file in to_remove:
                newn = file+"_NOK"
                os.rename(file, newn)
                show_menu = False
        elif user_input == "4":
            shutil.rmtree(snakemake.params[0]+"/primer_removed")
            exit(1)
        elif user_input == "5":
            exit(1)

exit(0)
Python Cutadapt QIIME2.0 benchmark-utils From line 7 of Scripts/removePrimersDemultiplex.py 
  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
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
import subprocess
from snakemake.utils import report
import benchmark_utils 
from benchmark_utils import countTxt
from benchmark_utils import readBenchmark
from benchmark_utils import countFasta
from benchmark_utils import countFastaGZ
from benchmark_utils import readSampleDist
from benchmark_utils import make_table
from countData import parseCounts
from seqsChart import createChart

#Parse the total number of counts
#countTxt = parseCounts(snakemake.input.counts)

################################################################################
#                         TOOLS VERSION SECTION                          #
################################################################################
#--fastq
fqv = subprocess.run([snakemake.config["fastQC"]["command"], '--version'], stdout=subprocess.PIPE)
fqVersion = "**" + fqv.stdout.decode('utf-8').strip() + "**"

if snakemake.config["demultiplexing"]["demultiplex"] !=  "F":
   #--qiime extract_barcodes
   ebv = subprocess.run([snakemake.config["qiime"]["path"]+'extract_barcodes.py', '--version'], stdout=subprocess.PIPE)
   ebVersion = ebv.stdout.decode('utf-8')
   ebVersion = "**" + ebVersion[ebVersion.find(":")+1:].strip() + "**"
   #--qiime split_libraries
   spVersion = "**TBD**"
   spv = subprocess.run([snakemake.config["qiime"]["path"]+'split_libraries_fastq.py', '--version'], stdout=subprocess.PIPE)
   spVersion = spv.stdout.decode('utf-8')
   if "Version" in spVersion:
       spVersion = "**" + spVersion[spVersion.find(":")+1:].strip() + "**"
else:
   ebVersion = "**NA**"
   SPvERSION = "**NA**"

vsearchVersion = "**TBD**"
vsearchV = subprocess.run(['vsearch', '--version'], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
vsearchVersion = "**" + vsearchV.stdout.decode('utf-8').split('\n', 1)[0].strip() + "**"


#--qiime identify_chimeric_seqs
icVersion = "**TBD**"
icv = subprocess.run([snakemake.config["qiime"]["path"]+'identify_chimeric_seqs.py', '--version'], stdout=subprocess.PIPE)
icVersion = icv.stdout.decode('utf-8')
if "Version" in icVersion:
    icVersion = "**" + icVersion[icVersion.find(":")+1:].strip() + "**"
#--pear
try:
    pearv = subprocess.run( [snakemake.config["pear"]["command"]+" -h | grep 'PEAR v'"], stdout=subprocess.PIPE, shell=True)
    pearversion = "**" + pearv.stdout.decode('utf-8').strip() + "**"
except Exception as e:
    pearversion = "Problem reading version"

#--cutadapt
cutVersion = "**TBD**"
if snakemake.config["primers"]["remove"].casefold() == "metadata" or snakemake.config["primers"]["remove"].casefold() == "cfg" or snakemake.config["primers"]["remove"].lower() != "f":
    cutv = subprocess.run(['cutadapt', '--version'], stdout=subprocess.PIPE)
    cutVersion = "**cutadapt v" + cutv.stdout.decode('utf-8').strip() + "**"
    #cutVersion = "cutadapt v TBD"

################################################################################
#                          Chimera check                                       #
################################################################################
removeChimeras = False
if snakemake.config["chimera"]["search"] == "T":
    ################################################################################
    #                       Read log file from remove_chimera.py                   #
    # After search for chimera, user have the option to remove them or not. If the #
    # user decides to remove the chimera, the executed command is stored on the log#
    # file, otherwise it stores a message indicating the user decision.            #
    ################################################################################
    chimera_log = ""
    try:
        with open(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/chimera/chimera.log") as chimlog:
            for line in chimlog:
                chimera_log += line
            chimlog.close()
    except FileNotFoundError:
        chiemra_log = "No Log for identify_chimeric_seqs.py"
    if "The chimeric sequences were removed" in chimera_log:
        removeChimeras = True


################################################################################
#                         Benchmark Section                                    #
# This section is to generate a pre-formatted text with the benchmark info for #
# All the executed rules.                                                      #
################################################################################
fqBench = readBenchmark(snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/qc/fq.benchmark")
pearBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/pear.benchmark")
if snakemake.config["demultiplexing"]["demultiplex"] != "F":
    barBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.benchmark")
    splitLibsBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/splitLibs.benchmark")
    #splitLibsRCBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/splitLibs.benchmark")
   # combineBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/combine_seqs_fw_rev.benchmark")
else:
    combineBench=pearBench #THIS IS ONLY FOR TESTING REMOVE!!! 
rmShorLongBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/filter.benchmark")
demultiplexFQBench=""
if snakemake.config["demultiplexing"]["demultiplex"] == "T" and snakemake.config["demultiplexing"]["create_fastq_files"] == "T":
    demultiplexFQBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/demultiplex_fq.benchmark")

################################################################################
#                           Compute Counts                                     #
################################################################################
if snakemake.config["gzip_input"] == "F":
    rawCounts = countFasta(snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/rawdata/fw.fastq", True);
else:
    rawCounts = countFastaGZ(snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/rawdata/fw.fastq.gz", True);
#rawCountsStr= '{0:g}'.format(float(rawCounts))
rawCountsStr= str(int(rawCounts))
#-peared
pearedCounts = 0
if snakemake.config["UNPAIRED_DATA_PIPELINE"] != "T":
    pearedCounts = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.fastq", True);
else:
    pearedCounts = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.UNPAIRED.fastq", True);

#pearedCountsStr='{0:g}'.format(float(pearedCounts))
pearedCountsStr=str(int(pearedCounts))
prcPeared = "{:.2f}".format(float((pearedCounts/rawCounts)*100))
#-dumultiplex
if snakemake.config["demultiplexing"]["demultiplex"] != "F": #starting to test this  and snakemake.config["demultiplexing"]["bc_mismatch"]>0:
    #in the past we had two files fw and reverse nos everything is on one file
    #fwAssignedCounts = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.fna", False)
    #barcodes.fastq_corrected_toRC
    #rvAssignedCounts = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/seqs.assigned.fna", False)

    #prcFwAssigned = "{:.2f}".format(float((fwAssignedCounts/pearedCounts)*100))
    #prcRvAssigned = "{:.2f}".format(float((rvAssignedCounts/pearedCounts)*100))
    #totalAssigned = fwAssignedCounts + rvAssignedCounts
    #prcPearedAssigned = "{:.2f}".format(float((totalAssigned/pearedCounts)*100))
    #prcRawAssigned = "{:.2f}".format(float((totalAssigned/rawCounts)*100))
    #New implementation
    totalAssigned =  countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.fna", False)
    rvAssignedCounts = countTxt(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq_corrected_toRC")
    fwAssignedCounts = totalAssigned - rvAssignedCounts 
    prcFwAssigned = "{:.2f}".format(float((fwAssignedCounts/pearedCounts)*100))
    prcRvAssigned = "{:.2f}".format(float((rvAssignedCounts/pearedCounts)*100))
    prcPearedAssigned = "{:.2f}".format(float((totalAssigned/pearedCounts)*100))
    prcRawAssigned = "{:.2f}".format(float((totalAssigned/rawCounts)*100))

else: 
    totalAssigned = pearedCounts
    prcPearedAssigned = "{:.2f}".format(float((totalAssigned/pearedCounts)*100))
    prcRawAssigned = "{:.2f}".format(float((totalAssigned/rawCounts)*100))

#--cutadapt
cutSequences = False
if snakemake.config["primers"]["remove"].casefold() == "metadata" or snakemake.config["primers"]["remove"].casefold() == "cfg":
    sequenceNoAdapters = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna", False)
    if (totalAssigned - sequenceNoAdapters) > 0:
        cutSequences = True
        prcCut = "{:.2f}".format(float((sequenceNoAdapters/totalAssigned)*100))
        prcCutRaw = "{:.2f}".format(float((sequenceNoAdapters/rawCounts)*100))

if removeChimeras:
    sequenceNoChimeras = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered_nc.fasta", False)
    prcChim = "{:.2f}".format(float((sequenceNoChimeras/totalAssigned)*100))
    prcChimRaw = "{:.2f}".format(float((sequenceNoChimeras/rawCounts)*100))
    if cutSequences:
        prcChimCut = "{:.2f}".format(float((sequenceNoChimeras/sequenceNoAdapters)*100))
#out="{PROJECT}/runs/{run}/{sample}_data/"
trimmedCounts = countFasta(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered.fasta", False)
prcTrimmedSplit ="{:.2f}".format(float((trimmedCounts/totalAssigned)*100))
prcTrimmedRaw= "{:.2f}".format(float((trimmedCounts/rawCounts)*100))
if cutSequences:
    prcTrimmedCut="{:.2f}".format(float((trimmedCounts/sequenceNoAdapters)*100))
#if removeChimeras:
#    prcTrimmedChimera="{:.2f}".format(float((trimmedCounts/sequenceNoChimeras)*100))
try:
    samplesLib = subprocess.run( ["cat " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered.dist.txt | wc -l"], stdout=subprocess.PIPE, shell=True)
    samplesLibInt = int(samplesLib.stdout.decode('utf-8').strip())
except Exception as e:
    totalReads = "Problem reading outputfile"
################################################################################
#                         Generate sequence amounts chart                      #
################################################################################
numbers=[rawCounts,pearedCounts];
labels=["Raw", "Assembled"];
if snakemake.config["demultiplexing"]["demultiplex"] == "T":
    numbers.append(totalAssigned)
    labels.append("Demultiplexed")
if snakemake.config["primers"]["remove"].casefold() == "metadata" or snakemake.config["primers"]["remove"].casefold() == "cfg":
    numbers.append(sequenceNoAdapters)
    labels.append("Cutadapt")
numbers.append(trimmedCounts)
labels.append("Length filtering")
if snakemake.config["chimera"]["search"] == "T" and removeChimeras:
    numbers.append(sequenceNoChimeras)
    labels.append("No Chimera")
createChart(numbers, tuple(labels),snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files/sequence_numbers."+snakemake.wildcards.sample+".png")
################################################################################
#                          Chimera check                                       #
################################################################################
variable_refs=""
if snakemake.config["chimera"]["search"] == "T" and snakemake.config["chimera"]["method"] == "usearch61":
    variable_refs+= ".. [usearch61] Edgar RC. 2010. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26(19):2460-2461.\n\n"
else: 
    variable_refs+= ".. [uchime] Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R (2011) UCHIME improves sensitivity and speed of chimera detection. Bioinformatics, 27 (16): 2194-2200. doi:10.1093/bioinformatics/btr381. \n\n"
quimeraStr = ""
if snakemake.config["chimera"]["search"] == "T":
    quimeraStr="Identify Chimera\n-------------------\n\n"
    quimeraStr+="Identify possible chimeric sequences (sequences generated due to the PCR amplification of multiple templates or parent sequences).\n\n"
    if snakemake.config["chimera"]["method"] == "usearch61":
        quimeraStr += ":red:`Tool:` [QIIME]_ - identify_chimeric_seqs.py\n\n"
        quimeraStr += ":red:`Version:` "+ icVersion +"\n\n"
        quimeraStr += ":red:`Method:` [usearch61]_ \n\n"
    else:
        quimeraStr += ":red:`Tool:` [Vsearch]_ - vsearch\n\n"
        quimeraStr += ":red:`Version:` "+ vsearchVersion +"\n\n"        
        quimeraStr += ":red:`Method:` "+ str(snakemake.config["chimera"]["method"]) +" - uses [uchime]_ \n\n"    
    quimeraStr += "**Command:**\n\n"
    if snakemake.config["chimera"]["method"] == "usearch61":
         quimeraStr+=":commd:`identify_chimeric_seqs.py -m "+ str(snakemake.config["chimera"]["method"])+" -i "+ str(snakemake.wildcards.PROJECT)+"/runs/"+str(snakemake.wildcards.run)+"/"+str(snakemake.wildcards.sample)+"_data/seqs_fw_rev_accepted.fna "+str(snakemake.config["chimera"]["extra_params"])
         quimeraStr+=" -o "+ str(snakemake.wildcards.PROJECT)+"/runs/"+str(snakemake.wildcards.run)+"/"+str(snakemake.wildcards.sample)+"_data/chimera/` \n\n"
    else:
         quimeraStr+=":commd:`vsearch --"+ str(snakemake.config["chimera"]["method"])+" "+ str(snakemake.wildcards.PROJECT)+"/runs/"+str(snakemake.wildcards.run)+"/"+str(snakemake.wildcards.sample)+"_data/seqs_fw_rev_accepted.fna --threads "+ str(snakemake.config["chimera"]["threads"]) +" " +str(snakemake.config["chimera"]["extra_params"])   
         quimeraStr+=" --uchimeout "+ str(snakemake.wildcards.PROJECT)+"/runs/"+str(snakemake.wildcards.run)+"/"+str(snakemake.wildcards.sample)+"_data/chimera/chimeras.summary.txt` \n\n"
    quimeraStr+="**Output files:**\n\n"
    if snakemake.config["chimera"]["method"] == "usearch61":
        quimeraStr+=":green:`- File with the possible chimeric sequences:` "+str(snakemake.wildcards.PROJECT)+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/chimera/chimeras.txt\n\n"
    else:
        quimeraStr+=":green:`- File with the possible chimeric sequences:` "+str(snakemake.wildcards.PROJECT)+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/chimera/chimeras.summary.txt\n\n"
    identifyChimeraBench=readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/chimera/chimera.benchmark")
    quimeraStr+=identifyChimeraBench
    quimeraStr+=chimera_log
    if removeChimeras:
        quimeraStr+=":red:`Reads after remove chimeric sequences:` "+ str(sequenceNoChimeras)+"\n\n"
        quimeraStr+=":red:`Percentage of reads vs raw reads:` "+ str(prcChimRaw) + "%\n\n"
        quimeraStr+=":red:`Percentage of reads vs demultiplexed reads:` "+ str(prcChim) + "%\n\n"
        if cutSequences:
            quimeraStr+=":red:`Percentage of reads vs cutadapt:` "+ str(prcChimRaw) + "%\n\n"




################################################################################
#                           Peared FastQC                                     #
################################################################################
fastQCPearStr = ""
if snakemake.config["fastQCPear"] == "T":
    fastQCPearStr = "Peared FastQC Analysis\n------------------------\n\n" # title
    fastQCPearStr += "Check the quality of the reads after assembly.\n\n"
    fastQCPearStr += ":red:`Tool:` [FastQC]_\n\n"
    fastQCPearStr += ":red:`Version:` "+ fqVersion +"\n\n"
    fastQCPearStr += "**Command:**\n\n"
    fastQCPearStr += ":commd:`fastqc "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.fastq --extract -o  "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/qc`\n\n"
    fastQCPearStr += "**Output files:**\n\n:green:`- FastQC report:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/qc/seqs.assembled_fastqc.html FQ_Report_ \n\n"
    fastQCPearStr += ".. _FQ_Report: peared/qc/seqs.assembled_fastqc.html \n\n"
    fastQCPearStrBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/qc/fq.benchmark")
    fastQCPearStr += fastQCPearStrBench

################################################################################
#                           Extract Barcode                                    #
################################################################################
extractBCStr = ""
if snakemake.config["demultiplexing"]["demultiplex"] != "F":
    extractBCStr ="Extract barcodes\n-----------------\n\n"
    extractBCStr +="Extract the barcodes used to identify individual samples.\n\n"
    extractBCStr +=":red:`Tool:` [QIIME]_ - extract_barcodes.py\n\n"
    extractBCStr +=":red:`Version:` "+ebVersion+"\n\n"
    extractBCStr +="**Command:**\n\n"
    extractBCStr +=":commd:`extract_barcodes.py -f "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.fastq -c "+str(snakemake.config["ext_bc"]["c"])+ " " + str(snakemake.config["ext_bc"]["bc_length"])+ " " + snakemake.config["ext_bc"]["extra_params"] + " -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/`\n\n"
    extractBCStr +="**Output files:**\n\n"
    extractBCStr +=":green:`- Fastq file with barcodes:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq\n\n"
    extractBCStr +=":green:`- Fastq file with the reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/reads.fastq\n\n"
    extractBCStr +=barBench
################################################################################
#                           CORRECT Barcodes                                   #
################################################################################
correctBCStr = ""
bcFile="barcodes.fastq"
if snakemake.config["demultiplexing"]["demultiplex"] != "F": # and snakemake.config["demultiplexing"]["bc_mismatch"]:
    correctBCStr = "Correct Barcodes\n--------------------\n"
    correctBCStr += "Try to correct the barcode from unassigned reads and place reads in correct orientetion.\n\n"
    correctBCStr += "Maximum number of mismatches **"  + str(snakemake.config["demultiplexing"]["bc_mismatch"]) + "**.\n\n"
    correctBCStr +=":red:`Tool:` Cascabel Java application\n\n"
    correctBCStr +="**Command:**\n\n"
    correctBCStr += ":commd:`java -jar Scripts/BarcodeCorrector.jar -b "+snakemake.wildcards.PROJECT+"/metadata/sampleList_mergedBarcodes_"+snakemake.wildcards.sample+".txt -fb "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq -fr "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/reads.fastq  -m "  + str(snakemake.config["demultiplexing"]["bc_mismatch"]) + " -o  " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq_corrected -or  " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/reads.fastq_corrected -rc -x " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/sample_matrix.txt  >  " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/demux.log`\n\n"
    correctBCStr += "**Output files:**\n\n:green:`- Barcode corrected file:` "+snakemake.wildcards.PROJECT+ "/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq_corrected\n\n"
    correctBCStr += ":green:`- Reads corrected file:` "+snakemake.wildcards.PROJECT+ "/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/reads.fastq_corrected\n\n"
    correctBCStr += ":green:`- Error correction summary:` "+snakemake.wildcards.PROJECT+ "/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/demux.log\n\n"
    correctBarBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes_corrected.benchmark")
    correctBCStr += correctBarBench
    bcFile="barcodes.fastq_corrected"

splitStr = ""
if snakemake.config["demultiplexing"]["demultiplex"] != "F":
    splitStr+="Demultiplexing\n"
    splitStr+="----------------\n"
    splitStr+="For library splitting, also known as demultiplexing, Cascabel performs several steps to assign fragments in the original as well as reverse orientation to the correct sample.\n\n"
    splitStr+="Split samples from Fastq file\n"
    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
    splitStr+=":red:`Tool:` [QIIME]_ - split_libraries_fastq.py\n\n"
    splitStr+=":red:`version:` "+ spVersion+"\n\n"
    splitStr+="**Command:**\n\n"
    splitStr+=":commd:`split_libraries_fastq.py -m "+snakemake.wildcards.PROJECT+"/metadata/sampleList_mergedBarcodes_"+snakemake.wildcards.sample+".txt -i "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/reads.fastq -o  "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs -b "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/"+bcFile+" -q "+str(snakemake.config["split"]["q"])+" -r "+str(snakemake.config["split"]["r"])+" --retain_unassigned_reads "+str(snakemake.config["split"]["extra_params"])+" --barcode_type "+str(snakemake.config["split"]["barcode_type"])+"`\n\n"
    splitStr+=splitLibsBench

    splitStr+="Retain assigned reads\n"
    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
    splitStr+="**Command:**\n\n"
    splitStr+=":commd:`cat "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.fna | grep -P -A1 \"(?!>Unass)^>\" | sed '/^--$/d' > "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.fna`\n\n"

    splitStr+="Create file with only unassigned reads\n"
    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
    splitStr+="**Command:**\n\n"
    splitStr+=":commd:`cat "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.fna | grep \"^>Unassigned\" |  sed 's/>Unassigned_[0-9]* /@/g' | sed 's/ .*//' | grep -F -w -A3  -f - "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.fastq |  sed '/^--$/d' >"+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/unassigned.fastq`\n\n"

#    splitStr+="Reverse complement unassigned reads\n"
#    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
#    splitStr+=":red:`Tool:` [Vsearch]_\n\n"
#    splitStr+=":red:`version:`  "+vsearchVersion+"\n\n"
#    splitStr+="**Command:**\n\n"
#    splitStr+=":commd:`vsearch --fastx_revcomp "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/unassigned.fastq  --fastqout "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/unassigned.reversed.fastq`\n\n"


#    splitStr+="Barcode extraction for reverse complemented, unassigned reads\n"
#    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
#    splitStr +=":red:`Tool:` [QIIME]_ - extract_barcodes.py\n\n"
#    splitStr +=":red:`Version:` "+ebVersion+"\n\n"
#    splitStr+="**Command:**\n\n"
#    splitStr +=":commd:`extract_barcodes.py -f "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/unassigned.reversed.fastq -c "+str(snakemake.config["ext_bc"]["c"])+" "+str(snakemake.config["ext_bc"]["bc_length"])+" "+snakemake.config["ext_bc"]["extra_params"]+" -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes_unassigned/`\n\n"

#    if snakemake.config["bc_mismatch"]:
#        splitStr += "Correct reverse complemented barcodes \n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
#        splitStr += "Maximum number of mismatches **"  + str(snakemake.config["bc_mismatch"]) + "**.\n\n"
#        splitStr +=":red:`Tool:` Cascabel Java application\n\n"
#        splitStr +="**Command:**\n\n"
#        splitStr += ":commd:`java -cp Scripts/BarcodeCorrector/build/classes/  barcodecorrector.BarcodeCorrector -b "+snakemake.wildcards.PROJECT+"/metadata/sampleList_mergedBarcodes_"+snakemake.wildcards.sample+".txt -f "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes_unassigned/barcodes.fastq_corrected -m "  + str(snakemake.config["bc_mismatch"]) + "`\n\n"
#        splitStr += "**Output file:**\n\n:green:`- Barcode corrected file:` "+snakemake.wildcards.PROJECT+ "/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes/barcodes.fastq_corrected\n\n"
#        splitStrBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes_unassigned/barcodes_corrected.benchmark")
#        splitStr += splitStrBench+"\n\n"

#    splitStr +="Split reverse complemented reads\n"
#    splitStr+="~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n"
#    splitStr +=":red:`Tool:` [QIIME]_ - extract_barcodes.py\n\n"
#    splitStr +=":red:`Version:` "+ebVersion+"\n\n"
#    splitStr+="**Command:**\n\n"
#    splitStr +=":commd:`split_libraries_fastq.py -m "+snakemake.wildcards.PROJECT+"/metadata/sampleList_mergedBarcodes_"+snakemake.wildcards.sample+".txt -i "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/barcodes_unassigned/reads.fastq -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC -b "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+str(snakemake.wildcards.sample)+"_data/barcodes_unassigned/"+bcFile+" -q "+str(snakemake.config["split"]["q"])+" -r "+str(snakemake.config["split"]["r"])+" "+str(snakemake.config["split"]["extra_params"])+" --barcode_type "+str(snakemake.config["split"]["barcode_type"])+"`\n\n"
#    splitStr +=splitLibsBench+"\n\n"

    splitStr +="**Output files:**\n\n"
#   # splitStr +=":green:`- FW reads fasta file with new header:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.fna\n\n"
    splitStr +=":green:`- Text histogram with the length of the fw reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/histograms.txt\n\n"
    splitStr +=":green:`- Log file for the fw reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/split_library_log.txt\n\n"
#   # splitStr +=":green:`- RV reads fasta file with new header:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/seqs.assigned.fna\n\n"
#    splitStr +=":green:`- Text histogram with the length of the rv reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/histograms.txt\n\n"
#    splitStr +=":green:`- Log file for the rv reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/split_library_log.txt\n\n"
#    splitStr +=":green:`- Fasta file with unassigned reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/seqs.unassigned.fna\n\n"
    splitStr +=":red:`Number of reads assigned on FW:` "+str(fwAssignedCounts)+" = "+str(prcFwAssigned)+"% of the peared reads\n\n"
    splitStr +=":red:`Number of reads assigned on RVC:` "+str(rvAssignedCounts)+" = "+str(prcRvAssigned)+"% of the peared reads\n\n"

################################################################################
#                           Single FastQ creation                              #
################################################################################
demultiplexFQ = ""
if snakemake.config["demultiplexing"]["demultiplex"] == "T" and snakemake.config["demultiplexing"]["create_fastq_files"] == "T":
    demultiplexFQ = "Generate single sample fastq files\n------------------------------------------\n\n" # title
    demultiplexFQ += "Create single fastq files per samples (based on the raw data without applying any filtering).\n\n"
    demultiplexFQ +=":red:`Tool:` Cascabel Java program\n\n"
    demultiplexFQ += "**Command:**\n\n"
    demultiplexFQ += ":commd:`"+snakemake.config["java"]["command"] + " -cp Scripts DemultiplexQiime --txt -a rv -b "+ str(snakemake.config["demultiplexing"]["bc_mismatch"]) + " -d "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.ori.txt -o "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/ "
    ext=".gz"
    if snakemake.config["gzip_input"].casefold() == "f":
        ext=""
    demultiplexFQ += "-r1 "+snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/rawdata/fw.fastq"+ext+" -r2 "+snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/rawdata/fw.fastq"+ext+"`\n\n"
    if snakemake.config["demultiplexing"]["remove_bc"]:
        demultiplexFQ +=":red:`Barcodes removed:` "+ str(snakemake.config["demultiplexing"]["remove_bc"]) + " first bases\n\n"
#Now only for ASV workflow
   # if snakemake.config["primers"]["remove"].lower() == "cfg":
   #     demultiplexFQ +=":red:`Primers removed:` **FW** " + snakemake.config["primers"]["fw_primer"] + " **RV** " +snakemake.config["primers"]["rv_primer"]+"\n\n"
   # elif snakemake.config["primers"]["remove"].lower() == "metadata":
   #     demultiplexFQ +=":red:`Removed primers` were obtained from the metadata file.\n\n" 
    demultiplexFQ += "**The demultiplexed fastq files are located at:**\n\n:green:`- Demultiplexed directory:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/\n\n"
    demultiplexFQ += ":green:`- Summary file:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/summary.txt\n\n"
    demultiplexFQ += demultiplexFQBench
   # also this only for the ASV workflow
   # if (snakemake.config["primers"]["remove"].lower() == "cfg" or snakemake.config["primers"]["remove"].lower() == "metadata"):
   #     demultiplexFQ += "**Remove primers:**\n\nFollowing, primers were removed from the fastq files\n\n"
   #     demultiplexFQ +=":red:`Tool:` [Cutadapt]_\n\n"
   #     demultiplexFQ += ":red:`Version:` "+cutVersion+"\n\n"
   #     demultiplexFQ += "**Command:**\n\n"
   #     if snakemake.config["primers"]["remove"].lower() == "cfg":
   #         if snakemake.config["LIBRARY_LAYOUT"].casefold()=="pe":
   #             demultiplexFQ += ":commd:`cutadapt -g "+ snakemake.config["primers"]["fw_primer"]  + " -G " + snakemake.config["primers"]["rv_primer"]  + " " +snakemake.config["primers"]["extra_params"]+" -O "+ snakemake.config["primers"]["min_overlap"]  +" -m " +snakemake.config["primers"]["min_length"]+ " -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_1.fastq.gz -p "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_2.fastq.gz "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_1.fq.gz  "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_2.fq.gz  >> "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/"+snakemake.wildcards.sample+".cutadapt.log`\n\n"
   #         else:
   #             demultiplexFQ += ":commd:`cutadapt -g "+ snakemake.config["primers"]["fw_primer"]  + " " +snakemake.config["primers"]["extra_params"]+" -O "+ snakemake.config["primers"]["min_overlap"]  +" -m " +snakemake.config["primers"]["min_length"]+ " -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_1.fastq.gz "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_1.fq.gz  >> "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/"+snakemake.wildcards.sample+".cutadapt.log`\n\n"
   #         demultiplexFQ += "The above command ran once for each single sample fastq file(s) using the mentioned primers\n\n"
   #     else: #is from metadata
   #         if snakemake.config["LIBRARY_LAYOUT"].casefold()=="pe":
   #             demultiplexFQ += ":commd:`cutadapt -g sample_FW_primer  -G sample_RV_primer " +snakemake.config["primers"]["extra_params"]+" -O "+ snakemake.config["primers"]["min_overlap"]  +" -m " +snakemake.config["primers"]["min_length"]+ " -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_1.fastq.gz -p "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_2.fastq.gz "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_1.fq.gz "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_2.fq.gz  >> "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/"+snakemake.wildcards.sample+".cutadapt.log`\n\n"
   #         elif snakemake.config["LIBRARY_LAYOUT"].casefold()=="se":
   #             demultiplexFQ += ":commd:`cutadapt -g sample_FW_primer "+ " " +snakemake.config["primers"]["extra_params"]+" -O "+ snakemake.config["primers"]["min_overlap"]  +" -m " +snakemake.config["primers"]["min_length"]+ " -o "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/SAMPLE_1.fastq.gz "+ snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/SAMPLE_1.fq.gz  >> "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed/"+snakemake.wildcards.sample+".cutadapt.log`\n\n"
   #         demultiplexFQ += "The above command ran once for each single sample fastq file(s) and primers were obtained from the mapping file accordingly to its sample\n\n"    
   #     demultiplexFQ += ":green:`- Reads without primers:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/primer_removed\n\n"
   #     if "--discard-untrimmed" in snakemake.config["primers"]["extra_params"]:
   #         demultiplexFQ += ":green:`- Discarded reads (no primer):` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/demultiplexed/reads_discarded_primer\n\n"
   #     else:
   #         demultiplexFQ += ":red:`- Given the options, reads without primers where not removed!`\n\n"
   #     demultiplexFQ += ":green:`- Primer removal results by sample:` primers_removal_\n\n"
   #     demultiplexFQ +=" .. _primers_removal: report_files/cutadapt."+snakemake.wildcards.sample+".fastq_summary.tsv\n\n"

################################################################################
#                           Combine FW and Reverse reads                       #
################################################################################

combineFR = ""
#if snakemake.config["demultiplexing"]["demultiplex"] != "F":
#    combineFR = "Combine reads\n---------------------------------\n\n" # title
#    combineFR += "Concatenate forward and reverse reads.\n\n"
#    combineFR += "**Command:**\n\n"
#    combineFR += ":commd:`cat "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibs/seqs.assigned.fna "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/splitLibsRC/seqs.assigned.fna > "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna`\n\n"
#    combineFR +="**Output files:**\n\n"
#    combineFR +=":green:`- Fasta file with combined reads:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna\n\n"
#    combineFR +=":red:`- Total number of acepted reads:` " +str(totalAssigned)+ " = "+ str(prcPearedAssigned)+ "% of the peared reads or "+str(prcRawAssigned)+"% of the raw reads.\n\n"
#    combineFR += combineBench

################################################################################
#                          Cut adapters                                        #
################################################################################
cutAdaptStr = ""
if snakemake.config["primers"]["remove"].casefold() == "metadata" or snakemake.config["primers"]["remove"].casefold() == "cfg":
    cutAdaptStr = "Remove sequence primers\n------------------------\n\n" # title
    cutAdaptStr +="Remove the adapters / primers from the reads.\n\n"
    cutAdaptStr +=":red:`Tool:` [Cutadapt]_\n\n"
    cutAdaptStr += ":red:`Version:` "+cutVersion+"\n\n"
    cutAdaptStr += "**Command:**\n\n"
    primer_lines=0
    if snakemake.config["primers"]["remove"].lower() == "cfg":
        #cutAdaptStr += ":commd:`cutadapt "+ str(snakemake.config["cutadapt"]["adapters"])+" " + str(snakemake.config["cutadapt"]["extra_params"]) + " -o " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna\n\n"
        #cutAdaptStr +=  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna > " +  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.log`\n\n"
        cutAdaptStr += ":commd:`cutadapt -g "+ str(snakemake.config["primers"]["fw_primer"])+"..."+str(snakemake.config["primers"]["rv_primer"])+" "+ str(snakemake.config["primers"]["extra_params"]) + " -o " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna > " +  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.log`\n\n"
    elif snakemake.config["primers"]["remove"].lower() == "metadata":
        primers=""
        try:
            #with open(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/primers.txt") as pfile:
            with open(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/report_files/primers."+snakemake.wildcards.sample+".txt") as pfile:
                primers=pfile.read()
                #primer_lines=len(pfile.readlines())
                primer_lines=len(primers.split("\n"))
                if primer_lines > 1:
                    if snakemake.config["LIBRARY_LAYOUT"].casefold()=="pe":
                        primers="-g sample_FW_primer...sampleRV_primer"
                    else:
                        primers="-g sample_FW_primer"

        except FileNotFoundError:
            primers="-ERROR reading primer file-"
        #cutAdaptStr += ":commd:`cutadapt "+primers +" " + str(snakemake.config["cutadapt"]["extra_params"]) + " -o " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna\n\n"
        #cutAdaptStr += snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna > " +  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.log`\n\n"
        cutAdaptStr += ":commd:`cutadapt "+primers +" " + str(snakemake.config["primers"]["extra_params"]) + " -o " + snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna "+  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna > " +  snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.log`\n\n"

#cutAdaptStr += "*PRIMERS: primer sequences were obtained from the metadata file\n\n"
    if primer_lines > 1:
        cutAdaptStr += ":green:`- Primers used by sample:` primers_sample_\n\n"
        cutAdaptStr +=  ".. _primers_sample: report_files/primers."+snakemake.wildcards.sample+".txt\n\n"
    cutAdaptStr += "**Output files:**\n\n:green:`- Reads without adapters:` "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna\n\n"
    if cutSequences:
        cutAdaptStr += ":red:`Total number of reads after cutadapt:` "+ str(sequenceNoAdapters) + " = " + str(prcCut) + "% of the assigned reads or "+ str(prcCutRaw)+"% of the total reads\n\n"
    #cutAdaptStr+=":\n\n"
    cutAdaptStr+=":green:`- Primer removal results by sample:` primers_OTU_\n\n"
    cutAdaptStr+=" .. _primers_OTU: report_files/cutadapt."+snakemake.wildcards.sample+".summary.tsv\n\n"

    cutAdaptBench =readBenchmark(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/cutadapt.benchmark")
    cutAdaptStr += cutAdaptBench+"\n\n"
################################################################################
#                          Counts for too long too shorts                      #
################################################################################
#trimmedStr =  ":red:`Total number of reads after trimming:` "+str(trimmedCounts)+ "="+ str(prcTrimmedSplit)+"% of the demultiplexed reads or " + str(prcTrimmedRaw) + "% of the raw reads\n\n"
trimmedStr =  ":red:`Total number of reads after length filtering:` "+str(trimmedCounts)+ "\n\n"
trimmedStr += ":red:`Percentage of reads vs raw reads:` "+str(prcTrimmedRaw)+"%\n\n"
trimmedStr+=":red:`Percentage of reads vs demultiplexed reads:` " + str(prcTrimmedSplit) + "%\n\n"
if cutSequences:
    trimmedStr+=":red:`Percentage of reads after cutadapt:` "+ str(prcTrimmedCut) + "%\n"
#if removeChimeras:
#    trimmedStr+=":red:`Percentage of reads after remove chimeras vs trimmed reads:` "+ str(prcTrimmedChimera) + "%\n"


#bcValidationBench =readBenchmark(snakemake.wildcards.PROJECT+"/metadata/bc_validation/"+snakemake.wildcards.sample+"/validation.benchmark")
################################################################################
#                     Remove too short and too long reads                      #
#  This rule creates a temporary file with the short and long values choosed   #
#  by the user in order to remove the reads. The file filter.log contains the  #
#  minimun expected length for the reads followed by the maximun length tab    #
#  separated (shorts <TAB> longs)                                              #
################################################################################
shorts = str(snakemake.config["rm_reads"]["shorts"])
longs = str(snakemake.config["rm_reads"]["longs"])
with open(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/filter.log") as trimlog:
    for line in trimlog:
        tokens = line.split("\t")
        if len(tokens)>2:
            shorts = tokens[1]
            longs = tokens[2]
################################################################################
#                       FInal Counts                              #
################################################################################
countTxt="Following you can see the final read counts: \n\n"
fileData = []
headers = []
data =[]
headers.append("File description")
headers.append("Location")
headers.append("Number of reads")
headers.append("Prc(%) vs raw")
fileData.append(headers)
#raw
data.append("Raw reads")
data.append(snakemake.wildcards.PROJECT+"/samples/"+snakemake.wildcards.sample+"/rawdata/\*.fq")
data.append(str(rawCounts))
data.append("{:.2f}".format(float((rawCounts/rawCounts)*100))+"%")
fileData.append(data)
data=[]
#pear
data.append("Assembled reads")
data.append(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/peared/seqs.assembled.fastq")
data.append(str(pearedCounts))
data.append("{:.2f}".format(float((pearedCounts/rawCounts)*100))+"%")
fileData.append(data)
data=[]
#splitted
if snakemake.config["demultiplexing"]["demultiplex"] == "T":
    data.append("Demultiplexed reads")
    data.append(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted.fna")
    data.append(str(totalAssigned))
    data.append("{:.2f}".format(float((totalAssigned/rawCounts)*100))+"%")
    fileData.append(data)
    data=[]
#adapters
if snakemake.config["primers"]["remove"].casefold() == "metadata" or snakemake.config["primers"]["remove"].casefold() == "cfg":
    data.append("Adapter removed")
    data.append(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_accepted_no_adapters.fna")
    data.append(str(sequenceNoAdapters))
    data.append("{:.2f}".format(float((sequenceNoAdapters/rawCounts)*100))+"%")
    fileData.append(data)
    data=[]
#length filtered
data.append("Length filtered")
data.append(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered.fasta")
data.append(str(trimmedCounts))
data.append("{:.2f}".format(float((trimmedCounts/rawCounts)*100))+"%")
fileData.append(data)
data=[]
#chimera
if snakemake.config["chimera"]["search"] == "T" and removeChimeras:
    data.append("Non chimeric reads")
    data.append(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered_nc.fasta")
    data.append(str(sequenceNoChimeras))
    data.append("{:.2f}".format(float((sequenceNoChimeras/rawCounts)*100))+"%")
    fileData.append(data)
    data=[]
countTxt += make_table(fileData)
################################################################################
#                       Sample distribution chart                              #
################################################################################

sampleDistChart = ""
if snakemake.config["demultiplexing"]["demultiplex"] == "T":
    dist_table = readSampleDist(snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered.dist.txt",trimmedCounts,samplesLibInt)
    sampleDistChart = "Sample distribution\n--------------------------------------\n\n" # title
    sampleDistChart += dist_table + "\n\n"
    sampleDistChart += ".. image:: report_files/seqs_fw_rev_filtered."+snakemake.wildcards.sample+".dist.png\n\n"
    sampleDistChart +="The previous chart shows the number of clean reads per sample. The bars are sorted from left to right, according to the metadata input file.\n\n"
    sampleDistChart +="**To see more details about the number of reads per sample in this library, please refer to the file:** "+snakemake.wildcards.PROJECT+"/runs/"+snakemake.wildcards.run+"/"+snakemake.wildcards.sample+"_data/seqs_fw_rev_filtered.dist.txt\n\n"


################################################################################
#                       User description section                               #
################################################################################
desc = snakemake.config["description"]
txtDescription = ""
if len(desc) > 0:
    txtDescription = "\n**User description:** "+desc+"\n"


################################################################################
#                       controls warning section                               #
################################################################################
"""
We want to include a small section to warn the user about the use of controls. This could be
the case if they are demultiplexing a complete library. 
"""
ctrlWarning =""
if snakemake.config["demultiplexing"]["demultiplex"] == "T":
    ctrlWarning="\n:warn:`Note: Library demultiplexing has been carried out, if you have controls among your samples, please be aware that Cascabel won't perform any special operation with them. They are treated as any other sample within this workflow. Please make sure to analyze your controls with other tools, and correct your sample counts for potential contamination.`\n"
################################################################################
#                                Report                                        #
################################################################################

report("""
Amplicon Analysis Report for Library: {snakemake.wildcards.sample}
=====================================================================
    .. role:: commd
    .. role:: red
    .. role:: green
    .. role:: warn

**CASCABEL** is designed to run amplicon sequence analysis across single or multiple read libraries.

The objective of this pipeline is to create different output files which allow the user to explore data in a simple and meaningful way, as well as facilitate downstream analysis, based on the generated output files.

Another aim of **CASCABEL** is also to encourage the documentation process, by creating this report in order to assure data analysis reproducibility.

{txtDescription}

{ctrlWarning}

Following you can see all the steps that were taken in order to get the final results of the pipeline.

Raw Data
---------
The raw data for this library can be found at:

:green:`- FW raw reads:` {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/fw.fastq

:green:`- RV raw reads:` {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/rv.fastq

:red:`Number of total reads:` {rawCountsStr}

Quality Control
------------------
Evaluate quality on raw reads.

:red:`Tool:` [FastQC]_

:red:`Version:` {fqVersion}

**Command:**

:commd:`fastqc {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/fw.fastq {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/rv.fastq --extract -o {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/qc/`

You can follow the links below, in order to see the complete FastQC report:

:green:`- FastQC for sample {snakemake.wildcards.sample}_1:` FQ1_

    .. _FQ1: ../../../samples/{snakemake.wildcards.sample}/qc/fw_fastqc.html

:green:`- FastQC for sample {snakemake.wildcards.sample}_2:` FQ2_

    .. _FQ2: ../../../samples/{snakemake.wildcards.sample}/qc/rv_fastqc.html

{fqBench}


Read pairing
----------------
Align paired end reads and merge them into one single sequence in case they overlap.

:red:`Tool:` [PEAR]_

:red:`version:` {pearversion}

**Command:**

:commd:`pear -f {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/fw.fastq -r {snakemake.wildcards.PROJECT}/samples/{snakemake.wildcards.sample}/rawdata/rv.fastq -t {snakemake.config[pear][t]} -v {snakemake.config[pear][v]} -j {snakemake.config[pear][j]} -p {snakemake.config[pear][p]} -o {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/peared/seqs > {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/peared/seqs.assembled.fastq`

**Output files:**

:green:`- Merged reads:` {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/peared/seqs.assembled.fastq

:green:`- Log file:` {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/peared/pear.log

:red:`Number of peared reads:` {pearedCountsStr} =  {prcPeared}%

{pearBench}

{fastQCPearStr}

{extractBCStr}

{correctBCStr}

{splitStr}

{demultiplexFQ}

{combineFR}

{cutAdaptStr}


Remove too long and too short reads
------------------------------------
Remove very short and long reads, with lengths more than some standard deviation below or above the mean to be short or long respectively

:green:`- Minimun length expected (shorts):` {shorts}

:green:`- Maximun length expected (longs):` {longs}

**Command:**

:commd:`awk '!/^>/ {{ next }} {{ getline seq }} length(seq) > shorts  && length(seq) < longs {{ print $0 \"\\n\" seq }}'  {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/seqs_fw_rev_accepted.fna  >  {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/seqs_fw_rev_filtered.fasta`

**Sequence distribution before remove reads**

.. image:: report_files/seqs_dist_hist.{snakemake.wildcards.sample}.png
    :height: 400px
    :width: 400px
    :align: center


**Output file:**

:green:`- Fasta file with correct sequence length:` {snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/{snakemake.wildcards.sample}_data/seqs_fw_rev_filtered.fasta

{trimmedStr}

{rmShorLongBench}


{quimeraStr}


{sampleDistChart}


Final counts
-------------

{countTxt}

.. image:: report_files/sequence_numbers.{snakemake.wildcards.sample}.png

OTU report
---------------------------

Cascabel report on downstream analyses in combination with multiple libraries (if supplied), can be found at the following link: otu_report_ ({snakemake.wildcards.PROJECT}/runs/{snakemake.wildcards.run}/otu_report_{snakemake.config[assignTaxonomy][tool]}.html)

    .. _otu_report: otu_report_{snakemake.config[assignTaxonomy][tool]}.html

References
------------------

.. [FastQC] FastQC v0.11.3. Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data

.. [PEAR] PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Zhang et al (2014) Bioinformatics 30(5): 614-620 | doi:10.1093/bioinformatics/btt593

.. [QIIME] QIIME. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, Fierer N, Gonzalez Pena A, Goodrich JK, Gordon JI, Huttley GA, Kelley ST, Knights D, Koenig JE, Ley RE, Lozupone CA, McDonald D, Muegge BD, Pirrung M, Reeder J, Sevinsky JR, Turnbaugh PJ, Walters WA, Widmann J, Yatsunenko T, Zaneveld J, Knight R. 2010. QIIME allows analysis of high-throughput community sequencing data. Nature Methods 7(5): 335-336.

.. [Cutadapt] Cutadapt v1.15 .Marcel Martin. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.Journal, 17(1):10-12, May 2011. http://dx.doi.org/10.14806/ej.17.1.200

.. [Vsearch] Rognes T, Flouri T, Nichols B, Quince C, Mahé F. (2016) VSEARCH: a versatile open source tool for metagenomics. PeerJ 4:e2584. doi: 10.7717/peerj.2584


{variable_refs}


""", snakemake.output[0], metadata="Author: J. Engelmann & A. Abdala ")
Python Cutadapt QIIME2.0 VSEARCH PEAR benchmark-utils From line 1 of Scripts/report_v2.py

A Snakemake pipeline for variant calling of genomic FASTQ data using GATK 

34
35
36
37
38
39
40
41
42
43
44
45
shell: """
    # ADAPTERS=$(gunzip -c {input.reads[0]} | head -n 1 | cut -d " " -f 2 | cut -d ":" -f 4);
    # i7=$(echo $ADAPTERS | cut -d "+" -f 1);
    # i5=$(echo $ADAPTERS | cut -d "+" -f 2);
    R1_END_ADAPTER="{params.pre_i7}";  # ${{i7}}{params.post_i7}";
    R2_END_ADAPTER="{params.pre_i5}";  # ${{i5}}{params.post_i5}";
    cutadapt {input.reads} \
        -a $R1_END_ADAPTER \
        -A $R2_END_ADAPTER \
        --cores {threads} \
        -o {output.trimmed[0]} \
        -p {output.trimmed[1]}"""
SnakeMake Cutadapt From line 34 of rules/variant_calling.smk
tool / biotools

Cutadapt

cutadapt

Find and remove adapter sequences, primers, poly-A tails and other types of unwanted sequence from your high-throughput sequencing reads.

biotools cutadapt
email [email protected]
edam data_3495
edam format_1930
edam format_1929
URL https://galaxy.pasteur.fr/tool_runner?tool_id=toolshed.pasteur.fr/repos/mvalade/cutadapt_v2/cutadapt/1.6
doi 10.7490/f1000research.1114334.1
edam operation_3192
URL https://pypi.python.org/pypi/cutadapt
URL https://cutadapt.readthedocs.io/en/stable/
doi 10.14806/ej.17.1.200
topic: Sequencing
toolType: Command-line tool
operation: Sequence trimming
signature: cutadapt
input:data: RNA sequence
input:format: FASTA
201 311 0 0
Are you the author of this tool? with a link to this page.
Insert the following code to your website to add tag icon with a link to this page.
    HTML Markdown