Tenth-R-Practice-exercise-merging-annotation-data-into-a-gene-expression-analysis-results-data-frame.Rmd

---
 output:
  md_document:
    variant: markdown_github
# output: pdf_document
---

## Tenth R Practice exercise merging annotation data into a gene expression analysis results data frame.Rmd

### Alan E. Berger  Feb 17, 2020

### available at https://github.com/AlanBerger/Practice-programming-exercises-for-R 

## Introduction

This is the tenth in a sequence of programming exercises in "composing" an R function 
to carry out a particular task. Several of these "exercise files" likely 
will take several sessions to master the content.  The material below practices composing a logical
sequence of steps to program a function that will accomplish a specified task, and 
preparing a corresponding data frame. 

The idea of this set of exercises is to practice correct use of R constructs and 
built in functions (functions that "come with" the basic R installation), while learning how 
to "put together" a correct sequence of blocks of commands that will obtain the desired result.  
Note these exercises are quite cumulative - one should do them in order. 

In these exercises, there will be a statement of what your function should do 
(what are the input variables and what the function should return) and a sequence of "hints". 
To get the most out of these exercises, try to write your function using as few hints as possible.  
Note there are often several ways to write a function that will obtain the correct result. 
For these exercises the directions and hints may point toward a particular approach intended to 
practice particular constructs in R and a particular line of reasoning, 
even if there is a more efficent way to obtain the same result. 
There may also be an existing R function or package that will do what is stated for a given 
practice exercise, but here the point is to practice formulating a logical sequence of steps, 
with each step a section of code, to obtain a working function, not to find an existing 
solution or a quick solution using a more powerful R construct that is better addressed later on.

## Motivation for this exercise

In some cases, such as with a gene expression data set, one will want to combine analysis results as 
obtained in the previous exercise with annotation information on the probes and on the genes that is in a 
separate file that can be read in as a data frame.

In the R code below we repeat the analysis, done in the previous exercise, of a small subset of gene expression 
data comparing expression levels in PBMC samples from patients with Wegener's granulomatosis (WG) with samples 
from normal controls (NC). We then also read in a small subset of the annotation file for the Illumina microarray 
platform used to measure these expression levels. The web site containing the full expression data set is: 
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE18885 
and the web site containing the full annotation data for the microarray platform used in obtaining this data is: 
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6104  

What we will want to do is, conceptually, for each row **r** of the analysis results data frame, "find" the row **ra** 
of the annotation data frame that has the same Illumina Probe_ID and in effect append selected columns of row ra 
from the annotation data frame to row r of the analysis results data frame. R has a function **merge** that will 
do this, but for the first exercise we will practice using basic R constructs to compose R code that will do this - 
the second exercise here will use the merge function. First read through the code below, that provides the data frames 
one will use.

