##############################################################
# QUALITY FILTERING AND ADAPTER TRIMMING
##############################################################

for i in {1..6};
	do
		cd ./BM_$i # path to raw fastq files
		files=(*.gz)
		fastp -w 16 -q 30 -u 70 -l 30 -r -W 4 -M 20 -i ${files[0]} -I ${files[1]} -o Trimmed.${files[0]} -O Trimmed.${files[1]}
		cd ..
	done 


##############################################################
# PREPARE REFERENCES
##############################################################

# mRNA reference: WormBase C. elegans WS279 reference
rsem-prepare-reference -p 48 --star-sjdboverhang 150 --gtf c_elegans.PRJNA13758.WS279.canonical_geneset.gtf --star c_elegans.PRJNA13758.WS279.genomic.fa star_150

# transposon reference: RepBase C. elegans repetitive elements
STAR --runMode genomeGenerate --runThreadN 48 --genomeDir star_index_transposons/ --genomeFastaFiles transposable_elements_all.txt 


##############################################################
# ALIGNMENT AND mRNA FEATURE COUNTS
##############################################################

touch rsem_counts.txt
touch rsem_log

for i in {1..6};
	do
		cd ./BM_$i # path to trimmed fastq files
		files=(Trimmed*)
		rsem-calculate-expression -p 32 --star --paired-end --star-gzipped-read-file \
		--strandedness 'reverse' --estimate-rspd --append-names --output-genome-bam \
		${files[0]} ${files[1]} star_150 ws279_${files[0]} &>>../rsem_log
		printf "${files[0]}\t" >>../rsem_counts.txt
		awk '{sum += $5} END {print sum}' <${files[0]}.genes.results >>../rsem_counts.txt
		cd ..
	done 

##############################################################
# ALIGNMENT AND TRANSPOSON FEATURE COUNTS
##############################################################

for i in {1..6};
	do
		cd ./BM_$i # path to trimmed fastq files
		files=(Trimmed*)
		STAR --runThreadN 48 --genomeDir ~/star_index_transposons/  --readFilesCommand zcat --readFilesIn ${files[0]} --outFileNamePrefix transposons_${files[0]} --outSAMtype BAM SortedByCoordinate
		samtools index  transposons_${files[0]}Aligned.sortedByCoord.out.bam
		
# count reads for each transposon using samtools view
# merge counts with mRNA counts from RSEM
# Use DESeq2 for normalization and statistical analysis (see DESeq2 script)