Welcome to the blog

Posts

My thoughts and ideas

Team Assignment - Expression and DE | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Team Assignment - Expression and DE

Course

In the previous exercise, teams have aligned their RNAseq data and performed QC evaluations. Using their aligned data, students will now apply the concepts they have learned regarding expression estimation and differential expression analysis. To complete this assignment, students will need to review commands we performed in earlier sections.

Before starting this team exercise, first find the folder containing your 6 aligned bam files (along with the index files). Note: In the previous exercise, you merged bams files for easy visualization in IGV, we will not be using that for expression and de analysis.

Expression Estimation

Goals:

  • Familiarize yourself with Stringtie options
  • Run Stringtie to obtain expression values
  • Run provided stringtie helper perl script to combine results into a single file

Teams can now use Stringtie to estimate the gene expression levels in their sample and answer the following questions:

### Remember to do this in a new directory under team_exercises
mkdir -p $RNA_HOME/team_exercise/expression/stringtie/ref_only
cd $RNA_HOME/team_exercise/expression/stringtie/ref_only

<L4> Q1. Based on your stringtie results, what are the top 5 genes with highest average expression levels across all knockout samples? What about in your rescue samples? How large is the overlap between the two sets of genes? (Hint: You can use R for this analysis)

Differential Expression

Goals:

  • Perform differential analysis between the knockout and rescued samples
  • Check which genes are differentially expressed with statistical significance
  • Visualize DE results

Teams will now use ballgown to perform differential analysis followed by visualization of their results.

<L3> Q2. Follow through the ballgown differential expression section by making modifications using your respective sample names. Hint: You will need to create a separate directory under your team_exercises folder for your ballgown outputs. You will also need to change the respective sample names and paths following the printf command.

<L3> Q3. How many significant differentially expressed genes do you observe?

<L5> Q4. By referring back to the supplementary tutorial in the DE Visualization Module, can you construct a heatmap showcasing the significantly de genes? Try playing around with the de filter to include more/less genes in your heatmap. Try to determine the best cutoff for your specific dataset.

<L5> (OPTIONAL) Q5. Pick one of the significantly differentially expressed genes and visualize gene expression levels across the 6 samples as well as individual transcript expression levels for those corresponding to your gene of interest. (Hint: How can you modify the transcript expression plot in the DE Visualization section to showcase gene expression levels instead of transcript expression levels?)

Additionally, students should feel free to explore other visualization methods, including those they may have used in past research experiences and share with the class.

OPTIONAL: Bonus questions

  • Run htseq to get raw counts and then use edgeR for differential expression
  • Compare results between ballgown de and edgeR

<L4> Q6. After obtaining your edgeR results, how does it agree with your previously obtained de results using ballgown?

Presenting Your Results

At the end of this team exercise, students will show how they visualized their differential expression results to the class.

Course
Team Assignment - Alignment | Griffith Lab

RNA-seq Bioinformatics

Introduction to bioinformatics for RNA sequence analysis

Team Assignment - Alignment

Course

The goal of the following team assignment is for students to gain hands-on experience by working on recently published RNA-seq data and apply the concepts they have learned up to RNA alignment. To complete this assignment, students will need to review commands we performed in earlier sections.

Background on Dataset used

In this assignment, we will be using subsets of the GSE136366 dataset (Roczniak-Ferguson A, Ferguson SM. Pleiotropic requirements for human TDP-43 in the regulation of cell and organelle homeostasis. Life Sci Alliance 2019 Oct;2(5). PMID: 31527135). This dataset consists of 6 RNA sequencing files of human cells that either express or lack the TDP-43 protein.

Experimental Details

Obtaining the dataset & reference files

Goals:

As mentioned previously, we have subsetted the 6 RNA-seq samples into 5 different chromosome regions. Each team can download their corresponding dataset using the following commands.

mkdir -p $RNA_HOME/team_exercise/data
cd $RNA_HOME/team_exercise/data


### TEAM A
wget -c http://genomedata.org/seq-tec-workshop/read_data/rna_alignment-de_exercise/dataset_A/dataset.tar.gz
tar -xzvf dataset.tar.gz

### TEAM B
wget -c http://genomedata.org/seq-tec-workshop/read_data/rna_alignment-de_exercise/dataset_B/dataset.tar.gz
tar -xzvf dataset.tar.gz

