Welcome to the blog


My thoughts and ideas

Proposed Improvements | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Proposed Improvements

Install docs improvements

  • consolidate all install/configuration documentation to the AWS setup page. Change the student Tool Installation section to be a more didactic section. Demonstrate different kinds of tool installations and have some exercises but it isn’t about the students install all the tools they need.

DE analysis improvements

  • incorporate ideas from this tutorial and make note of it in the resources: https://dockflow.org/workflow/rnaseq-gene/#content
  • New paper comparing many RNA-seq tools: https://www.nature.com/articles/s41467-017-00050-4.epdf
  • Great paper on number of replicates to use (incorporate into lecture): http://rnajournal.cshlp.org/content/22/6/839.long

Update data sets on the MGI FTP site

  • Use a larger data set that has not been so heavily downsampled
  • Put the original UHR and HBR instrument data on the FTP site
  • Put the instrument data for the new hcc1395 data on the FTP site
  • MGI notes on this data will be here: https://confluence.gsc.wustl.edu/display/CI/Cancer+Informatics+Test+Data

More independent exercises, group exercises, and integrated assignments

  • Each module should have at least two exercises where the students are not copying/pasting anything. One could be at 1/2 way point, and the other at the end of the module. The one at the end could be a group exercise.
  • Each day should end with at least one hour of integrated assignment.

Install ‘pip’ command into the AMI

In order for htseq-count to use bam files directly it needs pysam. This can be installed with pip but that is not available by default.

Then, On AMI install and test htseq-count with bam files:

sudo apt-get install python-pip
sudo pip install pysam

An alternative install procedure that has been tested and worked is as follows. The above procedure is preferred and should be tried first.

cd ~/bin/
wget https://pysam.googlecode.com/files/pysam-0.7.5.tar.gz
tar -zxvf pysam-0.7.5.tar.gz
cd pysam-0.7.5/
python setup.py build
sudo python setup.py install

Get X11 support working on AMI

For R and other applications it would be nice if X11 worked. Note the install instructions for R would need to change as well.

Create a trimming section

Create a wiki section and exercise that summarizes read trimming concepts. Start with some raw data, including aligned reads. Align these reads without any trimming and assess alignment statistics using Picard, FastQC, etc. Now take these same reads and perform both adaptor trimming and quality trimming. Re-align the trimmed reads and assess the effect of trimming on alignment metrics.

Note: we do have an okay trimming section now, but it should really be starting with raw reads, rather than those that had already been proven to align to chr 22. This makes it harder to see the value of trimming by comparing alignment pre/post trimming. Instead we should perhaps get chr22 matching reads by k-mer analysis (kallisto), then trim (adapter and/or quality), then align trimmed and untrimmed, and summarize the difference.

Create a batch effect section

We should add a section about batch effects. Both detecting the presence of batch effects as well as correcting for them during analysis.

Use the Snyder lab mouse/human tissue expression data as an example?

Add documentation to detailed cloud tutorial to provide some better security practices

For convenience the cloud instances have been set up with very permissive security. Some better practices should be documented.

Add a fusion detection section

We previously had a fusion detection module but it was difficult to complete in time frames appropriate for a workshop. Further optimization is required. Another challenge is the lack of well engineered fusion detection software. This publication State-of-the-art fusion-finder algorithms sensitivity and specificity does a decent job of summarizing the current options available. Another caveat of this topic is that is mostly of interest to cancer researchers so it might only be included where there are sufficient students with this interest.

Improve alignment QC section

