---
# System prepended metadata

title: rnaseq.sh

---

# rnaseq.sh

#!/bin/bash

# Número de processadores a ser utilizado durante o processo. Criando a variável, basta substituir nos comandos.

num_threads=8

# É necessário passar um diretório input de onde vai estar contido os dados de entrada para o alinhamento, na linha de
# comando ao chamar o script para ser executado.

indir=$1

# SE(if) não (!) for passado um diretório de entrada, print na tela a mensagem "Missing input directory."
if [ ! ${indir} ]; then
	echo "Missing input directory."
	exit
fi
# A função echo, mostra na tela a mensagem que etá entre aspas.

# SE for passado um argumento na linha de comando, mais esse não for o diretório, print na tela a mensagem:"Wrong input directory (${indir})."
if [ ! -d ${indir} ]; then
	echo "Wrong input directory (${indir})."
	exit
fi

# É necessário passar um diretório de saída, ou seja, que vai conter os dados gerados.

outdir=$2

# SE  não for passado o argumento 2, o diretório de saída na linha de comando, print a mensagem:"Missing output directory."
if [ ! ${outdir} ]; then
	echo "Missing output directory."
	exit
fi

# SE  o argumento passado não é um diretório, print na tela a mensagem dizendo:"Wrong output directory (${outdir})."

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

# É necessário passar um terceiro argumento na linha de comando, o arquivo de anotação do genoma no formato gtf.

refgtf=$3

# SE ${refgtf} NÃO EXISTE, SE NÃO FOI PASSADO ARGUMENTO 3 NA LINHA DE COMANDO
if [ ! ${refgtf} ]; then
	echo "Missing GTF file."
	exit
fi

if [ ! -e "${refgtf}" ]; then
	echo "Not found GTF file (${refgtf})."
	exit
fi

# É necessário passar um quarto argumento na linha de comando: Arquivo FASTA referente ao genoma a ser utilizado como referência.

refseq=$4

# SE ${refseq} NÃO EXISTE, SE NÃO FOI PASSADO ARGUMENTO 4 NA LINHA DE COMANDO
if [ ! ${refseq} ]; then
	echo "Missing GENOME fasta file."
	exit
fi

if [ ! -e "${refseq}" ]; then
	echo "Not found GENOME fasta file (${refseq})."
	exit
fi

# Primeiramente é chamado para a execução o script preprocess.sh, que realizará o processamento dos dados e controle de
# qualidade. O mesmo necessita de 2 argumentos na linha de comando: O diretório de entrada de dados (Arquivos fastq) e
# um diretório para saída dos dados processados, que por sua vez serão o input do presente script.

./preprocess.sh "${indir}" "${outdir}"

# O comando mkdir cria diretórios. Assim os diretórios abaixo serão criados dentro do diretório passado como saída de dados
# na linha de comando ao chamar o script.

mkdir -p ${outdir}/star_index
mkdir -p ${outdir}/star_out_pe
mkdir -p ${outdir}/star_out_se
mkdir -p ${outdir}/star_out_final
mkdir -p ${outdir}/cufflinks
mkdir -p ${outdir}/cuffmerge
mkdir -p ${outdir}/stringtie
mkdir -p ${outdir}/stringmerge
mkdir -p ${outdir}/cuffcompare
mkdir -p ${outdir}/cuffquant

# A estrutura "for" serve como um loop, facilitando o script, pois assim não é necessário fazer para cada amostra
# separadamente.  A estrutura abaixo diz: para todos os dados contidos no diretório prinseq que finalizam com
# .atropos_final.prinseq_1.fastq, faça o alinhamento.
# "ls" é pra procurar dentro do diretório todos os arquivos que terminam com .atropos_final.prinseq_1.fastq.
# a função "sed" é utilizada para substituição.
# a função "touch" cria arquivos vazios.

for r1 in `ls ${outdir}/processed/prinseq/*.atropos_final.prinseq_1.fastq`; do
	r1_singletons=`echo ${r1} | sed 's/prinseq_1.fastq/prinseq_1_singletons.fastq/'`
	if [ ! -e "${r1_singletons}" ]; then
		touch ${r1_singletons}
	fi

	r2=`echo ${r1} | sed 's/prinseq_1.fastq/prinseq_2.fastq/'`

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

	r2_singletons=`echo ${r2} | sed 's/prinseq_2.fastq/prinseq_2_singletons.fastq/'`
	if [ ! -e "${r2_singletons}" ]; then
		touch ${r2_singletons}
	fi

	name=`basename ${r1} | sed 's/.atropos_final.prinseq_1.fastq//'`

