Sequence Mapping, Part. 2

In this article, I will walk you through how to conduct a de novo transcriptome analysis. If you remember from my previous article, de novo analysis is used when there is no genome reference to map the sequence reads against. Therefore you will need to create a transcriptome assembly from several sequences of read data usually from multiple samples or multiple tissues.

Before we start, we need to set up some tools for this exercise. Click on the hyperlinks below and get the latest version of the tools installed in your Linux/Ubuntu. You will need:
Trinity, Jellyfish, Salmon
Bowtie2, Samtools

Note: You should already have Bowtie2 and Samtools installed, if not install the tools and make sure you have added these tools to your PATH environment.

The example data we will be using is from this article here. Go ahead and download the data, 2 files, from the NCBI database. In that article, two data files are Pair-end reads therefore you will need to 'split' the files, resulting in 4 files of .fastq format. Do not compress the files using the --gzip command when processing.

The command should look something like this to split, if you need a refresher:
fastq-dump [insert SRR/SRA Accession ID] --split-3

Step 0. Trimming Raw Data

After you have downloaded and split the files, run the Trimmomatic to trim the reads. using the following command format:

java -jar [insert full pathway of Trimmomatic directory]/trimmomatic-0.39.jar PE -phred33 [input_forward.fq.gz] [input_reverse.fq.gz] [output_forward_paired.fq.gz] [output_forward_unpaired.fq.gz] [output_reverse_paired.fq.gz] [output_reverse_unpaired.fq.gz] ILLUMINACLIP:[pathway to Trimmomatic adapters directory]/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

Once this process has been completed you should have forward and reverse un/paired fastq files for each sample.

Step 1. Sequence data Concatenation

To create a comprehensive transcriptome assembly, you would need to concatenate the data. You will be using processed .fastq files when concatenating.
So for example if you have sample data Testis_1.fastq Testis_2.fastq, Liver_1.fastq, Liver_2.fastq, Skin_1.fastq Skin_2.fastq, etc. you will need to group the files like below to concatenate the data:
Testis_1.fastq Liver_1.fastq Skin_1.fastq Testis_2.fastq Liver_2.fastq Skin_2.fastq

Now using our example data you should have 4 fastq files that you have processed split and trimmed:
SRR15013283_1.fastq SRR15013283_2.fastq SRR15013284_1.fastq SRR15013284_2.fastq

In your terminal where your .fastq files are stored, use the cat command to concatenate the data. The command format is below:

cat [sample1_1.fastq][sample2_1.fastq][sample3_1.fastq] >> [name_your_file]_1.fastq

cat [sample1_2.fastq][sample2_2.fastq][sample3_2.fastq] >> [name_your_file]_2.fastq

To give you an example the command should look something like this:
cat SRR15013283_1.fastq SRR15013284_1.fastq >> merged_tissue_1.fastq

I named my concatenated file as 'merged' for this demonstration.
As a result, you should have 2 .fastq files that are concatenated from this demo.

Step 2. De Novo Assembly

We will be using the Trinity tool to start the transcriptome assembly. Use the following command format below:

$TRINITY_HOME/Trinity --seqType fq --left [concatenated_fastq_file_1] --right [concatenated_fastq_file_2] --output [name_your_output_directory] --max_memory 100G --CPU 8

# Note:
# --seqType: assigning read file format e.g. fq: fastq, fa : fasta)
# -- left, right: assigning forward, reverse reads
# --output: creating output directory but must have the word 'trinity' included in your output directory name
# --max_memory: the max memory when processing assembly
# --CPU: the number of CPU to use

Note: make sure your Trinity PATH is set as export TRINITY_HOME=/path/to/trinity/installation/dir as stated in the manual otherwise the tool will not run properly. Yours should look similar to mine below:
export TRINITY_HOME=/home/compio/de_novo/trinityrnaseq-v2.15.0

When running the command, it should look like this:
$TRINITY_HOME/Trinity --seqType fq --left merged_tissue_1.fastq --right merged_tissue_2.fastq --output trinity_out --max_memory 100G --CPU 8

Depending on the file size and number of files the processing time varies. This will take a while to finish so run this process on screen if you are doing this analysis on a remote server.

Once the process has finished there will be several files in the output directory. The Trinity.fasta file is the assembled transcriptome. This file will be used for later analysis.

If you want to conduct some statistics you will use the TrinityStats.pl which is located in the util directory inside your Trinity tool directory. To run the Trinity statistical analysis you will need to create a PATH environment. So in your path editor, you should have something like this:
export TRINITY_HOME=/home/compio/de_novo/trinityrnaseq-v2.15.0/util

To run the statistical analysis using the following command format:

$TRINITY_HOME/TrinityStats.pl [name of the output file].Trinity.fasta