In particular we should add use of Picard CollectRnaSeqMetrics (https://broadinstitute.github.io/picard/command-line-overview.html) and RNA-SeQC (http://www.broadinstitute.org/cancer/cga/rna-seqc). It would also be good to include use of splicing metrics calculated from the TopHat junctions files. A standalone version of the TGI tool that does this would need to be created for this purpose.

Also use MultiQC to produce a combined report of the QC results: http://multiqc.info/

Improve Expression/Differential expression lectures

There a some nice slides/concepts that we could borrow from the BaseSpace Demo slides (see Obis ~/Dropbox Teaching/CSHL/2015/Workshop-CSHL-RNA-Seq-Metagenomics.pdf).

Create a genome reference free analysis module.

  • Expand on the current Kallisto exercise.
  • Add Sleuth for DE
  • Consider adding Salmon as an alternate to Kallisto?

Create a new gene fusion discovery module

  • Try Pizzly

Add more content on downstream analysis

For example, pathway analysis of RNA-seq data, clustering, etc.

Identify more interesting data sets to use for the alternative splicing module


Gray lab breast cancer cell line dataset:

  • http://www.ncbi.nlm.nih.gov/pubmed/24176112
  • https://www.biostars.org/p/111040/ (biostars tutorial on downloading data)
  • https://github.com/genome/gms/wiki/Guide-to-Importing-and-Analyzing-External-Data (another guide on downloading and reformatting this data)

Create lecture materials that introduce k-mer based approaches

  • Intro to Kallisto, Sailfish, Salmon
  • Integrate Kallisto content into the expression lecture and expression modules of the tutorial. That is where it makes sense anyway. Drop module 5.

Update the tutorial to take into account recent developments in RNA-seq analysis methods, best practices, and new tools

Update the alternative splicing module to provide an alternative analysis workflow to StringTie/Ballgown

Introduction to RegTools functionality. After HISAT2 alignment, QoRTs (written in Scala) performs QC but also processes RNA-seq data to produce count files needed for splicing analysis by DEXSeq and JunctionSeq.

Integrated Assignment Answers | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Integrated Assignment Answers

Integrated Assignment answers

Background: The use of cell lines are often implemented in order to study different experimental conditions. One such kind of study is the effects of shRNA on expression profiles, to determine whether these effects target specific genes. Experimental models for these include using control shRNA to account for any expression changes that may occur from just the introduction of these molecules.

Objectives: In this assignment, we will be using a subset of the GSE114360 dataset, which consists of 6 RNA sequence files on the SGC-7901 gastric cancer cell line, (3 transfected with tcons_00001221 shRNA, and 3 control shRNA), and determine the number of differentially expressed genes.

Experimental information and other things to keep in mind:

##PART 0 : Obtaining Data and References


Create a working directory ~/workspace/rnaseq/integrated_assignment/ to store this exercise. Then create a unix environment variable named RNA_ASSIGNMENT that stores this path for convenience in later commands.

export RNA_HOME=~/workspace/rnaseq
mkdir -p ~/workspace/rnaseq/integrated_assignment/
export RNA_ASSIGNMENT=~/workspace/rnaseq/integrated_assignment/

You will also need the following environment variables througout the assignment:

export RNA_REF_INDEX=$RNA_REFS_DIR/Homo_sapiens.GRCh38
export RNA_REF_FASTA=$RNA_REF_INDEX.dna.primary_assembly.fa
export RNA_REF_GTF=$RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf

Obtain reference, annotation, adapter and data files and place them in the integrated assignment directory Note: when initiating an environment variable, we do not need the $; however, everytime we call the variable, it needs to be preceeded by a $.

ln -s ~/CourseData/RNA_data/Integrative_Assignment/reference/
ln -s ~/CourseData/RNA_data/Integrative_Assignment/raw_reads/top_1mil/ raw_reads
ln -s ~/CourseData/RNA_data/Integrative_Assignment/adapter

Q1.) How many items are there under the “reference” directory (counting all files in all sub-directories)? What if this reference file was not provided for you - how would you obtain/create a reference genome fasta file. How about the GTF transcripts file from Ensembl?

A1.) The answer is 19. Review these files so that you are familiar with them. If the reference fasta or gtf was not provided, you could obtain them from the Ensembl website under their downloads > databases.

cd $RNA_ASSIGNMENT/reference/
find *
find * | wc -l

