####!!!! It is neccessary to have at least 10gb of free space to run this exercise

#### In this exercise, your are going to have a look at our data and perform basic quality control ####

# Please, define variable datafile as the name of one of your fastq files for demonstrative purposes. If you do not have your own datafile, you can download one here:

# Download both, forward and reverse reads of ERR2022829
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR014/SRR014849/SRR014849_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR014/SRR014849/SRR014849_2.fastq.gz

#More information about this sequencing experiment can be found at following links:
# https://www.ebi.ac.uk/ena/browser/view/PRJEB21619
# https://www.ebi.ac.uk/ena/browser/view/ERR2022829

gunzip *.fastq.gz

datafile="SRR014849_1.fastq"
datafile_re="SRR014849_2.fastq"

## Let's have a look at the file header ##
 head $datafile

## Show number of rows ##
 wc -l $datafile

## Calculate number of reads (there are 4 rows per read) ##
 wc -l $datafile | awk '{print $1/4}'

## What kind of coding is used for Phred? ##

 head -n 40 $datafile | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding\!";}'

## Let's get Phred scores and their probabilities ##

# We can extract each 4th row (those with quality scores) to a separate file

 sed -n '4~4p' $datafile > QCs.fastq

# Open R-studio

rstudio

# The following commands are in R console:

# read first 10 lines (first 10 quality scores) from the quality score line

 QCs <- read.table("QCs.fastq", nrows=10, as.is=TRUE)

 dim(QCs)


# Let's have a look at first read

# changing ascii to phred scores (we know it is phred+33, therefore we remove 33)
 phred <- as.integer(charToRaw(QCs[1,]))-33
 phred

# changing phred to probabilites of error
 prob <- 10^(-phred/10)
 prob

# plotting the probabilities
 plot(prob, pch=19, las=1)

# Let's plot qualities of all 10 reads

 plot(prob, pch=19, las=1, type="l", ylim=c(0,1), xlab="Cycle")

 for (i in 2:10)
{
phred <- as.integer(charToRaw(QCs[i,]))-33
prob <- 10^(-phred/10)
lines(prob, las=1, col=i)
}

 legend("topleft", legend=paste("read", c(1:10)), col=1:10, pch=19)


# Based on the graph, would you say that these reads are of good quality?


### Pre-alignment quality control ###

# back in terminal
## Using FastQC ##

# running the program in interactive version

 fastqc

# running the program in command line, output in the FastQC directory
mkdir FastQC

fastqc -o FastQC $datafile
fastqc -o FastQC $datafile_re


# If multiqc command is not availible, you can install it in two steps
# 1. Install pip Python 3 package manager by running: sudo apt install python3-pip
# 2. Install multiqc using pip package manager: pip install multiqc
# OR INSTALL LOCALLY on the WOLF computer:
#module add anaconda3
#pip install --user multiqc

cd FastQC

#now call the multiqc
#multiqc .
# on WOLF:
~/.local/bin/multiqc .

# for multiple files, build a for cycle

for file in *.fastq
do
fastqc -o FastQC $file
done


### checking for adapters ###

# next, we will check for the presence of illumina adapters (see https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/experiment-design/illumina-customer-sequence-letter.pdf for the most common illumina adapters)

# searching for TruSeq Universal Adapter:

 grep -c 'AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT' $datafile

# for multiple files, build a for cycle:

for file in *.fastq
do
cmd=`grep -c 'AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT' ${file}`

h1=${cmd}

echo 'AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT'

done


# Please, study the cutadapt options and decipher command below. Adjust selected values for your data if necessary!
module add cutadapt

#cutadapt -q 10 -b 'ATCTCGTATGCCGTCTTCTGCTTG' -m 50 ${datafile} -o ${datafile}-trimmed.fastq > ${datafile}-trimmed.info

# we decided based on the evaluation to make these adjustments:
# - no adapter trimming
# - shortening the minimal length to 10
# - raising the quality to 20

cutadapt -q 20 -m 10 ${datafile} -o ${datafile}-trimmed.fastq > ${datafile}-trimmed.info

# creating set of commands and saving to .sh for a number of *fastq files and calling it
for file in *.fastq
do
fil=`basename ${file} .fastq`
cmd="cutadapt -q 20 -m 10 ${file} -o ${fil}-trimmed.fastq > ${fil}-trimmed.info"
echo ${cmd} >> cutadapt-trim.sh
done

# you have to call the command!
bash cutadapt-trim.sh

# repeat filtering with fastp
#conda install -c bioconda fastp
# or
wget http://opengene.org/fastp/fastp
chmod a+x ./fastp

./fastp --in1 $datafile --in2 $datafile_re --out1 $datafile.trimmed.fq --out2 $datafile_re.trimmed.fq

# now the data are trimmed and filtered, lets do the quality control again on the trimmed data!

fastqc -o FastQC SRR014849_1-trimmed.fastq
fastqc -o FastQC SRR014849_2-trimmed.fastq
cd FastQC

# on all qc files in the folder
~/.local/bin/multiqc .

# only for trimmed 
~/.local/bin/multiqc *trimmed*





