.. _trinity:

=======
Trinity
=======

.. contents:: In this section:
   :local:
   :depth: 2

Description
===========

According to the `page of Trinity <https://github.com/trinityrnaseq/trinityrnaseq/wiki>`_, Trinity is a method for transcriptome reconstruction from RNA-seq data, developed at the Broad Institute and the Hebrew University of Jerusalem.

Available Versions
==================

* trinity/2.8.5 (default)
* trinity/2.15.1

Loading the Module
==================

.. code-block:: bash

   # Load Trinity
   module load trinity/2.15.1

   # Verify installation
   Trinity --version

Basic Assembly
==============

.. code-block:: bash
   :caption: submit_trinity_basic.sh

   #!/bin/bash
   #SBATCH -J trinity_basic
   #SBATCH -N 1
   #SBATCH -c 8
   #SBATCH -t 48:00:00
   #SBATCH --mem=64G

   export INPUT="reads.left.fq.gz reads.right.fq.gz"
   export OUTPUT="trinity_out_dir/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --left reads.left.fq.gz \
                     --right reads.right.fq.gz \
                     --SS_lib_type RF \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output trinity_out_dir

Strand-specific Assembly
========================

.. code-block:: bash
   :caption: submit_trinity_strand.sh

   #!/bin/bash
   #SBATCH -J trinity_strand
   #SBATCH -N 1
   #SBATCH -c 16
   #SBATCH -t 72:00:00
   #SBATCH --mem=128G

   export INPUT="reads_1.fq.gz reads_2.fq.gz"
   export OUTPUT="trinity_strand/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --left reads_1.fq.gz \
                     --right reads_2.fq.gz \
                     --SS_lib_type RF \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output trinity_strand \
                     --full_cleanup

Single-end Assembly
===================

.. code-block:: bash
   :caption: submit_trinity_single.sh

   #!/bin/bash
   #SBATCH -J trinity_single
   #SBATCH -N 1
   #SBATCH -c 8
   #SBATCH -t 48:00:00
   #SBATCH --mem=64G

   export INPUT="reads.fq.gz"
   export OUTPUT="trinity_single/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --single reads.fq.gz \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output trinity_single

Assembly with Multiple Samples
==============================

.. code-block:: bash
   :caption: submit_trinity_multiple.sh

   #!/bin/bash
   #SBATCH -J trinity_multiple
   #SBATCH -N 1
   #SBATCH -c 16
   #SBATCH -t 96:00:00
   #SBATCH --mem=128G

   export INPUT="samples.txt"
   export OUTPUT="trinity_merged/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --samples_file samples.txt \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output trinity_merged

.. code-block:: text
   :caption: samples.txt

   cond_A    sample1    sample1_R1.fq.gz    sample1_R2.fq.gz
   cond_A    sample2    sample2_R1.fq.gz    sample2_R2.fq.gz
   cond_B    sample3    sample3_R1.fq.gz    sample3_R2.fq.gz
   cond_B    sample4    sample4_R1.fq.gz    sample4_R2.fq.gz

Including Unpaired Reads
========================

.. code-block:: bash
   :caption: submit_trinity_with_unpaired.sh

   #!/bin/bash
   #SBATCH -J trinity_unpaired
   #SBATCH -N 1
   #SBATCH -c 16
   #SBATCH -t 72:00:00
   #SBATCH --mem=128G

   export INPUT="reads_1.fq.gz reads_2.fq.gz unpaired.fq.gz"
   export OUTPUT="trinity_with_unpaired/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --left reads_1.fq.gz \
                     --right reads_2.fq.gz \
                     --single unpaired.fq.gz \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output trinity_with_unpaired

Abundance Quantification
========================

.. code-block:: bash
   :caption: submit_abundance_estimation.sh

   #!/bin/bash
   #SBATCH -J abundance
   #SBATCH -N 1
   #SBATCH -c 8
   #SBATCH -t 24:00:00
   #SBATCH --mem=32G

   export INPUT="trinity_out_dir.Trinity.fasta reads_1.fq.gz reads_2.fq.gz"
   export OUTPUT="abundance/"

   module load trinity/2.15.1
   module load bowtie2/2.3.5.1
   module load samtools/1.15.1

   # Align reads with Bowtie2
   job-nanny bowtie2-build trinity_out_dir.Trinity.fasta trinity_index
   job-nanny bowtie2 -x trinity_index \
                     -1 reads_1.fq.gz \
                     -2 reads_2.fq.gz \
                     -p $SLURM_CPUS_PER_TASK \
                     | samtools sort -o aligned.bam

   # Estimate abundance
   job-nanny samtools index aligned.bam
   job-nanny samtools idxstats aligned.bam > counts.txt