Q2.) How many exons does the gene SOX4 have? How about the longest isoform of PCA3?

A2.) SOX4 only has 1 exon, while the longest isoform of PCA3 has 7 exons. Review the GTF file so that you are familiar with it. What downstream steps will we need this gtf file for?

cd $RNA_ASSIGNMENT/reference/gtf
grep -w "SOX4" Homo_sapiens.GRCh38.86.chr9.gtf
grep -w "PCA3" Homo_sapiens.GRCh38.86.chr9.gtf | grep "exon_number"

Q3.) How many samples do you see under the data directory?

A3.) The answer is 6. The samples are paired per file, and are named based on their accession number.

cd $RNA_ASSIGNMENT/raw_reads/
ls -l
ls -1 | wc -l

NOTE: The fastq files you have copied above contain only the first 1000000 reads. Keep this in mind when you are combing through the results of the differential expression analysis.

##Part 1 : Data preprocessing


Now create a new folder that will house the outputs from FastQC. Use the -h option to view the potential output on the data to determine the quality of the data.

mkdir raw_fastqc
fastqc $RNA_DATA_DIR/* -o raw_fastqc/

Q4.) What metrics, if any, have the samples failed? Are the errors related?

A4.) The per base sequence content of the samples don’t show a flat distribution and do have a bias towards certain bases at particular positions. The reason for this is the presense of adapters in the reads, which also shows a warning if not a failure in the html output summary.

Now based on the output of the html summary, proceed to clean up the reads and rerun fastqc to see if an improvement can be made to the data. Make sure to create a directory to hold any processed reads you may create.

mkdir trimmed_reads
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RNA_ILL_ADAPT/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155055_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155055_2.fastq.gz --target trimmed_reads/SRR7155055
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RNA_ILL_ADAPT/adapter/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155056_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155056_2.fastq.gz --target trimmed_reads/SRR7155056
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RNA_ILL_ADAPT/adapter/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155057_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155057_2.fastq.gz --target trimmed_reads/SRR7155057
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RNA_ILL_ADAPT/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155058_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155058_2.fastq.gz --target trimmed_reads/SRR7155058
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RNA_ILL_ADAPT/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155059_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155059_2.fastq.gz --target trimmed_reads/SRR7155059
flexbar --adapter-min-overlap 7 --adapter-trim-end RIGHT --adapters $RRNA_ILL_ADAPT/illumina_multiplex.fa --pre-trim-left 13 --max-uncalled 300 --min-read-length 25 --threads 8 --zip-output GZ --reads $RNA_DATA_DIR/SRR7155060_1.fastq.gz --reads2 $RNA_DATA_DIR/SRR7155060_2.fastq.gz --target trimmed_reads/SRR7155060

Q5.) What average percentage of reads remain after adapter trimming? Why do reads get tossed out?

A5.) Around 99% of reads still survive after adapter trimming. The reads that get tossed are due to being too short after trimming (they fall below our threshold of minimum read length of 25.

Q6.) What sample has the largest number of reads after trimming?

A6.) The control sample 2 (SRR7155059) has the most reads (1999678/2 = reads). An easy way to figure out the number of reads is to check the output log file from the trimming output. Looking at the “remaining reads” row, we see the reads (each read in a pair counted individually) that survive the trimming. Alternatively, you can make use of the command ‘wc’. This command counts the number of lines in a file. Since fastq files have 4 lines per read, the total number of lines must be divided by 4.

zcat $RNA_ASSIGNMENT/trimmed_reads/SRR7155059_1.fastq.gz | head -n 4
@SRR7155059.1 1 length=150
Running this command only give you the total nunmber of lines in the fastq file (Note that because the data is compressed, we need to use zcat to unzip it and print it to the screen, before passing it on to the wc command):
$RNA_ASSIGNMENT/trimmed_reads/SRR7155059_1.fastq.gz | wc -l

##PART 2: Data alignment


A useful option to add to the end of your commands is 2>, which redirects the stdout from any command into a specific file. This can be used to redirect your stdout into a summary file, and can be used as follows: My_alignment_script 2> alignment_metrics.txt. The advantage of this is being able to view the alignment metrics later on.

To create HISAT2 alignment commands for all of the six samples and run alignments:

mkdir -p $RNA_ALIGN_DIR
hisat2 -p 8 --rg-id=Transfect1 --rg SM:Transfect --rg LB:Transfect1_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155055_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155055_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155055.sam
hisat2 -p 8 --rg-id=Transfect2 --rg SM:Transfect --rg LB:Transfect2_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155056_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155056_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155056.sam
hisat2 -p 8 --rg-id=Transfect3 --rg SM:Transfect --rg LB:Transfect3_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155057_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155057_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155057.sam
hisat2 -p 8 --rg-id=Control1 --rg SM:Control --rg LB:Control1_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155058_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155058_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155058.sam
hisat2 -p 8 --rg-id=Control2 --rg SM:Control --rg LB:Control2_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155059_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155059_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155059.sam
hisat2 -p 8 --rg-id=Control3 --rg SM:Control --rg LB:Control3_sub --rg PL:ILLUMINA -x $RNA_REFS_DIR/Homo_sapiens.GRCh38 --dta --rna-strandness RF -1 $RNA_ASSIGNMENT/trimmed_reads/SRR7155060_1.fastq.gz -2 $RNA_ASSIGNMENT/trimmed_reads/SRR7155060_2.fastq.gz -S $RNA_ALIGN_DIR/SRR7155060.sam
#convert sam alignments to bam..how much space did you save by performing this conversion?

samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155055.bam $RNA_ALIGN_DIR/SRR7155055.sam
samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155056.bam $RNA_ALIGN_DIR/SRR7155056.sam 
samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155057.bam $RNA_ALIGN_DIR/SRR7155057.sam 
samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155058.bam $RNA_ALIGN_DIR/SRR7155058.sam
samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155059.bam $RNA_ALIGN_DIR/SRR7155059.sam
samtools sort -@ 8 -o $RNA_ALIGN_DIR/SRR7155060.bam $RNA_ALIGN_DIR/SRR7155060.sam

Q7.) How else could you obtain summary statistics for each aligned file?

A7.) There are many RNA-seq QC tools available that can provide you with detailed information about the quality of the aligned sample (e.g. FastQC and RSeQC). However, for a simple summary of aligned reads counts you can use samtools flagstat. You can also look for the logs generated by TopHat. These logs provide a summary of the aligned reads.

samtools flagstat SRR7155055.bam > SRR7155055.flagstat.txt
samtools flagstat SRR7155056.bam > SRR7155056.flagstat.txt
samtools flagstat SRR7155057.bam > SRR7155057.flagstat.txt
samtools flagstat SRR7155058.bam > SRR7155058.flagstat.txt
samtools flagstat SRR7155059.bam > SRR7155059.flagstat.txt
samtools flagstat SRR7155060.bam > SRR7155060.flagstat.txt
grep "mapped (" *.flagstat.txt

Q8.) Approximatly how much space is saved by converting the sam to a bam format?

A8.) We get about a 5.5x compression by using the bam format instead of the sam format. This can be seen by adding the -lh option when listing the files in the aligntments directory.

ls -lh $RNA_ALIGN_DIR/

In order to make visualization easier, we’re going to merge each of our bams into one using the following commands. Make sure to index these bams afterwards to be able to view them on IGV.

#merge the bams for visulization purposes
java -Xmx2g -jar ~/CourseData/RNA_data/Integrative_Assignment/picard.jar MergeSamFiles OUTPUT=transfected.bam INPUT=SRR7155055.bam INPUT=SRR7155056.bam INPUT=SRR7155057.bam
java -Xmx2g -jar /usr/local/picard/picard.jar MergeSamFiles OUTPUT=control.bam INPUT=SRR7155058.bam INPUT=SRR7155059.bam INPUT=SRR7155060.bam

Try viewing genes such as TP53 to get a sense of how the data is aligned. To do this:

Q9.) What portion of the gene do the reads seem to be piling up on? What would be different if we were viewing whole-genome sequencing data?

A9.) The reads all pile up on the exonic regions of the gene since we’re dealing with RNA-Sequencing data. Not all exons have equal coverage, and this is due to different isoforms of the gene being sequenced. If the data was from a whole-genome experiment, we would ideally expect to see equal coverage across the whole gene length.

Right-click in the middle of the page, and click on “Expanded” to view the reads more easily.

Q10.) What are the lines connecting the reads trying to convey?

A10.) The lines show a connected read, where one part of the read begins mapping to one exon, while the other part maps to the next exon. This is important in RNA-Sequencing alignment as aligners must be aware to take this partial alignment strategy into account.

##PART 3: Expression Estimation


mkdir -p $RNA_ASSIGNMENT/expression
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/transfect1/transcripts.gtf -A Tumor1/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155055.bam
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/transfect2/transcripts.gtf -A Tumor2/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155056.bam
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/transfect3/transcripts.gtf -A Tumor3/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155057.bam
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/control1/transcripts.gtf -A Normal1/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155058.bam
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/control2/transcripts.gtf -A Normal2/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155059.bam
stringtie -p 8 -G $RNA_REFS_DIR/Homo_sapiens.GRCh38.92.gtf -e -B -o $RNA_ASSIGNMENT/expression/control3/transcripts.gtf -A Normal3/gene_abundances.tsv $RNA_ALIGN_DIR/SRR7155060.bam

Q11.) How do you get the expression of the gene SOX4 across the transfect and control samples?

A11.) To look for the expression value of a specific gene, you can use the command ‘grep’ followed by the gene name and the path to the expression file

grep ENSG00000124766 $RNA_ASSIGNMENT/expression/*/transcripts.gtf | cut -f1,9 | grep FPKM