Once finished, in your terminal it will show the result like this:

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':    85109
Total trinity transcripts:    141368
Percent GC: 37.97

########################################
Stats based on ALL transcript contigs:
########################################

    Contig N10: 4662
    Contig N20: 3383
    Contig N30: 2635
    Contig N40: 2111
    Contig N50: 1678

    Median contig length: 483
    Average contig: 930.85
    Total assembled bases: 131591818


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

    Contig N10: 4246
    Contig N20: 2945
    Contig N30: 2218
    Contig N40: 1674
    Contig N50: 1202

    Median contig length: 376
    Average contig: 713.62
    Total assembled bases: 60735113

Step 3. Gene Prediction

In this step, we will be using the assembled transcript to find protein-coding genes using the TransDecoder tool.

1. Extract transcripts that are >= 100 amino acids using the following command below:

    TransDecoder.LongOrfs -t [output Trinity file name].Trinity.fasta
   

   As a result, you will have longest_orfs.pep, longest_orfs.gff3, longest_orfs.cds, etc. Of those files, you will be using longest_orfs.pep for later work.

   Homology-based search

   We will be conducting a homology search and gene prediction using the NCBI BLAST program and SwissProt DB (protein database).

2. Download uniprot_sprot.fasta file from here (under the 'Download Links' section) and convert it into blast db format using the following format:

    makeblastdb -in uniprot_sprot.fasta -dbtype prot
   
    # -2.x.x: type in your tool version (if you are using old version)
    # -in: DB format you are trying to create
    # -dbtype: stating whether your fasta file is nucleotide or protein sequence. Nucleotide: nucl, Protein: prot
   

   As a result, you should have 8 files ending in: .fasta.pdb, .fasta.phr, .fasta.pin, .fasta.pjs, .fasta.pot, .fasta.psq, .fasta.ptf, .fasta.pto

   Once the blast db file has been created, locate your longest_orfs.pep file and use the following command to start the homology search:

    blastp -query [pathway to file]/longest_orfs.pep -db  uniprot_sprot.fasta -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 8  > blastp.outfmt6
   
    #### Note ####
    # -query:  sequence file
    # -db:  database
   
    # -max_target_seqs: number of query match you want to see
    # -outfmt:  output format information in the scale of 0~11. 6 is  tabular format.
   
    # -evalue:  blast evalue cutoff
    # -num_threads: number of cpu to use
   

   You should see a blastp.outfmt file format created once it has finished.

3. Gene prediction

   Using the blast result from the above step input the following command format:

    TransDecoder.Predict  -t  [Trinity.fasta file]   --retain_blastp_hits blastp.outfmt6  --cpu 8  --single_best_only
   
    ### Note: ###
    # -t: Trinity assembly file
    # --retain_blastp_hits: blastp output file from  previous step above
    # --cpu: number of cpu to use
    # --single_best_only: selecting the best orf of the generated transcripts
   

   The resulting files from this step are: .Trinity.fasta.transdecoder.pep, .Trinity.fasta.transdecoder.gff3, .Trinity.fasta.transdecoder.cds, .Trinity.fasta.transdecoder.bed

   Note: make sure the Trinity.fasta and blastp.outfmt6 files are in the same directory where trinity_out.Trinity.fasta.transdecoder_dir is located.

4. Removing Redundant Transcripts

   When conducting de novo assembly you have to consider the transcriptomic isoforms, therefore, you will need to remove the redundant transcripts using the CD-hit tool. Use the following command format to execute:

    cd-hit -i [Trinity.fasta.trasndecoder.pep file] -o NRCDS_Trinity.fasta.transdecoder.pep.cdhit -c 0.99 -T 8
   
    ### Note: ###
    # -i: gene prediction fasta file generated via TransDecoder
    # -o: output file name
    # -c: identity cutoff; 0.99 means of the 99% transcript cluster similarity, pick the longest transcript
    # -T: number of CPU to use
   

   Once the process has finished, you should see 2 files generated: Trinity.fasta.transdecoder.pep.cdhit, Trinity.fasta.transdecoder.pep.cdhit.clstr . The .cdhit file contains protein sequence information, .cdhit.clstr file contains cluster grouping information.

   As a result of these assembly files, we have generated Non-redundant protein Coding Sequences (NRCDS).

Step 4. Expression Level Quantification

In this step, you will be mapping the sequence reads onto the NRCDS and basing the mapped read counts on to estimate of the expression level.