Differential Expression Analysis
================================

.. code-block:: bash
   :caption: submit_expression_analysis.sh

   #!/bin/bash
   #SBATCH -J expression
   #SBATCH -N 1
   #SBATCH -c 8
   #SBATCH -t 24:00:00
   #SBATCH --mem=32G

   export INPUT="trinity_out_dir.Trinity.fasta"
   export OUTPUT="expression/"

   module load trinity/2.15.1
   module load R/4.0.2

   # Align all samples and quantify
   for sample in sample1 sample2 sample3 sample4; do
       /usr/bin/time -v bowtie2 -x trinity_index \
                     -1 ${sample}_1.fq.gz \
                     -2 ${sample}_2.fq.gz \
                     -p $SLURM_CPUS_PER_TASK \
                     | samtools sort -o ${sample}.bam
       
       samtools index ${sample}.bam
       samtools idxstats ${sample}.bam > ${sample}.counts
   done

   # Create count matrix
   cat > create_matrix.R << 'EOF'
   library(edgeR)
   library(limma)

   # Combine counts
   samples <- c("sample1", "sample2", "sample3", "sample4")
   counts <- data.frame(row.names=1:1000)  # adjust

   for (s in samples) {
       data <- read.table(paste0(s, ".counts"), 
                          col.names=c("ref", "len", "mapped", "unmapped"))
       counts[[s]] <- data$mapped
   }

   # Differential expression analysis
   group <- factor(c("A", "A", "B", "B"))
   dge <- DGEList(counts=counts, group=group)
   dge <- calcNormFactors(dge)
   design <- model.matrix(~group)
   v <- voom(dge, design)
   fit <- lmFit(v, design)
   fit <- eBayes(fit)

   # Results
   results <- topTable(fit, coef=2, n=Inf)
   write.csv(results, "diff_expression.csv")
   EOF

   Rscript create_matrix.R

Job Array for Multiple Samples
==============================

.. code-block:: bash
   :caption: submit_trinity_array.sh

   #!/bin/bash
   #SBATCH -J trinity_array
   #SBATCH --array=1-5
   #SBATCH -N 1
   #SBATCH -c 8
   #SBATCH -t 48:00:00
   #SBATCH --mem=64G

   SAMPLES=(
       "sample1"
       "sample2"
       "sample3"
       "sample4"
       "sample5"
   )

   SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
   export INPUT="${SAMPLE}_1.fq.gz ${SAMPLE}_2.fq.gz"
   export OUTPUT="${SAMPLE}_trinity/"

   module load trinity/2.15.1

   job-nanny Trinity --seqType fq \
                     --max_memory $(($SLURM_MEM_PER_NODE / 1024))G \
                     --left ${SAMPLE}_1.fq.gz \
                     --right ${SAMPLE}_2.fq.gz \
                     --CPU $SLURM_CPUS_PER_TASK \
                     --output ${SAMPLE}_trinity \
                     --full_cleanup

Results Analysis
================

.. code-block:: bash
   :caption: analyze_trinity.sh

   #!/bin/bash
   #SBATCH -J analyze_trinity
   #SBATCH -N 1
   #SBATCH -n 1
   #SBATCH -t 04:00:00
   #SBATCH --mem=8G

   module load trinity/2.15.1

   # Assembly statistics
   $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta > assembly_stats.txt

   # Extract long contigs
   awk '/^>/ {if (seqlen) print seqlen; print; seqlen=0; next} 
        {seqlen+=length($0)} END {print seqlen}' \
       trinity_out_dir/Trinity.fasta > contig_lengths.txt

   # Count contigs by size
   echo "Contigs > 1kb: $(awk '$1>1000' contig_lengths.txt | wc -l)" >> assembly_stats.txt
   echo "Contigs > 2kb: $(awk '$1>2000' contig_lengths.txt | wc -l)" >> assembly_stats.txt
   echo "Contigs > 5kb: $(awk '$1>5000' contig_lengths.txt | wc -l)" >> assembly_stats.txt

   # BUSCO (if available)
   if command -v busco &> /dev/null; then
       busco -i trinity_out_dir/Trinity.fasta \
             -l eukaryota_odb10 \
             -o busco_out \
             -m transcriptome
   fi

References
==========

- Documentation: https://github.com/trinityrnaseq/trinityrnaseq/wiki
- Tutorial: https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity
- Publication: https://www.nature.com/articles/nbt.1883

.. seealso::

   - :ref:`spades` - Genome assembly
   - :ref:`rsem` - Expression quantification
   - :ref:`salmon` - Transcriptome quantification
   - :ref:`running_simulations` - How to submit jobs