##PART 4: Differential Expression Analysis


mkdir -p $RNA_ASSIGNMENT/ballgown/
cd $RNA_ASSIGNMENT/ballgown/

Perform transfect vs. control comparison, using all samples, for known (reference only mode) transcripts: First create a file that lists our 6 expression files, then view that file, then start an R session where we will examine these results:

printf "\"ids\",\"type\",\"path\"\n\"transfect1\",\"Transfect\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/transfect1\"\n\"transfect2\",\"Transfect\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/transfect2\"\n\"transfect3\",\"Transfect\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/transfect3\"\n\"control1\",\"Control\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/control1\"\n\"control2\",\"Control\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/control2\"\n\"control3\",\"Control\",\"/home/ubuntu/workspace/rnaseq/integrated_assignment/expression/control3\"\n" > Transfect_vs_Control.csv

cat Transfect_vs_Control.csv


*Adapt the R tutorial file that has been provided in the github repo for part 1 of the tutorial: Tutorial_Part1_ballgown.R. Modify it to fit the goals of this assignment then run it.

Q12.) Are there any significant differentially expressed genes? How many in total do you see? If we expected SOX4 to be differentially expressed, why don’t we see it in this case?

A12.) Yes, there are about 523 significantly differntially expressed genes. Due to the fact that we’re using a subset of the fully sequenced library for each sample, the SOX4 signal is not significant at the adjusted p-value level. You can try re-running the above exercise on your own by using all the reads from each sample in the original data set, which will give you greater resolution of the expression of each gene to build mean and variance estimates for eacch gene’s expression.

Q13.) What plots can you generate to help you visualize this gene expression profile

A13.) The CummerBund package provides a wide variety of plots that can be used to visualize a gene’s expression profile or genes that are differentially expressed. Some of these plots include heatmaps, boxplots, and volcano plots. Alternatively you can use custom plots using ggplot2 command or base R plotting commands such as those provided in the supplementary tutorials. Start with something very simple such as a scatter plot of transfect vs. control FPKM values.