# all the tools used in this pipeline were download through conda

# Step 1: Raw sequencing quality check using fastqc and multiqc 

for filename in $(ls *'.fastq.gz')
do
fastqc $filename
done

multiqc .

# Step 2: quality filtering and trimming using bbduk 

## note: Be aware that Novaseq instrument reduces data footprints by binning Q-scores in only 4 bins (2,11,25,37).
## As a consequence, quality filtering at Q30 will remove a lot of high quality reads. 

export OUTDIR=path_to_outputDIR

for filename1 in $(ls *'_R1.fastq.gz')
do
fileID=`expr match "$filename1" '\(RA_[A-Za-z0-9._]*_L2\)'`
echo $fileID
filename2=$fileID'_R2.fastq.gz'
newname1=$fileID'_R1_no_adapter.fastq.gz'
newname2=$fileID'_R2_no_adapter.fastq.gz'

## adapter trimming
bbduk.sh in1=$filename1 in2=$filename2 out1=$OUTDIR/$newname1 out2=$OUTDIR/$newname2 ref=/media/bioinfoteam/SCRATCH_1TB/temp_metatranscriptomic/truseq_adapter.fasta ktrim=r k=23 mink=11 hdist=1 tpe tbo stats=$OUTDIR/${fileID}stats.txt
done

for newname1 in $(ls *'_R1_no_adapter.fastq.gz')
do
fileID=`expr match "$newname1" '\(RA_[A-Za-z0-9._]*_L2\)'`
echo $fileID
newname2=$fileID'_R2_no_adapter.fastq.gz'
trimmed1=$fileID'_R1_QC.fastq.gz'
trimmed2=$fileID'_R2_QC.fastq.gz'

## Quality trimming and filtering
bbduk.sh in1=$newname1 in2=$newname2 out1=$trimmed1 out2=$trimmed2 qtrim=rl trimq=20 minlen=50 maq=20 maxns=1 stats=${fileID}_stats.txt
done

### check quality after quality trimming and filtering

for filename in $(ls *'_QC.fastq.gz')
do
fastqc $filename
done

multiqc .

# Step 2: Sort rRNA using sortmeRNA 

for filename in $(ls *'.fastq.gz')
do
gzip -d $filename
done

export DATABASE=path_to_/rna_families

for filename1 in $(ls *'_R1_QC.fastq')
do
fileID=`expr match "$filename1" '\(RA_[A-Za-z0-9._]*_L2\)'`
filename2=$fileID'_R2_QC.fastq'
echo $fileID

export OUTFILE1=path_to_output/${fileID}_aligned
export OUTFILE2=path_to_output/${fileID}_non_aligned

sortmerna -ref $DATABASE/rfam-5.8s-database-id98.fasta -ref $DATABASE/rfam-5s-database-id98.fasta -ref $DATABASE/silva-arc-16s-id95.fasta -ref $DATABASE/silva-arc-23s-id98.fasta -ref $DATABASE/silva-bac-16s-id90.fasta -ref $DATABASE/silva-bac-23s-id98.fasta -ref $DATABASE/silva-euk-18s-id95.fasta -ref $DATABASE/silva-euk-28s-id98.fasta -reads $filename1 -reads $filename2 -fastx --aligned $OUTFILE1 --other $OUTFILE2 -paired_in True