``` {r}
##############  analyze the gene expression data

# the url for reading the little gene expression data file into an R data frame using 
# read.delim  (for reading in tab delimited text files) is given in the next 3 lines 
url.for.data.file <- "https://raw.githubusercontent.com/AlanBerger/
Practice-programming-exercises-for-R/master/tiny-subset-of-GSE18885-
gene-expression-data-9-genes-WG-5-samples-Normal-Control-4-samples.tab.txt" 
 
# read in the data as a data frame
ma <- read.delim(url.for.data.file, nrows = 9, check.names = FALSE, 
                 stringsAsFactors = FALSE)

# display ma
ma

# now, in a for loop, get the p-values and fold changes 
num.genes <- nrow(ma)  # the number of genes in this data frame

gene <- ma$gene  # the column of gene names
probe.vec <- ma[[1]]  # the column of Illumina Probe_IDs
# get vectors to hold the p-value and fold change values
p.value <- numeric(num.genes)
fold.change <- numeric(num.genes)

for (i in 1:num.genes) {
# get the vector for the WG expression values and the vector 
# for the NC expression values for the ith gene
   NCvec <- unlist(ma[i, 3:6])
   WGvec <- unlist(ma[i, 7:11])

# calculate the p-value and fold change
   pval <- t.test(NCvec, WGvec)$p.value  # two-sided unequal variance (Welch) t-test
   p.value[i] <- pval 
   WG.over.NC.fold.change <- 2^(mean(WGvec) - mean(NCvec))
   fold.change[i] <- WG.over.NC.fold.change  
}

# Construct the desired data frame. 
analysis.results <- data.frame(probe.vec, gene, p.value, fold.change, 
                               stringsAsFactors = FALSE, check.names = FALSE)
colnames(analysis.results) <- c("Illumina PROBE_ID", "gene", "two-sided p-value", 
                                "WG/NC fold change")
analysis.results 

##############  read in the annotation data file (a small subset of the full annotation)

# read in the short edited Illumina microarray annotation data file called
# GPL6104-Illumina-microarray-platform-annotation-from-GEO-repository-
# small-subset-edited-example-Feb12.tab.txt

url.for.annotation.file <- 
"https://raw.githubusercontent.com/AlanBerger/Practice-programming-exercises-for-R/
master/GPL6104-Illumina-microarray-platform-annotation-from-GEO-repository-small-
subset-edited-example-Feb12.tab.txt"

annotation.df  <- read.delim(url.for.annotation.file, nrows = 15, check.names = FALSE, 
                 stringsAsFactors = FALSE)

# Note the use of nrows = 15 since there is information on the 
# source of this data in later rows that should not be read in as data.
# The choice check.names = FALSE "tells" R to leave the column headers as is 

annotation.df

# We see that this annotation data file has data for more Illumina probes than 
# are in the analysis results data frame, and that the probe IDs are not in the 
# same order as in the analysis results data frame.  For the purpose of the 
# practice exercise below we will only append the columns containing the gene name,  
# Chromosome number (which chromosome the gene is located on) 
# and the short desciption of the protein encoded by the gene. 
# Repeating the gene name gives a indicator to use to double check that the "merge"
# was done correctly.

# Note this data frame has an example of more than 1 probe for a given 
# gene (ATP2B1), where different parts of the same gene are "queried". 
 
# To keep things simpler with this example, for each probe ID in the analysis results, 
# there is 1 row of the annotation data frame with the same probe ID. 
 
# The merge function can handle the case where there is no matching probe ID in 
# the annotation file for a probe ID in the analysis results data frame, in which case 
# we would want to append NA's indicating that information is not available in the 
# annotation file being used.

columns.to.keep <- c(1, 2, 4, 6)  # keep just these columns of the annotation data frame 
# to have print outs easy to see
# we need to keep the probe IDs in column 1 to be able to match rows

annotation.df <- annotation.df[, columns.to.keep]
# from now on annotation.df will refer to this shortened version of the
# annotation data frame

annotation.df 

```

## Programming Exercise:  Append to analysis.results information from the annotation data frame

Approach: Form a vector of row numbers, call it annot.rows, such that for each row r of the analysis.results 
data frame; annot.rows[r] will contain the row of the annotation.df data frame that has the same 
Illumina Probe_ID as does row r of analysis.results  
Then column binding  annotation.df[annot.rows, ] to analysis.results (using **cbind**) will yield the desired result.

Hint: Use a for loop, and use the **which** function to find, for each row r of analysis.results the
row number ra of annotation.df that has the same probe ID as does row r of analysis.results 

A working version of R code which does this is given below.

``` {r}
# use the analysis.results and annotation.df obtained in the R code above.
nrows <- nrow(analysis.results) 
# create the integer vector annot.rows of length nrows to hold the 
# row numbers of annotation.df matching (with respect to the probe ID) the 
# analysis.results rows 
annot.rows <- vector(mode = "integer", length = nrows)  

# get the Illumina probe IDs vector from the annotation data frame
annotation.df.probe.ids <- annotation.df[[1]]  

for (r in 1:nrows) {
   probe.id <- analysis.results[r, 1]
#  find the row ra of annotation.df whose Illumina probe ID matches probe.id
ra <- which(annotation.df.probe.ids == probe.id)
if (length(ra) != 1) stop("did not find unique matching probe id row")
annot.rows[r] <- ra
}

# append the correct rows (correctly lined up) of annotation.df to analysis.results 
analysis.results.with.annotation <- cbind(analysis.results, annotation.df[annot.rows, ])
analysis.results.with.annotation  # display it

# Note the row numbers are from the rows of annotation.df whose probe IDs 
# matched up with the those in analysis.results
```    

## Second exercise: use the R **merge** function to append matching annotation lines to analysis.results 

The R merge function can combine two data frames in various ways. See for example the web page by Joachim Schork 
which is a page in https://statisticsglobe.com/ titled "Merge Data Frames by Column Names in R (3 Examples)": 
https://statisticsglobe.com/r-merging-data-frames-by-column-names-merge-function   
See also the R help on the merge function (via ? merge). 

While it is good practice to use basic R constructs until they are easy for you to use, using an available 
R function can greatly simplify code which then makes it easier to keep free of bugs. Code that uses the 
merge function is given below. The merge function is capable of a number of types of merging in addition to the 
example below.