1. NRCDS_Trinity.fasta Generation

   You would need the corresponding nucleotide sequence of NRCDS. Therefore you will use the sequence id from NRCDS file to match onto Trinity.fasta file to extract the nucleotide sequence using the following command:

    samtools faidx [Trinity.fasta.transdecoder.pep.cdhit file]
   

   As a result, you should have a .Trinity.fasta.transdecoder.pep.cdhit.fai file generated (NRCDS sequence file). From there we will extract only the sequence ID using this command:

    awk '{print$1}' filename.Trinity.fasta.fai | grep TRINITY > NCRDS_id.txt
    #and then
    sed 's/...$//' NCRDS_id.txt 
    #this gets rid of the '.pxx' ending since that signifies protein which we don't want
   

   When finished you should have an output file generated as NCRDS_id.txt which contains only the sequence IDs.

   We will now find the matching nucleotide using the sequence ID file and .Trinity.fasta file . To extract the sequence and its ID you will need to use samtools using the following command:

    xargs samtools faidx [Trinity.fasta file] <  NRCD_id.txt > NRCDS_sequence.fasta
   
    #you can name yours NRCDS_Trintiy if you want. I added '_sequence' to note that this is nucleotide sequence since there are so many files with 'Trinity' names.
   

   Once this process has been completed, your file content should look like this containing the nucleotide sequences :

    >TRINITY_DN36095_c0_g1_i1
    CCTTCTTCTTATCTAAAGCCACAACCCTAAACCCAACATTTCTGAGTGAAAATGAGCACC
    TGATCTGACGAGACAAACCGCTGTAGCACTCCTTGGTATAAATCGGTAGATCTGACTTGT
    AGGTCTCGTTGAGCAAAATTCATGAAAAGATCCGCTAACATTTGAAAAGACCTTATCAGA
    CGATGACTCGGGTCTTGTTTTGGGATCAGGATAGTTGAAGAAATCGTTCTGGAAGTTGAA
    GTCCCAAGTTTCGTCTGAAATTTCTTCGATTTCTTGAGACATAATGGACGACGTGATGCT [....]
    >TRINITY_DN36088_c0_g1_i1
    ATCATTTTTAATGTGAACAACATCCACGAGCGTGTCAGTTTCTAATTCGGAAAAAGTTGG
    GATGTTTGTTGTCTCTTCGAAATCAATTTTTATTAACTGCGAGTCTTCATCGTGGACTGC
    TGTGCAAAAAACTTGTTCATCATTAGTGAGATCTGCTGTGTCTGTACTTTCGCAGTTACA
    TTCGGAACTCTTGAAATCGGTCTGCGACGAGCGTGAAACAAGTTGGAGTCGATCGTTTAA
    CATTTGCAGACGCAGTGAAGATAATTCAATTTCGTATTGCTCTTCTTTGGATTTAAATGC
    TGAAATGGTCTCTTGCAGGCTG [....]
    >TRINITY_DN36057_c0_g1_i1
    AAAATATTTTGTGAAAACCAATCCATCTACCAGCAAAAAGAGATACTTGCCTTAAGAGAG
    TCTTTGAAGATTGCCAGAGCCGAAAATGAAAAACTTAAAAGATCTCTAGAGAACGAGATG
    AAAGAGAAAGAGAAACGGACGAGCGCCATGCAGCAGCAAATCATGTCTGCGAGACAGATG
    GAGGAAGCGAGGGAGAGCAAGATCAAGGAACTGCATGCTAAGGTGGAATGCAAAGAAGAA [....]
   

2. Read Alignment & Abundance Estimation

   In this step, we will be mapping the reads from the NRCDS_Trinity.fasta file. We will be using each sample file to concatenate for the reads using the following command:

    $TRINITY_HOME/align_and_estimate_abundance.pl   --transcripts  NRCDS_sequence.fasta  --seqType  fq  --left  [rawdata file 1].fastq   --right  [rawdata file 2].fastq   --est_method   RSEM  --aln_method  bowtie2  --prep_reference  --output_dir  [name output directory]  --thread_count  8
   
    ### Note: ###
    # the raw data files should be from the same sample
    # --transcripts: the assembly file you want to map (NRCDS_Trinity.fasta)
    # --seqType: read format (fq: fastq, fa:fasta)
    # --left, right: forward, reverse reads
    # --est_method: abundacne estimation method (RSEM)
    # --aln_method: read alignment method (bowtie2)
    # --prep_reference: assembly file indexing
    # --output_dir: name of output directory
    # --thread_count: number of CPU to use
   

   So to give you an example my command will look like this:

   $TRINITY_HOME/align_and_estimate_abundance.pl --transcripts NRCDS_sequence.fasta --seqType fq --left /pathway/to/SRR15013284_1_P.fastq --right /pathway/to/SRR15013284_2_P.fastq --est_method RSEM --aln_method bowtie2 --prep_reference --output_dir SRR15013284_rsem --thread_count 8

   As a result, you should have a .RSEM.isoforms.results file for each sample. For this demonstration, you should have 2 isoform result files. Each transcript files contain read count, TPM, FPKM information. Note that you are using RAW DATA PAIRED fastq files.