rm /home/bioinfoteam/sortmerna/run/kvdb/*

done

## check quality after sortmeRNA

for filename in $(ls *'non_aligned.fq')
do
fastqc $filename
done

multiqc .

# de novo assembly using RNA_spades

spades.py --rna --pe1-1 RA_Ae_1_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_2_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_3_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_1_111022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_2_121022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_3_121022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_1_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_2_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_3_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_1_290323_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_2_290323_L2_R1_non_aligned.fastq --pe1-1 RA_Ae_3_290323_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_1_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_2_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_3_090123_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_1_111022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_2_121022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_3_121022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_1_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_2_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_3_241022_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_1_290323_L2_R1_non_aligned.fastq --pe1-1 RA_Fe_2_290323_L2_R1_non_aligned.fastq --pe1-2 RA_Ae_1_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_2_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_3_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_1_111022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_2_121022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_3_121022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_1_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_2_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_3_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_1_290323_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_2_290323_L2_R2_non_aligned.fastq --pe1-2 RA_Ae_3_290323_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_1_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_2_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_3_090123_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_1_111022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_2_121022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_3_121022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_1_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_2_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_3_241022_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_1_290323_L2_R2_non_aligned.fastq --pe1-2 RA_Fe_2_290323_L2_R2_non_aligned.fastq \
    -o RAgran_rnaspades_out

# mapping back the read onto the assembly using bowtie 2
## build index using the assembly
bowtie2-build -f path_to_RNAspadeoutput/hard_filtered_transcripts.fasta path_to_outputINDEX/deNovo_rnaSpades

## mapping
export DATA=path_to_rRNAsortedReads
export DB=path_to_outputINDEX/deNovo_rnaSpades
export OUT=path_to_outputBOWTIE2

cd $DATA
for filename1 in $(ls *'R1_non_aligned.fastq')
do
fileID=`expr match "$filename1" '\(RA_[A-Za-z0-9._]*_L2\)'`
echo $fileID
filename2=$fileID'_R2_non_aligned.fastq'
bowtie2 -p 10 -q --no-unal -x $DB -k 20 -1 $filename1 -2 $filename2 2>$OUT/${fileID}_align_stats.txt| samtools view -@10 -Sb -o $OUT/${fileID}_bowtie2.bam
done

## example in one sample
10310053 reads; of these:
  10310053 (100.00%) were paired; of these:
    771178 (7.48%) aligned concordantly 0 times
    3171512 (30.76%) aligned concordantly exactly 1 time
    6367363 (61.76%) aligned concordantly >1 times
    ----
    771178 pairs aligned concordantly 0 times; of these:
      67872 (8.80%) aligned discordantly 1 time
    ----
    703306 pairs aligned 0 times concordantly or discordantly; of these:
      1406612 mates make up the pairs; of these:
        903552 (64.24%) aligned 0 times
        128282 (9.12%) aligned exactly 1 time
        374778 (26.64%) aligned >1 times
95.62% overall alignment rate

# BUSCO

busco -i path_to_RNAspadeoutput/hard_filtered_transcripts.fasta -l bacteria_odb10 -o rnaSpades_BUSCO --out_path path_to_output_BUSCO -m transcriptome

--------------------------------------------------
	|Results from dataset bacteria_odb10              |
	--------------------------------------------------
	|C:65.3%[S:18.5%,D:46.8%],F:7.3%,M:27.4%,n:124    |
	|81	Complete BUSCOs (C)                       |
	|23	Complete and single-copy BUSCOs (S)       |
	|58	Complete and duplicated BUSCOs (D)        |
	|9	Fragmented BUSCOs (F)                     |
	|34	Missing BUSCOs (M)                        |
	|124	Total BUSCO groups searched               |

# Clustering using CD-HIT
cd-hit-est -i path_to_RNAspadeoutput/hard_filtered_transcripts.fasta -o path_to_RNAspadeoutput/hard_filtered_transcripts_cluster.fasta -c 0.95 -n 10 -d 0 -M 16000 -T 8

================================================================
                            Output                              
----------------------------------------------------------------
total seq: 1712386
longest and shortest : 110348 and 179
Total letters: 1326663130
Sequences have been sorted

Approximated minimal memory consumption:
Sequence        : 1603M
Buffer          : 8 X 61M = 492M
Table           : 2 X 44M = 88M
Miscellaneous   : 26M
Total           : 2211M

Table limit with the given memory limit:
Max number of representatives: 4000000
Max number of word counting entries: 1723528033

# mapping of teh reads after clustering
## build index 
bowtie2-build -f path_to_RNAspadeoutput/hard_filtered_transcripts_cluster.fasta path_to_index/deNovo_rnaSpades_cluster

export DATA=path_to_rRNAsortedReads
export DB=path_to_index/deNovo_rnaSpades_cluster
export OUT=path_to_output/bowtie2

cd $DATA
for filename1 in $(ls *'R1_non_aligned.fastq')
do
fileID=`expr match "$filename1" '\(RA_[A-Za-z0-9._]*_L2\)'`
echo $fileID
filename2=$fileID'_R2_non_aligned.fastq'
bowtie2 -p 10 -q --no-unal -x $DB -k 20 -1 $filename1 -2 $filename2 2>$OUT/${fileID}_align_stats.txt| samtools view -@10 -Sb -o $OUT/${fileID}_bowtie2.bam
done

# assign taxonomy usinh Kaiju
export OUTDIR=path_to_output_Kaiju
export INPUT=path_to_RNAspadeoutput
export DB=path_to_kaiju_DB/kaiju_db
kaiju -z 16 -a mem -t $DB/nodes.dmp -f $DB/kaiju_db_nr_euk.fmi -i $INPUT/hard_filtered_transcripts_cluster.fasta > $OUTDIR/RAgran_RNAspadesAssem_taxo.out

## add taxon info into output file
kaiju-addTaxonNames -t $DB/nodes.dmp -n $DB/names.dmp -r superkingdom,phylum,class,order,family,genus -i $OUTDIR/RAgran_RNAspadesAssem_taxo.out -o $OUTDIR/RAgran_RNAspadesAssem_taxo_names.out

# gene prediction and annotation using DRAM
DRAM.py annotate -i 'path_to_RNAspadeoutput/hard_filtered_transcripts_cluster.fasta' --min_contig_size 1500 -o path_to_output/DRAM

# eggNOG
path_to_eggNOGmapper/emapper.py -i path_to_output/DRAM/genes.faa -o RAgran_eggNOG --output_dir path_to_output/eggNOG

# feature count
export WORKING_DIR=path_to_output/bowtie2
export DATABASE=path_to_output/DRAM/genes.gff
export OUTPUT=path_to_output_FeatureCount/RAgran_rnaspades_countFC.txt
export CORRESP=path_to_correspondingFile

featureCounts -p -t CDS -O -M -C -B -A $CORRESP --fraction -T 8 -g ID -a $DATABASE -o $OUTPUT $WORKING_DIR/*.bam