# Se não exister o index do genoma, então faça o index.
# Para a construção do index a ferramenta STAR será utilizada.
# O index facilitará o mapeamento das reads contra o genoma de referência.

	if [ ! -e "${outdir}/star_index/SAindex" ]; then
		echo "Indexing genome (${refseq}) ..."


		STAR 	--runThreadN        ${num_threads} \
			--runMode           genomeGenerate \
			--genomeFastaFiles  ${refseq} \
			--genomeDir         ${outdir}/star_index \
			--sjdbGTFfile       ${refgtf} \
#			--genomeSAindexNbases 12 \ Apenas para genomas muito pequenos (escalonamento)
			--sjdbOverhang      149 \
		 > ${outdir}/star_index/STAR.index.log.out.txt \
		2> ${outdir}/star_index/STAR.index.log.err.txt

	fi

# --runThreadN: Número de processadores a ser utilizado, foi criado uma variável no inicio do script.
# --runMode: Direciona o STAR a construir o index do genoma.
# --genomeFastaFiles: Arquivo FASTA referente ao genoma de referência.
# --genomeDir: Diretório de saída que irá conter o index gerado.
# --sjdbGTFfile: Arquivo gtf - anotação do genoma de referência.
# --sjdbOverhang: Especifica o tamanho da sequencia genomica levando em consideração a anotação para ser usado na construção
# do banco de dados de junções de splicing. Esse tamanho geralmente é igual ao tamanho das reads -1. Logo, 149.
#  > log.out.txt: saída que contém as etapas da criação do index.
# 2> log.err.txt: segunda saída que contem erros que venham a ocorrer durante o processo.

	echo "STAR alignment PE with sample ${name}: ${r1} & ${r2} ..."

	mkdir -p ${outdir}/star_out_pe/${name}

# Para o alinhamento das reads contra o genoma de referência também será utlizado a ferramenta STAR.

	STAR	--runThreadN        ${num_threads} \
		--genomeDir         ${outdir}/star_index \
		--readFilesIn       ${r1} ${r2} \
		--outSAMstrandField intronMotif \
		--outFilterIntronMotifs RemoveNoncanonical \
		--sjdbGTFfile       ${refgtf} \
		--outFilterMultimapNmax 20 \
		--outFileNamePrefix ${outdir}/star_out_pe/${name}/ \
		--outSAMtype        BAM Unsorted \
		 > ${outdir}/star_out_pe/${name}/STAR.alignment_pe.log.out.txt \
		2> ${outdir}/star_out_pe/${name}/STAR.alignment_pe.log.err.txt

# --readFilesIn: Reads (no caso já processadas) a serem utilizadas no mapeamento.
# --outSAMstrandField: Para dados unstranded,as ferramentas Cufflinks/Cuffdiff requerem alinhamento de splicing com o
# o atributo de fita XS, que será gerado com esse parâmetro e a opção intronMotif.
# --outFilterIntronMotifs: Remove introns não canônicos.
# --outFilterMultimapNmax: Número máximo de alinhamentos multiplos permitido por read, se exceder esse valor, a read é
# considerada não mapeada; Foi utilizado o valor default.
# --outFileNamePrefix: prefixo a ser dado para o arquivo gerado.
# --outSAMtype: arquivo de saída do tipo BAM sem estar ordenado por nome ou por coordenada.

	echo "STAR alignment SE with sample ${name}: ${r1_singletons} & ${r2_singletons} ..."

	mkdir -p ${outdir}/star_out_se/${name}

# Para as reads que não foram encontradas seu par, denomindas singletons, realizar o alinhamento:

	STAR	--runThreadN        ${num_threads} \
		--genomeDir         ${outdir}/star_index \
		--readFilesIn       ${r1_singletons},${r2_singletons} \
		--sjdbGTFfile       ${refgtf} \
		--outSAMtype        BAM Unsorted \
		--outFilterMultimapNmax 20 \
		--outSAMstrandField intronMotif \
		--outFileNamePrefix ./$outdir/star_out_se/${name}/ \
		 > ./${outdir}/star_out_se/${name}/STAR.alignment_se.log.out.txt \
		2> ./${outdir}/star_out_se/${name}/STAR.alignment_se.log.err.txt

	echo "Merging STAR alignment PE & SE (${name}) ..."

# STAR_out é a junção do alinhamento das reads paired-end e das singletons

	mkdir -p ${outdir}/star_out_final/${name}

