RNA-seq analysis in R

Gene Set Testing for RNA-seq - Solutions

Challenge 1

Another common way to rank the genes is to order by pvalue, but also, sorting so that upregulated genes are at start and downregulated at the other - you can do this combining the sign of the fold change and the pvalue. 1. Rank the genes by statisical significance - you will need to create a new ranking value using -log10({p value}) * sign({Fold Change}) 2. Load the “C2” pathways from the the data/mouse_c2_v5.RData file 3. Run fgsea using the new ranked genes and the C2 pathways 4. Run fgsea using the new ranked genes and the H pathways. How do these results differ from the ones we got when ranking by the fold change alone?

ranks <- -log10(gseaDat$pvalue) * sign(gseaDat$logFC)
names(ranks) <- gseaDat$Entrez  
load("../Robjects/mouse_c2_v5.RData")
pathwaysC2 <- Mm.c2
fgseaResC2 <- fgsea(pathwaysC2, ranks, minSize=15, maxSize = 500, nperm=1000)

There are ties in the preranked stats (0.04% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

head(fgseaResC2[order(padj, -abs(NES)), ], n=10)

fgseaResH <- fgsea(pathwaysH, ranks, minSize=15, maxSize = 500, nperm=1000)

There are ties in the preranked stats (0.04% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

head(fgseaResH[order(padj, -abs(NES)), ], n=20)

Challenge 2

1. Create a gene list for genes that are up-regulated by at least 4x (logFC>2) in lactating mice
2. Run a goseq analysis on this gene list
3. Plot the results
4. How is this result different to the previous GO analysis?

Look at genes that are up regulated by at least 4x (logFC>2) in lactating mices

isSigGeneUp <- shrinkLvV$FDR < 0.01 &
    !is.na(shrinkLvV$FDR) &
    shrinkLvV$logFC > 2
genesUp <- as.integer(isSigGeneUp)
names(genesUp) <- shrinkLvV$GeneID
pwf <- nullp(genesUp, "mm10", "ensGene", bias.data = shrinkLvV$medianTxLength)

initial point very close to some inequality constraints

goResultsUp <- goseq(pwf, "mm10","ensGene", test.cats=c("GO:BP"))
goResultsUp %>%
    top_n(10, wt=-over_represented_pvalue) %>%
    mutate(hitsPerc=numDEInCat*100/numInCat) %>% 
    ggplot(aes(x=hitsPerc, 
               y=term, 
               colour=over_represented_pvalue, 
               size=numDEInCat)) +
        geom_point() +
        expand_limits(x=0) +
        labs(x="Hits (%)", y="GO term", colour="p value", size="Count")

Challenge 3

1. Use pathview to export a figure for “mmu04060”, but this time only use genes that are statistically significant at FDR < 0.01

sigGenes <- shrinkLvV$FDR < 0.01 & !is.na(shrinkLvV$FDR)
logFC <- annotLvV$logFC[sigGenes]
names(logFC) <- annotLvV$Entrez[sigGenes]
pathview(gene.data = logFC, 
         pathway.id = "mmu04060", 
         species = "mmu", 
         limit = list(gene=5, cpd=1))

Info: Downloading xml files for mmu04060, 1/1 pathways..
Info: Downloading png files for mmu04060, 1/1 pathways..
'select()' returned 1:1 mapping between keys and columns
Info: Working in directory /Users/sawle01/Documents/RNAseq_course/CRUK_CI_RNAseq_course/Course_Materials/solutions
Info: Writing image file mmu04060.pathview.png

mmu04060.pathview.png:

mmu04060 -Cytokine-cytokine receptor interaction