Step 5. Differentially Expressed Genes (DEGs)

You will use each sample to conduct an analysis. In this step, we will check the DEG for each sample.

1. Gene Expression Matrix

   From each RSEM.isoforms.results file of each sample, located in your output directory from the above step, we will merge them into one count matrix using the following command format:

    $TRINITY_HOME/abundance_estimates_to_matrix.pl  --est_method  RSEM --gene_trans_map none --out_prefix ABC [RSEM.isoforms.results file_1]  [RSEM.isoforms.results file_2] [RSEM.isoforms.results file_3]
   

   To give you an example my command will look like this:

   $TRINITY_HOME/abundance_estimates_to_matrix.pl --est_method RSEM --gene_trans_map none --out_prefix ABC RSEM.isoforms.results

   After this process, you should have .isoform.counts.matrix and .isoform.TPM.not_cross_norm files, etc.

2. Differential Expression Analysis

   In this step, you will be conducting DEG analysis for each sample using the Trinity DE analysis script. Following, below, is the command format:

    $TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl  --matrix  [.isoform.counts.matrix file] --method  edgeR  --output  edgeR_dir  --dispersion  0.1
   
    ### NOTE: ###
    # --matrix: refers to the count matrix file generated from previous step
    # --dispersion: if the samples are from same species we use the dispersion value of 0.1
   

   Once this process has finished the resulting files are in edgeR_dir.

3. TMM Normalization

   In this step, we will be using normalizing the gene expression level from each sample.

   First, we need to obtain the Transcript Length using the following command:

    cut -f 1,3,4 [.isoform.results file] > Trinity.trans_lengths.txt
    ### Note:###
    #-f 1,3,4: we just want column 1,3, and 4 from the isoform result file containing id, length, effective_length.
   

   To give you an example my command will look like this:
   cut -f 1,3,4 84.RSEM.isoforms.results > 84.Trinity.trans_lengths.txt

   Next, use the TMM normalization method (Trimmed mean of M-values) using the following Trinity command:

    $TRINITY_HOME/Analysis/DifferentialExpression/run_TMM_normalization_write_FPKM_matrix.pl --matrix [counts.matrix file] --length [Trinity.trans_lengths.txt file]
   

   Your command should look similar to this:
   $TRINITY_HOME/Analysis/DifferentialExpression/run_TMM_normalization_write_FPKM_matrix.pl --matrix ABC.isoform.counts.matrix --length 83.Trinity.trans_lengths.txt

   As a result, you should have .counts.matrix.TMM_normalized.FPKM file generated. You only need one normalized FPKM file for the rest of the workflow.

4. Identifying DEGs

   Move the .counts.matrix.TMM_normalized.FPKM file into the edgeR_dir .

   Use the DE Trinity script to conduct DEG analysis:

    $TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix [.counts.matrix.TMM_normalized.FPKM file] -C 2 -P 0.001
   
    ### Note: ###
    # -C: Log2 fold change value; 2 means DEGs cutoff when expression level is x4 high or low
    # -P: DEGs cutoff if corrected p-value (FDR) is less than 0.001
   

   To give you an example, your code should look similar to this:
   $TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix 83.isoform.counts.matrix.TMM_normalized.FPKM

Step 6. Annotation

1. Homology Search

   This homology search is different from the previous step.

    blastp -query [Trinity.fasta.transdecoder.pep.cd.hit file] -db path/way/to/uniprot_sprot.fasta -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 8  > annotation.blastp.outfmt6
   

   Note that Uniprot fasta file should be coming from your protein database directory. Therefore your command should look like this:
   blastp -query TrinityTrinity.fasta.transdecoder.pep.cdhit -db /var2/compbio/de_novo/Data/02.transDecoder/protein_DB/uniprot_sprot.fasta -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 8 > anno.blastp.outfmt6

2. Annotation Data Parsing

   Once anno.blastp.outfmt6 files have been generated, next is to extract the gene and protein ID from the file using the command below:

    awk '{print$1,$2}' annotation.blastp.outfmt6 | grep TRINITY > output.txt
    # followed by
    sed 's/\(^.*\)\.p1/\1/' output.txt > annotation.txt
   

   Once this process has been completed you can export the edgeR_dir to the local PC to conduct further analysis e.g. correlation, heatmap, PCA (if you have multiple samples). To transfer the directory use the following command format below:

    scp -r -P [port number] [username]@[IP address]:/path/to/edgeR_dir /path/to/local/PC