NGS RNA-seq exploration short course

This brief tutorial will walk you through data analysis of an RNA-seq experiment.

In this experiment, E. coli was inoculated into culture and the culture was then sampled at 4 hours and 24 hours post inoculation.  The experiment was run in triplicate.

RNA was extracted from the 6 samples, fragmented, and sequenced.   All sequencing runs were of the paired-end 2x100 type, so each RNA fragment is read from both ends, 100 bp from each end.

Here is a table showing the data we have:

In class demo

In class, we will explore and characterize the raw data.  

cd scratch
mkdir rawdata
cd rawdata;
ln -s /home/scott/e0/* . ;
cd ..;
mkdir finaldata;
cd finaldata;
ln -s /home/scott/e3/* .;
cd ..;

Here are some elements (programs & techniques) we may use (you will need some of these for the homework):

 gunzip -c MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz | egrep [ACGT]\{101\} | grep -v GATCGGAAGAGCACAC | head -10000 | sort | uniq -c | sort -n -r | head

fastqc MURI_102_SA14008_CAGATC_L005_R1_001.fastq.gz

gunzip -c MURI_102_SA14008_CAGATC_L006_R1_001.fastq.gz | head -400000 > r1.fq
fastqc r1.fq

bwa index REL606.fna

bwa mem -t 16 REL606.fna MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz MURI_21_SA13172_GTGGCC_L007_R2_001.fastq.gz | awk '{if ($3!="*") {print $0}}' | samtools view -b -S - > MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.bam

samtools sort MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.bam MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted
samtools index MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam
samtools flagstat MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam > MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam.flagstat
samtools idxstats MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam > MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam.idxstats

samtools mpileup -Euf $3 $1.sorted.bam | bcftools view -bvcg - > $1.sorted.bam.raw.bcf &

#!/bin/bash
# $1 = input file base name fwd, 
# $2 = input file base name rev,
# $3 = reference file, already bwa indexed
echo "$0: Starting on `date`"
echo "BWA mem"
bwa mem -t 16 $3 $1 $2 | awk '{if ($3!="*") {print $0}}' | samtools view -b -S - > $1.bam
# echo "Samtools view: convert to binary"
# samtools view -b -S $1.sam > $1.bam
# rm $1.sam
echo "Samtools sort: sorting bam file"
samtools sort $1.bam $1.sorted
# rm $1.bam
echo "Samtools index: creating index of sorted bam file"
samtools index $1.sorted.bam
samtools flagstat $1.sorted.bam > $1.sorted.bam.flagstat
samtools idxstats $1.sorted.bam > $1.sorted.bam.idxstats
# echo "Samtools mpileup and bcftools"
samtools mpileup -Euf $3 $1.sorted.bam | bcftools view -bvcg - > $1.sorted.bam.raw.bcf &
echo "$0: Done on `date`"

bedtools coverage -abam  MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam -b REL606.gff > MURI_21_SA13172_GTGGCC_L007_R1_001.fastq.gz.sorted.bam.bedtools.out

for files in `ls *.sorted.bam`; do bedtools coverage -abam $files  -b REL606.gff > $files.bedtools.out; done

for file in `ls *.bedtools.out`; do cat $file | grep CDS | awk 'BEGIN {FS="\t"} {print $10}' > $file.counts ; done

cat MURI_17_SA13172_ATGTCA_L007_R1_001.fastq.gz.sorted.bam.bedtools.out | grep CDS | cut -f 9 | awk 'BEGIN {FS=";"} {print $3":"$4}' > genes.txt
paste genes.txt MURI_17*counts MURI_26*counts MURI_98*counts MURI_21*counts MURI_30*counts MURI_102*counts > all3x3.counts

cd finaldata
R

library(ggbiplot)
wall<-read.table(file="all3x3.counts",sep="\t",header=FALSE); load the raw data into variable "wall"
wallt<-t(wall[,2:7]); # This just transposes the read count data into a new variable, wallt

wallt.pca<-prcomp(wallt)
print(ggbiplot(wallt.pca, groups=c("4hr","4hr","4hr","24hr","24hr","24hr"), ellipse = TRUE, circle=TRUE, obs.scale = 1, var.scale = 1, var.axes=FALSE))

boxplot(wall)

boxplot(log(wall[,2:7]))

Homework

For your homework, you will investigate the validity of combining data files from different sequencing runs.  Only a few of these questions require working at a computer keyboard, but I encourage you to work in groups to solve the entire set of questions.

1. Based on what you learned about the T-test (that is, using terms associated with a T-test), explain what criteria you might use to consider it "invalid" to combine the multiple raw sequence data files from samples.  (5 points)
2. Outline the steps needed to reduce the raw data to numbers suitable for evaluation of your criteria in question #1. (5 points)
3. Perform the steps you outlined in #2 and tell whether or not it was valid to combine the data files. (20 points)
4. Starting with the raw "count" data, explore the effect on PCA of NOT normalizing.  Turn in a print out of the new PCA plot. (10 points)
5. Although we did not explore this in class, DNA mutations were automatically tallied during our mapping process.  These results are in the files ending in ".bcf".  Using the tool "bedtools" to view these results, test the hypothesis that transitions are more common than transversions.  Support your answer with data from this experiment. (20 points)
6. Continuing with the mutation analysis, examine whether the mutation frequency in this sample set differs between protein coding and non-protein coding regions of the genome. Support your answer with data from this experiment. (30 points).

Email your answers/PDFs to shunickesmith <at> gmail.com, cc: Prof. Matouschek no later than TBD, 10:00 am (BETTER: before Thanksgiving break).

Homework - Revised 9:00 pm Monday 11/24/14  Revised 5:00 pm Thursday 11/20/14

For your homework, you will investigate the validity of combining data files from different sequencing runs.  Only a few of these questions require working at a computer keyboard, but I encourage you to work in groups to solve the entire set of questions.

Before you begin, pull up this web site in a new window - it's an all-class, live group chat to which you can post questions, get answers and even answer questions of your fellow students.  Dr. Hunicke-Smith and Benni will be monitoring it periodically (largely during the daytime and early evenings).

By next Tuesday, 12/2 11/25, 10:00 am do the following:

1. Follow the steps listed above on this web page under the expanding section "This is your homework due 11/25..." to log into appsoma and setup access to the data you will need from here on.
2. Move into the "rawdata" directory, find the first four lines of the read 1 sequence file for the MURI 102 sample from lane 6 of sequencing run SA14008 - put them into a new file in that directory called "s1.fq" and copy it into an email.
3. Move into the "finaldata" directory and make sure you can see the gene expression data file "all3x3.counts"
4. Using the linux "sort" command, sort all3x3.counts 6 times, sorting on the expression values of each of the 6 samples separately from lowest to highest, redirecting the output of each sort into a separate file.
5. Using the linux command "tail -1" on each of these 6 files, copy the name of the most abundant gene from each sample into the same email.

Remember - use Etherpad to ask questions and get answers!  Email your answers to shunickesmith <at> gmail.com, cc: Prof. Matouschek no later than 11/25, 10:00 am.