# preprocess

#!/bin/bash

# É necessário passar um diretório de entrada de dados ao chamar o script para execução:

indir=$1

# ! Significa negação, no caso abaixo, se não foi passado diretório de entrada print a mensagem "Missing input directory"

if [ ! ${indir} ]; then
echo "Missing input directory."
exit
fi

# SE o que foi passado não é um diretório, print a mensagem: "Wrong input directory"

if [ !-d ${indir} ]; then
echo "Wrong input directory (${indir})."
exit
fi

# Também é necessário passar um diretório de saída, ou seja onde os arquivos gerados ficarão:

outdir=$2

# SE não foi passado o diretório de saída, print uma mensagem dizendo que o argumento está faltando:

if [ ! ${outdir} ]; then
echo "Missing output directory."
exit
fi

# SE o que foi passado não é um diretório, print uma mensagem dizendo que argumento está errado:
if [ ! -d ${outdir} ]; then
echo "Wrong output directory (${outdir})."
exit
fi

# A função mkdir serve para criação de diretórios,no caso dentro do diretório que foi passado como output será criada uma
# pasta chamada processed, e dentro da mesma, serão criadas outras 3 pastas: atropos, prinseq e fastqc.
# A última conterá outras duas pastas: pre e pos. Ou seja, antes do processamento as reads passarão por um controle de
# qualide(pre) e após o processamento, serão novamente submetidas ao controle de qualidade (pos), pondendo realizar uma
# comparação de ambos posteriormente:

mkdir -p ${outdir}/processed/fastqc/pre
mkdir -p ${outdir}/processed/atropos
mkdir -p ${outdir}/processed/prinseq
mkdir -p ${outdir}/processed/fastqc/pos

# A estrutura for é utilizada para facilitar o processo, com ela pode-se realizar o comando para todas as reads,
# não sendo necessário a execução para cada uma delas. Assim, no comando abaixo, para cada read no formato .fastq será
# executado o mesmo comando: Pegar a base de cada nome (o que difere) e realizar o controle de qualidade (comando fastqc),
# gerando 2 arquivos um de erro (log.err.txt), caso venha a ocorrer, e outro de saída que mostra o andamento do processo
# (log.out.txt).

for r1 in `ls ${indir}/*_R1.fastq`; do

r2=`echo ${r1} | sed 's/_R1.fastq/_R2.fastq/'`

if [ ! -e "${r2}" ]; then
echo "Read 2 (${r2}) paired with Read 1 ($r1) not found."
exit
fi

name=`basename ${r1} | sed 's/_R1.fastq//'`
echo -e "FastQC pre-evaluation using sample ${name}: ${r1} & ${r2} ...\n"

fastqc -t 2 \
${r1} \
-o ${outdir}/processed/fastqc/pre/ \
> ${outdir}/processed/fastqc/pre/${name}_R1.log.out.txt \
2> ${outdir}/processed/fastqc/pre/${name}_R1.log.err.txt

fastqc -t 2 \
${r2} \
-o ${outdir}/processed/fastqc/pre/ \
> ${outdir}/processed/fastqc/pre/${name}_R2.log.out.txt \
2> ${outdir}/processed/fastqc/pre/${name}_R2.log.err.txt

# A função echo é para printar na tela a mensagem que vem em seguida entre aspas. Isso contribui para saber em que etapa
# o processo se encontra.

echo -e "Running atropos (insert) for adapter trimming using sample ${name}: ${r1} & ${r2} ...\n"

# É necessário inicializar o ambiente pyenv dentro desta sessão, pois para a execução do atropos é necessário um
# ambiente Python, dado o seu desenvolvimento

eval "$(pyenv init -)"

# Para ativar o ambiente Python para execução do atropos:
pyenv activate atropos

# Execução do Atropos - trimagem das reads, eliminando os adaptadores e sequencias de baixa qualidade. O atropos funciona
# utilizando 2 algortimos de trimagem o "insert" e o "adapter". O insert realiza a trimagem após alinhamento das reads
# por match das reads, o que sobrar é adaptadore e então ocorre a trimagem. e as reads que não foram trimadas nessa etapa
# passa para o segundo algortimo o adapter, que realiza a trimagem baseando-se na correspondendia do alinhamento dos adaptadores.

# -e: Esse parâmetro fornece o valor máximo da taxa de erro permitida para matches de adaptadores para que o alinhamento
# seja considerado que obteve sucesso. É o número de erros dividido pelo comprimento da região de matching.

# -n: Nessa opção deve-se optar por até quantos adaptadores deve ser removidos de cada read.
# No caso do nosso comando, o parâmtro diz: remova até dois adaptadores de cada read.

# -m: Essa opção descarta reads trimadas que apresentem tamanho menor que X (valor que você opta).
# No comando abaixo, reads menores que 15 pb serão descartadas.

# Op-order: Esse parâmetro descreve a ordem que as operações de trimagem são aplicadas. As letras significam:
# -A: Trimagem do adaptador,C: Corte,G: Trimagem NextSeq,Q: Qualidade da trimagem,W: Sobrescrever leituras de baixa qualidade.
# O default é WCGQA entretanto no comando acima a ordem foi alterada para GAWCQ devido a experiências prévias.

