RNA-seq Analysis in R

Annotation and Visualisation of RNA-seq results - Solutions

# load data
load("../Robjects/DE.RData")
## set up connection to ensembl database
ensembl=useMart("ENSEMBL_MART_ENSEMBL")
# list the available datasets (species)
listDatasets(ensembl) %>% 
    filter(str_detect(description, "Mouse"))
# specify a data set to use
ensembl = useDataset("mmusculus_gene_ensembl", mart=ensembl)
# Set the filter type and values
filterType <- "ensembl_gene_id"
# annotated
load("../Robjects/Ensembl_annotations.RData")
annotLvV <- as.data.frame(resLvV) %>% 
    rownames_to_column("GeneID") %>% 
    left_join(ensemblAnnot, "GeneID") %>% 
    rename(logFC=log2FoldChange, FDR=padj)
# txdb object
txMm <- makeTxDbFromBiomart(dataset="mmusculus_gene_ensembl")

Challenge 1

That was just 1000 genes. We need annotations for the entire results table. Also, there may be some other interesting columns in BioMart that we wish to retrieve.

1. Search the attributes and add the following to our list of attributes:
   1. The gene description
      
   2. The genomic position - chromosome, start, end, and strand (4 columns)
   3. The gene biotype
      
2. Query BioMart using all of the genes in our results table (resLvV)
   
3. How many Ensembl genes have multipe Entrez IDs associated with them?
   
4. How many Ensembl genes in resLvV don’t have any annotation? Why is this?

filterValues <- rownames(resLvV)

# check the available "attributes" - things you can retreive
listAttributes(ensembl) %>%
    head(20)
attributeNames <- c('ensembl_gene_id',
                    'entrezgene',
                    'external_gene_name',
                    'description',
                    'gene_biotype')

# run the query
annot <- getBM(attributes=attributeNames,
               filters = ourFilterType,
               values = filterValues,
               mart = ensembl)

# dulicate ids
sum(duplicated(annot$ensembl_gene_id))

# missing gens
missingGenes <- !rownames(resLvV)%in%annot$ensembl_gene_id
rownames(resLvV)[missingGenes]

Challenge 2

Use the log2 fold change (logFC) on the x-axis, and use -log10(FDR) on the y-axis. (This >-log10 transformation is commonly used for p-values as it means that more significant genes have a >higher scale)

1. Create a column of -log10(FDR) values

2. Create a plot with points coloured by if FDR < 0.05

# first remove the filtered genes (FDR=NA) and create a -log10(FDR) column
filtTab <- shrinkLvV %>% 
    filter(!is.na(FDR)) %>% 
    mutate(`-log10(FDR)` = -log10(FDR))
ggplot(filtTab, aes(x = logFC, y=`-log10(FDR)`)) + 
    geom_point(aes(colour=FDR < 0.05), size=1)

Challenge 3 {.challenge} - In Supplementary Materials

Use the txMm to retrieve the exon coordinates for the genes: + ENSMUSG00000021604 + ENSMUSG00000022146 + ENSMUSG00000040118