# Combinar resultados do alinhamento com reads paired-end e alinhamento com reads single-end (singletons):
# Utilizar a função samtools merge para junção de ambos os alinhamentos.
# -f: caso tenha o arquivo, sobreescreva.
# -n: ordena o alinhamento por nome.
# @: número de processadores.

	 samtools merge -@ ${num_threads} -f -n  ${outdir}/star_out_final/${name}/Aligned.out.bam \
                                                ${outdir}/star_out_pe/${name}/Aligned.out.bam \
                                                ${outdir}/star_out_se/${name}/Aligned.out.bam \
	 > ${outdir}/star_out_final/${name}/samtools.merge.log.out.txt \
	2> ${outdir}/star_out_final/${name}/samtools.merge.log.err.txt

	echo "Sorting STAR alignment final (${name}) ..."

# Ordenar o resultado do alinhamento por coordenadas genômicas, pois é uma exigência para executar o cufflinks posteriormente.
# Para a ordenação utilzar a função samtools sort:

	 samtools sort -@ ${num_threads} -o      ${outdir}/star_out_final/${name}/Aligned.out.sorted.bam \
                                                ${outdir}/star_out_final/${name}/Aligned.out.bam \
	 > ${outdir}/star_out_final/${name}/samtools.sort.log.out.txt \
	2> ${outdir}/star_out_final/${name}/samtools.sort.log.err.txt

	echo "Collecting alignment statistics (${name}) ..."

# Gerando as estatisticas do alinhamento utilizando o script SAM_nameSorted_to_uniq_count_stats.pl.
# É necessário passar  o arquivo de entrada (BAM final, contendo todos os alinhamentos) e o nome para o arquivo de saída.

	SAM_nameSorted_to_uniq_count_stats.pl ${outdir}/star_out_final/${name}/Aligned.out.bam > ${outdir}/star_out_final/${name}/Aligned.stats.txt

# Após o mapeamento das reads contra o genoma de referência, o transcriptoma será então montado utilizando o cufflinks:

	echo "Running Cufflinks (${name}) ..."

	mkdir -p ${outdir}/cufflinks/${name}

	cufflinks --output-dir ${outdir}/cufflinks/${name} \
		  --num-threads ${num_threads} \
		  --GTF-guide ${refgtf} \
		  --frag-bias-correct ${refseq} \
		  --multi-read-correct \
		  --library-type fr-unstranded \
		  --frag-len-mean 250 \
		  --frag-len-std-dev 40 \
		  --total-hits-norm \
		  --compatible-hits-norm \
		  --max-frag-multihits 20 \
		  --min-isoform-fraction 0.25 \
		  --max-intron-length 9000 \
		  --min-intron-length 80 \
		  --overhang-tolerance 10 \
		  --max-bundle-frags 999999 \
		  --max-multiread-fraction 0.65 \
		  --overlap-radius 20 \
		  --3-overhang-tolerance 300 \
		  ${outdir}/star_out_final/${name}/Aligned.out.sorted.bam \
		 > ${outdir}/star_out_final/${name}/cufflinks.log.out.txt \
		2> ${outdir}/star_out_final/${name}/cufflinks.log.err.txt

# --output-dir: Diretório para a saída de dados.
# --num-threads: Número de processadores a ser utilizado.
# --GTF-guide: Arquivo de anotação referente ao genoma de referência no formato gtf.
# --frag-bias-correct: Arquivo FASTA do genoma de referência.
# --multi-read-correct: Esse parametro solicita uma estimaçao inicial das reads mapeadas em multiplas posições.
# --library-type: Indica o tipo de biblioteca.
# --frag-len-mean: Média do tamanho dos fragmentos.
# --frag-len-std-dev: Desvio padrão referente ao tamnho dos fragmentos.
# --overhang-tolerance: Número de pb permitidos a entrar em um intron, caso seja mapeado corretamente.
# --total-hits-norm: Conta todos os fragmentos, incluindo aqueles não compativeis com nenhum transcrito de referência,
# levando ao número de mapeamentos, utilizando o FPKM.
# –compatible-hits-norm: Conta somente os fragmentos compativeis com algums transcrito de referência, levando ao número
# de mapeamento, utilizando o FPKM.
# --min-isoform-fraction: Após calcular a abundancia de isoformas para um gene, Cufflinks elimina transcritos que apresentem
# baixa abundancia, pois esses transcritos podem não ser montados de maneira confiável, podendo ser artefatos.Também elimina
# introns que apresentem poucos alinhamentos de splicing.
# --max-intron-length: Tamanho máximo dos introns.
# --min-intron-length: Tamanho minimo dos introns.
# --max-bundle-frags: Número maximo de fragmentos que um locus deve apresentar para que então passe para o proximo.
# --max-multiread-fraction: Porção de reads que podem ser mapeadas em multiplas posições.
# --overlap-radius: Transcrito que estão separados por menos que essa distancia são unidos e o gap é preenchido.
# --3-overhang-tolerance: Número de pares de bases que se excedem na porção 3', para ser permitido que se una ao transcrito ja montado.
# Será também realizada uma montagem utilizando o software Stringtie:

	echo "Running StringTie (${name}) ..."

	mkdir -p ${outdir}/stringtie/${name}