### TEAM C
wget -c http://genomedata.org/seq-tec-workshop/read_data/rna_alignment-de_exercise/dataset_C/dataset.tar.gz
tar -xzvf dataset.tar.gz

### TEAM D
wget -c http://genomedata.org/seq-tec-workshop/read_data/rna_alignment-de_exercise/dataset_D/dataset.tar.gz
tar -xzvf dataset.tar.gz

### TEAM E
wget -c http://genomedata.org/seq-tec-workshop/read_data/rna_alignment-de_exercise/dataset_E/dataset.tar.gz
tar -xzvf dataset.tar.gz

Additionally, teams will need to create a separate directory and download the corresponding reference files needed for RNA alignment & further expression analysis:

mkdir -p $RNA_HOME/team_exercise/references
cd $RNA_HOME/team_exercise/references

## Adapter trimming
wget -c http://genomedata.org/seq-tec-workshop/references/RNA/illumina_multiplex.fa

## Reference fasta corresponding to your team's assigned chromosome (e.g. chr6)
wget -c http://genomedata.org/seq-tec-workshop/references/RNA/chr6.fa

## Obtain annotated reference gtf file corresponding to your team's assigned chromosome (e.g. chr6)
wget -c http://genomedata.org/seq-tec-workshop/references/RNA/chr6_Homo_sapiens.GRCh38.95.gtf

Data Preprocessing (QC & Trimming)

Goals:

Prior to aligning RNA-seq data, teams should perform adapter trimming using flexbar. Once the team has both the pre-trim and post-trim data, QC metrics should be evaluated using fastqc and a general report can be generated using multiqc.

<L1> Q1. What is the average percentage of reads that are trimmed?

<L2> Q2. Before looking at the multiqc report, how do you expect the sequence length distribution to look both prior to and after trimming? Is your answer confirmed by the multiqc report results?

<L1> Q3. Are there any metrics where the sample(s) failed?

Alignment Exercise

Goals:

A useful option to add to the end of your commands is 2>, which redirects the stderr from any command into a specific file. This can be used to redirect your stderr 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.

<L1> Q4. What were the percentages of reads that aligned to the reference for each sample?

<L1> Q5. By compressing your sam format to bam, approximately how much space is saved (fold change in size)?

Post-alignment QC & IGV Visualization

Goals:

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 in IGV.

# Merge the bams for visualization purposes
cd <path to dir with alignments>
java -Xmx16g -jar $PICARD MergeSamFiles OUTPUT=KO_merged.bam INPUT=SRR10045016.bam INPUT=SRR10045017.bam INPUT=SRR10045018.bam
java -Xmx16g -jar $PICARD MergeSamFiles OUTPUT=RESCUE_merged.bam INPUT=SRR10045019.bam INPUT=SRR10045020.bam INPUT=SRR10045021.bam

<L2> Q6. How does the information from your post-alignment QC report differ from pre-alignment QC?

<L3> Q7. IGV: Can you identify certain exons that have significantly more/less coverage in one of your KO/RESCUE samples compared to the other? What is happening here?

<L3> Q8. IGV: Can you identify regions where the RNAseq reads are mapping to unexpected regions? What do you think is the reason for this phenomenon?

<L3> Q9. IGV: Can you identify a gene region that has RNA sequencing support for multiple isoforms?

BONUS QUESTION for after you obtain your alignments

Upon obtaining the different reference files, explore the annotated reference gtf file and answer the following questions using your choice of commands. Hint: useful commands include cat, grep, cut, sort, uniq, awk

<L4> Q10. What are the different types of features contained in the gtf file (e.g. transcript, gene)? What are the frequencies of the different types of features? (This is referring to the third field/column of the data).

In order to get this answer, there are a series of commands that we can pipe together: cat <YOUR GTF FILE> | grep gene_name | cut -d$'\t' -f3 | sort | uniq -c | sort -r

Can you explain how this command works to one of the TAs?

<L5> Q11. Now that you have seen the example in Q1, can you construct a similar command that answers the questions: Which genes have the highest number of transcripts (either gene id or gene name)? How many?

Presenting Your Results

At the end of this team exercise, groups will present findings from their QC reports and IGV analysis to the class for specific questions listed below.

Team A: Present IGV findings regarding question 9.

Team B: Present multiqc report on pre- and post-alignment qc files (question 6).

Team C: Present IGV findings regarding question 7.

Team D: Present IGV findings regarding question 8.

Team E: Present multiqc report on pre- and post-trimming qc files (Data preprocessing section).

Course