# Match-read-wildcards: Modo de interpretação da IUPAC está habilitado, assim ele identifica caracteres coringas quando
# não se sabe qual base será inserida por exemplo.

# -O: Nesse parâmetro você fornece um valor minimo que indica que se a sobreposição entre a read e o adaptador for menor
# que esse valor, a read não será modificada.

# -q: Nesse parâmetro deve ser fornecido um valor de qualidade de bases, o qual abaixo desse as bases das extremidades
# serão removidas. Se quer que ambas as extremidades sejam trimadas, deve-se passar 2 valore, o primeiro para a extremidade
# 5’ e o segundo para a extremidade 3’. No comando acima apenas um valor foi passado: 25, assim apenas bases com score de
# qualidade abaixo de 25 da extreimidade 3’ serão trimadas.

# -T: Fornece o número de threads a ser utilizado para trimagem das reads.

# Correct-mismatches: Esse parâmetro indica como lidar com os mismatches durante o alinhamento.
# As opções são “Liberal, conservativa, e N”. O N significa substituir a base que deu mismatch por N.
# As opções de correção ‘Liberal e conservativa’ envolvem a susbtituição de bases por aquela que apresenta a maior
# qualidade. Elas diferem apenas quando os scores de qualidade são iguais. A correção liberal substitui de acordo com a
# melhor mediana das qualidades das bases, enquanto que a conservativa significa manter inalterada.

# Pair-filter: Nesse parâmetro deve-se escolher qual das reads (pair-end) deve corresponder ao critério de filtragem para
# que ocorra a filtragem. Você pode escolher que ambas as reads entrem no critério ou qualquer uma. No caso foi escolhido qualquer uma.

# -a: Nesse parametro deve-se fornecer a sequência do adaptador a ser removido na primeira read na região 3’.
# -A: Nesse parametro deve-se fornecer a sequência do adaptador a ser removido na segunda read na região 3’.

# -o: Local de saída da read 1 trimada.
# -p: Local de saída da read 2 trimada.
# -pe1: Primeiro arquivo de entrada
# -pe2: Segundo arquivo de entrada

# Untrimmed-output: Indica o local de saída das reads 1 que não foram trimadas.
# Untrimmed-paired-output: Indica o local de saída das reads 2 que não foram trimadas.

atropos trim --aligner insert \
-e 0.1 \
-n 2 \
-m 15 \
--op-order GAWCQ \
--match-read-wildcards \
-O 25 \
-q 28 \
-T 8 \
--correct-mismatches conservative \
--pair-filter any \
-a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-o ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \
-p ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \
-pe1 ${r1} \
-pe2 ${r2} \
--untrimmed-output ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq \
--untrimmed-paired-output ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq \
> ${outdir}/processed/atropos/${name}.atropos.log.out.txt \
2> ${outdir}/processed/atropos/${name}.atropos.log.err.txt

# As reads que não foram trimadas na etapa anterior são submetidas a esse outro algoritmo "adapter"

echo -e "Running atropos (adapter) for adapter trimming using sample ${name}: ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq & ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq ...\n"

# -pe1: É o arquivo de entrada da read 1 para a trimagem, entretanto são o arquivo de saída do comando anterior,
# que não foram trimadas pelo algoritmo ‘insert’.
# -pe2: É o arquivo de entrada da read 2 para a trimagem, entretanto são o arquivo de saída do comando anterior,
# que não foram trimadas pelo algoritmo ‘insert’.

# -g: Nesse parametro deve-se fornecer a sequência do adaptador a ser removido na primeira read na região 5’.
# -G: Nesse parametro deve-se fornecer a sequência do adaptador a ser removido na segunda read na região 5’.