# Deve-se passar primeiramente o arquivo de alinhamento

	stringtie ${outdir}/star_out_final/${name}/Aligned.out.sorted.bam \
		-G ${refgtf} \
		-o ${outdir}/stringtie/${name}/transcripts.gtf \
		-p ${num_threads} \
		-f 0.25 \
		-a 10 \
		-j 5 \
		-c 6 \
		-g 20 \
		-M 0.65 \
		-A ${outdir}/stringtie/${name}/gene_abundance.txt

# -G: Arquivo de anotação referente ao genoma de referência no formato gtf.
# -o: Arquivo de saída - formato gtf.
# -p: Número de processadores a serem utilizados no processo.
# -f: Fração minima de abundância de isoformas.
# -a: Junções que não apresentam reads que alinharam através delas em pelo menos essa quantidade de bases em ambos os lados são filtradas. 
# -c: Minimo de cobertura permitido para os transcritos preditos.
# -g: Transcrito que estão separados por menos que essa distancia são unidos e o gap é preenchido.
# -M: Porção de reads que podem ser mapeadas em multiplas posições.
# -A: Saida do arquivo de abudancia dos genes - formato txt.

done

# Junção de todos os arquivos de saída "transcripts.gtf" do cufflinks em um só - assembly_GTF_list.txt, utilizando a 
# ferramenta cuffmerge.

echo "Running cuffmerge ..."

# O comando find encontra todos os arquivos com esse nome no diretório dado.

find ${outdir}/cufflinks/ -name 'transcripts.gtf' > ${outdir}/cuffmerge/assembly_GTF_list.txt

cuffmerge -o ${outdir}/cuffmerge/ \
	--ref-gtf ${refgtf} \
	--ref-sequence ${refseq} \
	--min-isoform-fraction 0.20 \
	--num-threads ${num_threads} \
	   ${outdir}/cuffmerge/assembly_GTF_list.txt \
	 > ${outdir}/cuffmerge/cuffmerge.log.out.txt \
	2> ${outdir}/cuffmerge/cuffmerge.log.err.txt

# -o: Saída do arquivo gerado.
# -ref-gtf: Arquivo de anotação referente ao genoma no formato gtf.
# - ref-sequence: Arquivo no formato FASTA do genoma de referência.
# -min-isoform-fraction: minimo de abundancia de isoformas, abaixo disso são eliminadas.
# -num-threads: Número de processadores a serem utilizados.


# Similarmente ao processo anterior, se realizará a junção dos arquivos gerados do StringTie em um só arquivo, utilizando a
# ferramenta stringtie merge:

echo "Running stringtie merge ..."

find ${outdir}/stringtie/ -name 'transcripts.gtf' > ${outdir}/stringmerge/assembly_GTF_list.txt

stringtie --merge \
	-G ${refgtf} \
	-o ${outdir}/stringmerge/merged.gtf \
	-c 1 \
	-T 1 \
	-f 0.25 \
	-g 20 \
	-i \
	${outdir}/stringmerge/assembly_GTF_list.txt

# -T: minimo de transcritos TPM a serem utilizados na junção.
# -i: mantém os transcritos unidos com retenção dos introns.


# Em seguida será realizada uma comparação das montagens obtida com a referência utilizando a ferramenta cuffcompare:

cuffcompare	-r ${refgtf} \
		-s ${refseq} \
		-o ${outdir}/cuffcompare/stringmerge \
		${outdir}/stringmerge/merged.gtf \
		 > ${outdir}/stringmerge/cuffcompare.log.out.txt \
		2> ${outdir}/stringmerge/cuffcompare.log.err.txt

# Criação da variávek biogroup_label:

biogroup_label=()

for bamfile in `ls ${outdir}/star_out_final/*/Aligned.out.sorted.bam`; do
	name=`basename $(dirname ${bamfile})`

