RNA-seq analysis (practice)

    For this RNA-seq analysis, we will use the data from: Cheng et al, MicroRNA silencing for cancer therapy targeted to the tumour microenvironment. Nature 2015. The raw sequencing data is at GEO with six FASTQ files ：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE61851

    Normally for paired-end RNA-seq, each sample will have two separate FASTQ files, with line-by-line correspondence to the two reads from the same fragment. This data is a little unusual in that the two reads of the same fragment are written in the same file. We are going to use the RNA-seq workflow (http://www.bioconductor.org/help/workflows/rnaseqGene/) in BioConductor to analyze the RNA-seq data.

    For RNA-seq read mapping, a popular mapper is TopHat (Trapnell et al, Bioinformatics 2009): https://ccb.jhu.edu/software/tophat/index.shtml. However, if the data is paired-end sequencing with huge read count, STAR is much faster

1.STAR and TopHat : Use STAR (Dobin et al, Bioinformatics 2012)  and TopHat to map the reads to the mouse reference genome.

Build index:

    STAR --runMode genomeGenerate --runThreadN 2–-genomeDir /week_5/index/ --genomeFastaFiles /mnt/Storage/data/BWA/mm9.fa --sjdbGTFfile /mnt/Storage/home/huse/Data/RNA_pipe/mm9_refSeq.gtf --sjdb Overhang 49

Mapping:

    STAR runThreadN 2 --genomeDir RNASEQ/star/index/ --readFilesIn RNASEQ/SRR1592185_1.fastq  SRR1592185_2.fastq --outFileNamePrefix SRR1592185

    tophat -p 2 -G /mnt/Storage/home/huse/Data/RNA_pipe/mm9_refSeq.gtf /mnt/Storage/data/Bowtie/mm9  SRR1592185_1.fastq  SRR1592185_2.fastq > mapping.log

#–runThreadN 30

#–runMode genomeGenerate

#–genomeDir  /home/jmzeng/hoston/mouse/STAR-mouse #put builded indexes into this directory

#–genomeFastaFiles  fasta file of

#–sjdbGTFfile  #annotation file

#–sjdbOverhang 284   #length of reads sequenced - 1

tophat -p 2 -G /mnt/Storage/home/huse/Data/RNA_pipe/mm9_refSeq.gtf /mnt/Storage/data/Bowtie/mm9 GSM1515732_1.fastq GSM1515732_2.fastq > mapping.log

 2. The above test will allow you to learn how to use the read mappers well.After mapping using TopHat , we can get bam files. Use HTSeq to count the reads for each gene from the BAM files. For paired-end sequencing, the reads must be sorted by name. Input files of htseq-count are sam format. So we sort sam files by name and transfrom bam files into sam files. Below are Commands:

samtools sort -n GSM1515732.bam

samtools view -h GSM1515732.bam > GSM1515732.sam

htseq-count -s no -a 10 week_5/GSM1515732.sam RNA_pipe/mm9_refSeq.gtf > GSM1515732_count.txt

3. Use DESeq2 (Love et al, Genome Biol 2014) to look at differential expression. \

source("https://bioconductor.org/biocLite.R")

biocLite("DESeq2")

#count matrix

sampleNames <- c("DOX_1","DOX_2","DOX_3","NODOX_1","NODOX_2","NODOX_3")

mydata <- read.table("counts.txt", header = TRUE, quote = '\t',skip =1)

names(mydata)[7:12] <- sampleNames

countMatrix <- as.matrix(mydata[7:12])

rownames(countMatrix) <-mydata$Geneid

table2 <- data.frame(name = c("DOX_1","DOX_2","DOX_3","NODOX_1","NODOX_2","NODOX_3"),condition = ("DOX","DOX","DOX","NODOX","NODOX","NODOX"))

rownames(table2) <- sampleNames

head(countMatrix)

#transform count matrix into format of DESeq2

dds <- DESeqDataSetFromMatrix(countMatrix, colData=table2, design= ~ condition)

#filter data of count below 0.Because they are useless

dds <- dds[ rowSums(counts(dds)) > 1, ]

#PCA

rld <- rlog(dds)

plotPCA(rld, intgroup=c("name","condition"))

#PCA - ggplot2  note : don't library DESeq before PCA analysis  

library(ggplot2)

rld <- rlog(dds)

data <- plotPCA(rld, intgroup=c("condition", "name"), returnData=TRUE)

percentVar <- round(100 * attr(data, "percentVar"))

p<- ggplot(data, aes(PC1, PC2, color=condition, shape=name)) +

geom_point(size=3) +

xlab(paste0("PC1: ",percentVar[1],"% variance")) +

ylab(paste0("PC2: ",percentVar[2],"% variance"))

#differential expression

library(DESeq)

dds <- DESeq(dds)

res <- results(dds)

write.table(res,"result.csv", sep = ",", row.names = TRUE)

head(res)

#(1)rownames:gene ID (2)baseMean (3)log2FoldChange (4)pvalue (5)padj:adjusted p value

#summary　

summary(res)

#MA plot

library(geneplotter)

plotMA(res, main="DESeq2", ylim=c(-2,2))

plotDispEsts(dds)

#heatmap

sum(res$padj < 0.1, na.rm=TRUE)

library("pheatmap")

select <- order(rowMeans(counts(dds,normalized=TRUE)),decreasing=TRUE)[1:1000]

nt <- normTransform(dds) # defaults to log2(x+1)

log2.norm.counts <- assay(nt)[select,]

df <- as.data.frame(colData(dds)[,c("name","condition")])

pdf('heatmap1000.pdf',width = 6, height = 7)

pheatmap(log2.norm.counts, cluster_rows=TRUE, show_rownames=FALSE,

cluster_cols=TRUE, annotation_col=df)

dev.off()

( To be continued ... )