``` {r}
# Use the R merge function to append annotation to the analysis results

# recall that annotation.df is referring to the shortened version of the annotation
# merged.df <- merge(x = analysis.results, y = annotation.df, 
#              by.x = "Illumina PROBE_ID", by.y = "Illumina Probe_ID", 
#              all.x = TRUE, all.y = FALSE, sort = FALSE)

# What the above call to merge will do (once the comment symbols # are removed) is: 
# use the by.x = "Illumina PROBE_ID" column of analysis.results as the "guide" 
# and for each row r of analysis.results, the merge function will in effect search 
# to find the row ra of annotation.df such that the entry of row ra in column 
# by.y = "Illumina Probe_ID" of annotation.df  matches the entry of row r in 
# the column by.x = "Illumina PROBE_ID" of analysis.results
# Note these 2 column names are not exactly the same so we need to specify
# the column names in x and y to be used to do matching of rows, 
# using the arguments by.x and by.y
# The merge function will, in effect, append row ra of annotation.df to row r 
# of analysis.results 
# The choice all.x = TRUE means: if there is no match for the entry of row r in 
# the column "Illumina PROBE_ID" of analysis.results anywhere in the column 
# "Illumina Probe_ID" of annotation.df, then a row of NA's is appended to row r 
# of analysis.results
# The choice all.y = FALSE means don't include rows of annotation.df other than 
# those appended to analysis.results as desribed above. 
# The choice sort = FALSE means do not sort the resulting data frame  
# (any sorting would have been done for this call to merge using 
# the "Illumina PROBE_ID" column). 

merged.df <- merge(x = analysis.results, y = annotation.df, 
              by.x = "Illumina PROBE_ID", by.y = "Illumina Probe_ID", 
              all.x = TRUE, all.y = FALSE, sort = FALSE)
# display it
merged.df 

# Note the Illumina Probe_ID column of annotation.df is NOT included in merge.df
# Let's check that merged.df is the same as analysis.results.with.annotation obtained 
# above. First we need to remove the Illumina Probe_ID column from annotation.df 
# that is included in analysis.results.with.annotation before we check.
analysis.results.with.annotation <- analysis.results.with.annotation[, -5]

# check if they are the same
identical(analysis.results.with.annotation, merged.df)


# What happened? they looked the same -- so now a little "adventure" 
# in finding out what happened -- this sort of thing "comes with the territory" 
# when programming in any language (they each have their own quirks).

# Let's look closer
attributes(analysis.results.with.annotation)

attributes(merged.df)

# So the row names were different

# Looks like we can fix this by setting the row names of
# analysis.results.with.annotation to be those for merged.df
row.names(analysis.results.with.annotation) <- row.names(merged.df)
identical(analysis.results.with.annotation, merged.df)


# Now what ????   Let's look at the attributes again

attributes(analysis.results.with.annotation)

attributes(merged.df)

# So the row names for analysis.results.with.annotation are 1:9 as characters
# and the row names for merged.df are 1:9 as integers -
# As I said, every language has its quirks
row.names(analysis.results.with.annotation) <- 1:9
attributes(analysis.results.with.annotation)

# Now if they aren't identical we really do have problems
identical(analysis.results.with.annotation, merged.df) 

# So some semblance of order is restored.  The problem was 
# row.names(merged.df) returned a character vector
str(row.names(merged.df))

# One final verification: I'm going to remove the row of the annotation
# data frame corresponding the probe ID for the BPI gene 
# and then use the merge function
annotation.df <- annotation.df[-5, ]

merged.df <- merge(x = analysis.results, y = annotation.df, 
              by.x = "Illumina PROBE_ID", by.y = "Illumina Probe_ID", 
              all.x = TRUE, all.y = FALSE, sort = FALSE)
# display it
merged.df 

# Note merge filled in NA's for the annotation columns for the row for the probe ID 
# corresponding to BPI as expected.
# The merge function also placed the row for which there was no match for the 
# probe ID  in the annotation file at the bottom of the merged data frame.

# This illustrates the kind of "exploring" one should do when using a new R function,
# particularly one that has a somewhat complex range of options and for which the 
# output has a range of possibilities, in order to be confident about what it will 
# do when called a certain way
 
```

Hope this was informative and good practice. 
The next exercise will contain further practice in using data frames, and point out some 
types of logical mistakes that may result in actual output that, however, is incorrect, 
rather than an error message. This is the most dangerous type of mistake, in that if the 
incorrect output is not obviously wrong, the mistake might not be recognized until it 
causes serious consequences.  That is why it is always wise to do, whenever possible, test 
runs for cases where one knows or can independently calculate the true result. 

= = = = = = = = = = = = = = = = = = = = = = = = 

This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. 
To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter 
to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. There is a full version of this license at this web site: 
https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode   

Note the reader should not infer any endorsement or recommendation or approval for the material in this article from 
any of the sources or persons cited above or any other entities mentioned in this article.