# Quantificação dos transcritos utilizando a ferramenta cuffquant:

	echo "Running cuffquant using sample ${name} with ${outdir}/stringmerge/merged.gtf as reference ..."

	mkdir -p ${outdir}/cuffquant/${name}

	cuffquant 	--output-dir ${outdir}/cuffquant/${name} \
			--frag-bias-correct ${refseq} \
			--multi-read-correct \
			--num-threads ${num_threads} \
			--library-type fr-unstranded \
			--frag-len-mean 250 \
			--frag-len-std-dev 40 \
			--max-bundle-frags 9999999 \
			--max-frag-multihits 20 \
			${outdir}/stringmerge/merged.gtf \
			${bamfile} \
		 > ${outdir}/cuffquant/${name}/cuffquant.log.out.txt \
		2> ${outdir}/cuffquant/${name}/cuffquant.log.err.txt

# Os parâmetros são similares ao cufflinks.

	groupname=`echo ${name} | sed 's/[0-9]\+$//'`
	biogroup_label=($(printf "%s\n" ${biogroup_label[@]} ${groupname} | sort -u ))

done

biogroup_files=()

# Para realizar a análise de genes diferencialmente expressos,utilizaremos os arquivo de abundância gerados na etapa anterior.

echo "Running Differential Expression Analysis ..."

# Será utilizado novamente a estrutura for, para fazer com todas as amostras de uma só vez.
# Preparo dos dados para cuffdiff

for label in ${biogroup_label[@]}; do
	echo -e "\tCollecting .cxb files for ${label} ..."
	group=()
	for cxbfile in `ls ${outdir}/cuffquant/${label}*/abundances.cxb`; do
		echo -e "\t\tFound ${cxbfile}"
		group=(${group[@]} "${cxbfile}")
	done
	biogroup_files=(${biogroup_files[@]} $(IFS=, ; echo "${group[*]}") )
done

echo -e "\tRunning cuffnorm & cuffdiff ..."
echo -e "\t\tLabels.: " $(IFS=, ; echo "${biogroup_label[*]}")
echo -e "\t\tFiles..: " ${biogroup_files[*]}
echo -e "\t\t\tGenerating abundance matrices (cuffnorm) ..."

# Normalização dos resultados utilizando o cuffnorm:

mkdir -p ${outdir}/cuffnorm/

cuffnorm 	--output-dir ${outdir}/cuffnorm \
 		--labels $(IFS=, ; echo "${biogroup_label[*]}") \
 		--num-threads ${num_threads} \
		--library-type fr-unstranded \
 		--library-norm-method geometric \
		--output-format simple-table \
 		${outdir}/stringmerge/merged.gtf \
 		${biogroup_files[*]} \
 	 	> ${outdir}/cuffnorm/cuffdiff.log.out.txt \
 		2> ${outdir}/cuffnorm/cuffdiff.log.err.txt

# -output-format: formato da saída: Tabela simples
# Restante dos parametros já foram citados no cufflinks.

echo -e "\t\t\tAnalysing differential expression (cuffdiff) ..."

# Análise de genes diferencialmente expressos utilizando o cuffdiff

mkdir -p ${outdir}/cuffdiff/ 

cuffdiff 	--output-dir ${outdir}/cuffdiff \
 		--labels $(IFS=, ; echo "${biogroup_label[*]}") \
 		--frag-bias-correct ${refseq} \
 		--multi-read-correct \
 		--num-threads ${num_threads} \
 		--library-type fr-unstranded \
 		--frag-len-mean 250 \
 		--frag-len-std-dev 40 \
 		--max-bundle-frags 9999999 \
 		--max-frag-multihits 20 \
 		--total-hits-norm \
 		--min-reps-for-js-test 2 \
 		--library-norm-method geometric \
 		--dispersion-method per-condition \
 		--min-alignment-count 12 \
 		${outdir}/stringmerge/merged.gtf \
 		${biogroup_files[*]} \
 	 	> ${outdir}/cuffdiff/cuffdiff.log.out.txt \
 		2> ${outdir}/cuffdiff/cuffdiff.log.err.txt

# --min-reps-for-js-test: Não será testada a expressão genica diferencial se não houver esse número minimo de réplicas.
# --library-norm-method: Normalização do método. Geométrica: FPKMs e contagem dos fragmentos são escalados com base na
# mediana das médias geométricas da contagem de fragmentos em todas as bibliotecas.
# --dispersion-method: Método de dispersão. Por condição: Cada réplica recebe seu próprio modelo, não é realizado um pool.
# --min-alignment-count:Número minimo de alinhamentos em um loco para esse ser submetido ao teste.
# Restante dos parametros já foram citados anteriormente