atropos trim --aligner adapter \
-e 0.1 \
-n 2 \
-m 1 \
--match-read-wildcards \
-O 3 \
-q 28 \
--pair-filter both \
-pe1 ${outdir}/processed/atropos/${name}_R1.atropos_untrimmed.fastq \
-pe2 ${outdir}/processed/atropos/${name}_R2.atropos_untrimmed.fastq \
-a GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG \
-g AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
-G CAAGCAGAAGACGGCATACGAGATNNNNNNGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT \
-T 8 \
-o ${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
-p ${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}.atropos_adapter.log.out.txt \
2> ${outdir}/processed/atropos/${name}.atropos_adapter.log.err.txt

# Após esses comandos, ambas as saidas dos dois comandos são concatenadas em um só arquivo, utlizando o comando cat:
# -read1 trimada com algoritmo insert + read 1 trimada com algoritmo adapter.
# -read 2 trimada com algoritmo insert + read 2 trimada com algoritmo adapter.

echo -e "Merging atropos adapter trimming results using sample ${name}: ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq and ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq + ${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq and ${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq ...\n"

cat ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq

cat ${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq \
> ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq

# Em seguida, os arquivos gerados antes da concantenação foram removidos pois não serão utilizados.

echo -e "Removing useless atropos results ...\n"

rm -f ${outdir}/processed/atropos/${name}_R1.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R1.atropos_adapter.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_insert.fastq \
${outdir}/processed/atropos/${name}_R2.atropos_adapter.fastq

# Em seguida as reads são submetidas ao programa PrinSeq, que é uma ferramenta para garantir que os dados obtidos não estejam
# comprometidos por sequencias de baixa qualidade, artefatos ou contaminantes que podem levar a conclusões errôneas.
# Essa ferramenta, realiza tanto a trimagem, filtragem de sequencias de baixa qualidade, quanto gera um resumo estatistico
# das caracteristicas das sequencias, como um controle de qualidade, parecido com o fastqc.

echo -e "PrinSeq processing: ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq & ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq ...\n"

# Após as reads serem processadas pelo Atropos, servirão de entrada para o prinseq (paramentros:fastq e fastq2).
# -min_len: Filtrar as sequências menores que 25pb.
# -noniupac: Filtrar sequências que apresentem caracteres que não seja A,T,C,G ou N.

# -lc-method: Método a ser utilizado para filtragem de sequências de baixa complexidade, pode ser "dust" ou "entropy". No caso,
# foi selecionado o dust: Os scores são computados de acordo com a frequencia de diferentes trinucleotideos (0-100).
# Scores elevados implicam em baixa complexidade.

# -lc_threshold: Valor a ser usado para filtragem de sequencias de baixa complexidade, no caso, como escolhido o método
# "dust", é escolhido um valor máximo aceito. Acima desse score, as sequencias serão filtradas.

# -out_format: Utilizado para escolher o formato do output. No comando abaixo foi selecionado o formato fastq.

# -out_good: Por default a saida de dados é o diretorio em que se encontram as reads. Para alterar o diretório de
# saída e nome do arquivo de saída utiliza-se esse parametro. No comando abaixo o diretório prinseq como output e o nome
# dos arquivos de saída.

# -out_bad: Por default, a saída é o diretório input. Para alterar deve-se usar esse parametro. No caso foi utilizado o null
# para que esses arquivos bad não sejam gerados.

# -trim_qual_window: Tamnho da "janela" a ser usado para calcular o score de qualidade. No caso, está dizendo para parar na
# primeira base que falhar.

# -trim_qual_step: De quanto em quanto a janela desliza avaliando os scores de qualidade sem faltar nenhum.
# -trim_qual-right: Trima as sequências que apresentam score de qualidade menor que o threshold na extremidade 3'.
# -trim_qual_type: Tipo de calculo a ser usado para se obter o score de qualidade, foi selecionado a média.
# -trim_qual_rule: Regra a ser utilizada para comparação do score de qualidade gerado.Foi selecionado lt (less than).

prinseq-lite.pl -fastq ${outdir}/processed/atropos/${name}_R1.atropos_final.fastq \
-fastq2 ${outdir}/processed/atropos/${name}_R2.atropos_final.fastq \
-out_format 3 \
-trim_qual_window 1 \
-trim_qual_step 1 \
-trim_qual_right 30 \
-trim_qual_type mean \
-trim_qual_rule lt \
-out_good ${outdir}/processed/prinseq/${name}.atropos_final.prinseq \
-out_bad null \
-lc_method dust \
-lc_threshold 30 \
-min_len 25 \
-trim_tail_right 5 \
-trim_tail_left 5 \
-ns_max_p 80 \
-noniupac \
> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.out.log \
2> ${outdir}/processed/prinseq/${name}.atropos_final.prinseq.err.log

# Realização do controle de qualidade após ter sido realizado o pre-processamento dos dados

echo -e "FastQC pos-evaluation using sample ${name}: ${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1.fastq & ${outdir}/processed/prinseq/${name}_2.atropos_final.prinseq.fastq ...\n"

fastqc -t 2 \
${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1.fastq \
-o ${outdir}/processed/fastqc/pos/ \
> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_1.log.out.txt \
2> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_1.log.err.txt

fastqc -t 2 \
${outdir}/processed/prinseq/${name}_2.atropos_final.prinseq.fastq \
-o ${outdir}/processed/fastqc/pos/ \
> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_2.log.out.txt \
2> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_2.log.err.txt

# SE EXISTIR <SAMPLE_NAME>.atropos_final.prinseq_1_singletons.fastq
if [ -e "${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1_singletons.fastq" ]; then
fastqc -t 2 \
${outdir}/processed/prinseq/${name}.atropos_final.prinseq_1_singletons.fastq \
-o ${outdir}/processed/fastqc/pos/ \
> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_1_singletons.log.out.txt \
2> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_1_singletons.log.err.txt
fi

# SE EXISTIR <SAMPLE_NAME>.atropos_final.prinseq_2_singletons.fastq
if [ -e "${outdir}/processed/prinseq/${name}.atropos_final.prinseq_2_singletons.fastq" ]; then
fastqc -t 2 \
${outdir}/processed/prinseq/${name}.atropos_final.prinseq_2_singletons.fastq \
-o ${outdir}/processed/fastqc/pos/ \
> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_2_singletons.log.out.txt \
2> ${outdir}/processed/fastqc/pos/${name}.atropos_final.prinseq_2_singletons.log.err.txt
fi

done