Welcome to the blog

Posts

My thoughts and ideas

Introduction to Transcript Annotation | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Introduction to Transcript Annotation

Course

Module 7 - Key concepts

  • Introduction to transcript annotation

Module 7 - Learning objectives

  • Explore methods to determine biological function from transcript sequences.
  • Differentiate between homology-based and sequence composition-based functional inference.

Module 7 - Lecture

  • (RNASeq_Module7_Lecture.pdf)(https://github.com/griffithlab/rnabio.org/blob/master/assets/lectures/cbw/2019/full/RNASeq_Module7_Lecture.pdf)
Course
Trinity | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Trinity

Course

Intro to Genome-guided RNA-Seq Assembly

To make use of a genome sequence as a reference for reconstructing transcripts, we’ll use the Tuxedo2 suite of tools, including Hisat2 for genome-read mappings and StringTie for transcript isoform reconstruction based on the read alignments. We’ll also explore using Trinity in genome-guided mode, performing a de novo assembly for reads aligned and clustered along the reference genome.

The data set that we’ll be leveraging for this application are human RNA-Seq data corresponding to a subset of genes on the human X chromosome, all based on data provided as part of the Tuxedo2 protocol paper.

This exercise provides a quick introduction to the Tuxedo2 toolkit and leverages only a single small pair of fastq files and a ‘miniaturized’ version of the X chromosome containing ~100 genes. For more detailed exploration of genome-guided assembly, you are encouraged to read the Tuxedo2 protocol paper and leverage the complete data sets involving multiple sample replicates.

Data Overview

The data sets we’ll use for genome-guided assembly are located at:

~/CourseData/RNA_data/trinity_trinotate_tutorial/mini_humanX/

If we list this path, we can see that we have a ‘minigenome.fa’ genome fasta file and corresponding ‘minigenome.gtf’ file providing the corresponding gene structure annotations. A single pair of fastq files is provided as ‘reads_1.fq.gz’ and ‘reads_2.fq.gz’.

 ls ~/CourseData/RNA_data/trinity_trinotate_tutorial/mini_humanX/
 minigenome.fa  minigenome.gtf  reads_1.fq.gz  reads_2.fq.gz

How many reads are there?

Setting up your environment

For convenience, we’ll be making use of certain environmental variables, such as $TRINITY_HOME to define where the Trinity software is installed.

To configure your environment, simply run the following command:

source ~/CourseData/RNA_data/trinity_trinotate_tutorial/environment.txt

Now, if you type

env | grep TRINITY

you should see:

TRINITY_HOME=/usr/local/trinityrnaseq-Trinity-v2.8.4/

Setting up your workspace

In your workspace directory, create a new workspace called ‘workspace_GG’ (with the ‘GG’ short for genome-guided).

cd ~/workspace  # change to home directory
mkdir workspace_GG
cd workspace_GG

Create symbolic links (shortcuts) to our reference genome, annotation, and read files:

ln -s ~/CourseData/RNA_data/trinity_trinotate_tutorial/mini_humanX/* .

At this point, you’re familiar with FASTA (.fa) files and FASTQ (.fq) files. The .gtf file provides the genome annotation and it’s formatted like so:

head minigenome.gtf

    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_172246";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform b precursor"; orig_coord_info "chrX,NM_172246,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_172245";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform a precursor"; orig_coord_info "chrX,NM_172245,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_006140";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform a precursor"; orig_coord_info "chrX,NM_006140,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_001161529";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform a precursor"; orig_coord_info "chrX,NM_001161529,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_001161532";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform h"; orig_coord_info "chrX,NM_001161532,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_001161531";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform g precursor"; orig_coord_info "chrX,NM_001161531,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_172249";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform e precursor"; orig_coord_info "chrX,NM_172249,1268800,1268879,+";
    minigenome	hg38_refGene	exon	1001	1080	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NR_027760";  gene_name "CSF2RA"; orig_coord_info "chrX,NR_027760,1268800,1268879,+";
    minigenome	hg38_refGene	exon	3651	3738	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_172246";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform b precursor"; orig_coord_info "chrX,NM_172246,1274755,1274842,+";
    minigenome	hg38_refGene	exon	3651	3714	0.000000	+	.	gene_id "CSF2RA"; transcript_id "CSF2RA^NM_172245";  gene_name "CSF2RA"; product "granulocyte-macrophage colony-stimulating factor receptor subunit alpha isoform a precursor"; orig_coord_info "chrX,NM_172245,1274755,1274818,+";

GTF format is a popular way to store genome feature annotations.

Tuxedo2

The Tuxedo2 protocol involves first aligning reads to the genome using hisat2, followed by transcript reconstruction using StringTie.

Aligning reads to the genome using Hisat2

Before we can align reads to the genome, we must index for use with hisat2. This involves a few operations:

Extract splice sites from intron-containing transcript structures:

hisat2_extract_splice_sites.py minigenome.gtf > minigenome.gtf.ss

Extract the exon records:

hisat2_extract_exons.py  minigenome.gtf > minigenome.gtf.exons

Now build the hisat2 index of the genome, leveraging the splice sites and exon data extracted above:

hisat2-build --exon minigenome.gtf.exons --ss minigenome.gtf.ss minigenome.fa minigenome.fa

After the genome index has been created, we’re ready to align reads to the genome. The native output from hisat2 is SAM and so we pipe this into samtools to convert to BAM and coordinate-sort the output.

hisat2 --dta -x minigenome.fa --max-intronlen 5000 \
        -1 reads_1.fq.gz -2 reads_2.fq.gz \
        | samtools view -Sb \
        | samtools sort -o alignments.hisat2.bam

The above command limits the intron length to 5k only because our ‘minigenome’ was constructed to limit intron lengths to a max of 5kb. The minigenome is far more dense than the true human X chromosome.

Examine the contents of the bam file using ‘samtools’:

samtools view alignments.hisat2.bam | head

    ERR188245.78785_HWI-ST661:130:C037KACXX:6:1101:8916:31520	137	minigenome	2183	60	4S71M	=	2183	0	GTACGCCTGTAATCCCAGTGACTTGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGT	C@@FFFFFHHHHHJJJJJGHHIJJJJIIJJJJIIJJIIIJEGHJIIIIIJFGGIJJEHFABDBDD;BB2?CBCC@	AS:i:-14	XN:i:0	XM:i:2	XO:i:0	XG:i:0	NM:i:2	MD:Z:14C0T55	YT:Z:UP	NH:i:1
    ERR188245.78785_HWI-ST661:130:C037KACXX:6:1101:8916:31520	69	minigenome	2183	0	*	=	2183	0	GGATAATCATGTTCTTGATACACATTTCCTTTTTTTTTGAGATGGAGTCTCGCTGTTGCCCAGGCTGGAGTACAG	CCCFFFFFHHHHHJJJJIJJJJJJJJJJJJJJJJJJJJBC@FCHCCGC==E9@EAEH:@DFFDDDEDDABCDDDE	YT:Z:UP
    ERR188245.525513_HWI-ST661:130:C037KACXX:6:1101:18753:196335	73	minigenome	2183	60	1S74M	=	2183	0	CGCCTGTAATCCCAGTGACTTGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAG	CCCFFFFFHHHHHJJHIJJJJJJJIJJIJIIJJIJJJGGIJGGGGGICFHIJJIIHEDDBDDBBD<AC?CD:AAC	AS:i:-12	XN:i:0	XM:i:2	XO:i:0	XG:i:0	NM:i:2	MD:Z:14C0T58	YT:Z:UP	NH:i:1
    ERR188245.525513_HWI-ST661:130:C037KACXX:6:1101:18753:196335	133	minigenome	2183	0	*	=	2183	0	CGCCAACTCTTCTATATAATCAGTTTGATGATCTGAATTAGAAAATACCCCTGGATAATCATGTTCTTGATACAC	C@CFFFFFHHHHHJJJJJJJJIJHJJJJJJJJJJGIIGIIJJJJJCIJJJHIJJJJIJJJJJJJJJJIJJJJJIJ	YT:Z:UP
    ERR188401.42420_HWI-ST0814:100:D0BT5ACXX:4:1101:20732:26752	69	minigenome	2183	0	*	=	2183	0	ATCAAGCCAGCCTGACCAACATGGTGAAACCCCATCTGTACTAAACATAAAAAAATTAGCCTGGCATGGTGGTGT	@CCFFFFFHHHHHJJJIJJIJJJJHGIIJIIJJJIJJJHFHIIJICHIJIJJIIIJGICEHGHFFFFEE;CD@?B	YT:Z:UP
    ERR188401.66672_HWI-ST0814:100:D0BT5ACXX:4:1101:3195:41421	73	minigenome	2183	60	75M	=	2183	0	GCCTGTAATCCCAGTGACTTGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGC	CCCFFFFFHHHHHJFHJJJJJJJIJJJJJIJJJJIJIIIIGAFDHIGIIIIJJJHEBD?BDBBD8?CD@CACACD	AS:i:-10	XN:i:0	XM:i:2	XO:i:0	XG:i:0	NM:i:2	MD:Z:14C0T59	YT:Z:UP	NH:i:1
    ERR188401.66672_HWI-ST0814:100:D0BT5ACXX:4:1101:3195:41421	133	minigenome	2183	0	*	=	2183	0	TTTTGAGATGGAGTCTCGCTGTTGCCCAGGCTGGAGTACAGTGGCGCGATCTTAGCTCACTGCAACCTCCGCCTC	CCCFFFFFHHHHHJJJIJJIJHJJJJJJJJJIEBHIBFHIIECHIJJGIGHIIJJHHHEHFFFFFFFEEDDDDDD	YT:Z:UP
    ERR188401.42420_HWI-ST0814:100:D0BT5ACXX:4:1101:20732:26752	153	minigenome	2183	60	1S74M	=	2183	0	CGCCTGTAATCCCAGTGACTTGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAG	DDEEDCEEFFFFHHGHECFJIGIGFJHIIHGIGJJJIHGJJIJJJJJJHGJIJJJJJJJJJJHGHHHFFFFFCCC	AS:i:-11	XN:i:0	XM:i:2	XO:i:0	XG:i:0	NM:i:2	MD:Z:14C0T58	YT:Z:UP	NH:i:1
    ERR188428.130590_HWI-ST0814:100:D0BT5ACXX:4:1101:14424:73163	69	minigenome	2185	0	*	=	2185	0	AAAAAATTAATGTGTATCCAGGACATTTTAAAAACCTGTACACAGTGTTTATTGTGGTTAGGAAGCAATTTCCCA	CCCFFFFFHHHHHJIIJJJJJIJJJJJJJJJJJIHIJJJIJJJIJGGGIJJJJJIJJJJJJJJJJJJJJIHGHHG	YT:Z:UP
    ERR188428.130590_HWI-ST0814:100:D0BT5ACXX:4:1101:14424:73163	153	minigenome	2185	60	73M2S	=	2185	0	CTGTAATCCCAGTGACTTGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCGG	DDDDDEDDEFFFFFFFHHHHJJIJIJIIJJJJJJJJJJJJJJJJHF@IJJJJJJJJJJJJJJHHGHHFFFFFCCC	AS:i:-12	XN:i:0	XM:i:2	XO:i:0	XG:i:0	NM:i:2	MD:Z:12C0T59	YT:Z:UP	NH:i:1

What do the different fields represent? Which is the sam flag?

Use the https://broadinstitute.github.io/picard/explain-flags.html utility to decipher the different flag values

Decoding Sam Flags

Examine the reads aligned to the genome using IGV. Don’t forget to first index the .bam (to .bai) and .fa (to .fai) files before loading into IGV. (Do you remember the indexing commands for .bam and .fa files? hint: use ‘samtools’). Also load in the reference annotation file ‘minigenome.gtf’ to provide additional perspective.

Transcript reconstruction using StringTie

StringTie reconstructs transcripts from the aligned reads, leveraging the .bam file as input, and generating a GTF file containing transcript structures as output.

stringtie -o stringtie.gtf alignments.hisat2.bam

Examine the stringtie.gtf output:

head stringtie.gtf

# stringtie -o stringtie.gtf alignments.hisat2.bam
# StringTie version 1.3.5
minigenome	StringTie	transcript	36732	41779	1000	-	.	gene_id "STRG.1"; transcript_id "STRG.1.1"; cov "3.432668"; FPKM "1150.780029"; TPM "2182.717041";
minigenome	StringTie	exon	36732	37082	1000	-	.	gene_id "STRG.1"; transcript_id "STRG.1.1"; exon_number "1"; cov "3.341223";
minigenome	StringTie	exon	41667	41779	1000	-	.	gene_id "STRG.1"; transcript_id "STRG.1.1"; exon_number "2"; cov "3.716716";
minigenome	StringTie	transcript	47357	81576	1000	-	.	gene_id "STRG.2"; transcript_id "STRG.2.1"; cov "9.523120"; FPKM "3192.565186"; TPM "6055.428711";
minigenome	StringTie	exon	47357	47707	1000	-	.	gene_id "STRG.2"; transcript_id "STRG.2.1"; exon_number "1"; cov "4.632516";
minigenome	StringTie	exon	52409	52531	1000	-	.	gene_id "STRG.2"; transcript_id "STRG.2.1"; exon_number "2"; cov "11.336511";
minigenome	StringTie	exon	53825	53968	1000	-	.	gene_id "STRG.2"; transcript_id "STRG.2.1"; exon_number "3"; cov "14.663729";
minigenome	StringTie	exon	54834	54966	1000	-	.	gene_id "STRG.2"; transcript_id "STRG.2.1"; exon_number "4"; cov "15.149121";

Notice that in addition to reconstructing transcripts, StringTie also provides expression values in TPM and FPKM values along with read coverage stats for the transcript and individual exons.

Load the stringtie.gtf file into IGV. How do the stringtie transcript structures compare to the reference transcripts?

Genome-guided Trinity de novo assembly

If your RNA-Seq sample differs sufficiently from your reference genome and you’d like to capture variations within your assembled transcripts, you might consider performing a genome-guided de novo assembly. Trinity can accept a bam file containing genome-aligned rna-seq reads as input. Reads are partitioned into coverage groups along the reference genome and each read cluster is assembled using the standard Trinity de novo assembly. Here, de novo assembly is restricted to only those reads that map to the genome. The advantage is that reads that share sequence in common but map to distinct parts of the genome will be targeted separately for assembly. The disadvantage is that reads that do not map to the genome will not be incorporated into the assembly. Unmapped reads can, however, be targeted for a separate genome-free de novo assembly.

Trinity assembly from genome-aligned reads (bam file)

Run genome-guided Trinity leveraging our hisat2-aligned reads like so:

$TRINITY_HOME/Trinity --genome_guided_bam alignments.hisat2.bam \
           --CPU 2 --max_memory 1G --genome_guided_max_intron 5000 \
	   --output trinity_out_dir_GG

Once Trinity completes, you’ll once again a trinity_out_dir_GG/ in your new workspace, and in this case it’ll contain the resulting assembly as ‘trinity-GG.fasta’.

Examine this ‘trinity-GG.fasta’ file:

less trinity_out_dir_GG/Trinity-GG.fasta

    >TRINITY_GG_1_c0_g1_i1 len=495 path=[1:0-494] [-1, 1, -2]
    AGTTATTCAAGTTGTAAAAGGTTATACAATAATTTAACAACTACCTTTTTTATTCTGTCG
    GGTTACTGACCTCACTTTATGTAAATACTTCGCATGACAAATTCAGTAACTCGTCTATTT
    CAGCATGCATAAGACTTTTCACTAGGGAAACTGATAAAGCTTGAGTCAACTAAATCTGCC
    TTCATACTTTATCAAGGGGAACCAAGCCTGCTGTGCTTACATCAGCATCTGGAAGACTTT
    CCTCTCCTCTAATCTGTGTACACATCTCCAAGCAAGGAAGAAAAAACAAACTCTGCTCAG
    ACGCCTATGAAACACCTGAATGAACTTTGATGAAGTACAGTCTGAGTTACCATCATGCAC
    AAGTAGAACTGCTCTTGGACTTGTTTTCCTGTTGTTTGTGGAACCTACGCGTTTGAATGG
    CTTGAACGTTGCATCTTTTAAAGTTATTTTTTAAGGGTTCTTGGCATTTATCCTAGTTGT
    CCGTGTTTGGCAATG
    >TRINITY_GG_2_c0_g1_i1 len=227 path=[1:0-226] [-1, 1, -2]
    TAGAGGAGAAAATTTCTATGGTCTAGATATTACTTGTAAAGACATGAGAAACCTGAGGTT
    CGCTTTGAAACAGGAAGGCCACAGCAGAAGAGATATGTTTGAGATCCTCACGAGATACGC
    GTTTCCCCTGGCTCACAGTCTGCCATTATTTGCATTTTTAAATGAAGAAAAGTTTAACGT
    GGATGGATGGACAGTTTACAATCCAGTGGAAGAATACAGGATGCCGG

Aligning Trinity-assembled Transcripts to the Genome

We’ll use the GMAP software to align the Trinity transcripts to our reference genome. Trinity contains a utility that facilitates running GMAP, which first builds an index for the target genome followed by running the gmap aligner:

${TRINITY_HOME}/util/misc/process_GMAP_alignments_gff3_chimeras_ok.pl \
         --genome minigenome.fa \
         --transcripts trinity_out_dir_GG/Trinity-GG.fasta \
         --SAM | samtools view -Sb | samtools sort -o trinity-GG.gmap.bam

Index the bam file and import it into IGV to view alongside the aligned reads and the stringtie transcripts.

How do the Trinity-reconstructed transcripts compare to StringTie?

Genome-free de novo assembly

Next, try assembling the reads directly, without using the genome sequence:

$TRINITY_HOME/Trinity --left reads_1.fq.gz --right reads_2.fq.gz \
           --seqType fq --CPU 2 --max_memory 1G \
	   --output trinity_out_dir

You’ll find the Trinity assembly output file as:

trinity_out_dir/Trinity.fasta

How many transcripts did Trinity reconstruct?

Exploring transcript structures can be more challenging when you do not have a genome sequence to serve as a reference for orienting transcripts and defining intron/exon structures or comparing structures of alternatively spliced variants.

There is a tool called Bandage that allows you to explore your transcriptome assembly, and the various structures that may result from alternative splicing.

Try downloading/installing Bandage, importing your Trinity de novo assembly, and exploring the data. Do you find evidence of alternative splicing within your assembly?

Below is a view of a Trinity assembly from within Bandage.

Bandage_Trinity-GG

Do you find evidence for alternative splicing?

Course