Updated 2021/04/26

The Book of R: chapter 4.3 (factors), 8 (files)


Module 3 review

  • Variable types: list, array, data frame
  • List
    • x = listVariable = list(numericVector, stringVector, logicVector, 77)
    • listVariable[1] #access first list element
    • listVariable[[1]] #access content of first element
  • Array
    • arrayVariable = array(c(vector1,vector2),dim = c(3,3,2))
    • arrayVariable[,,1] #get matrix 1
  • Data frame
    • patientRecords = data.frame(name, age, smoker)
    • patientRecords[1,] #row 1
    • patientRecords[,1] #column 1
    • patientRecords[1,2] #element at row 1, column 2
    • patientRecords$age

Exercise: Querying data frames (10 minutes)

Using the same variable you read “gene_condition_source_id.txt” into for the previous exercise, find how many unique diseases are linked to TP53?

Hint: you can access vector elements matching a specific criteria using: vectorVariable$colName[vectorVariable$colName=="string"]

Hint: ?unique

Exercise solution (5 minutes)

db <- read.table("gene_condition_source_id.txt",
                 header=T,sep="\t",quote="",comment.char="",
                 colClasses="character")
sum(db$AssociatedGenes == "TP53")
## [1] 15
table(db$AssociatedGenes == "TP53")
## 
## FALSE  TRUE 
## 11046    15
table(db$SourceName)
## 
##                                       GeneReviews Human Phenotype Ontology 
##                      616                       40                      623 
##                    MONDO            NCBI curation                     OMIM 
##                     4917                     3815                      128 
##   OMIM phenotypic series                 Orphanet 
##                      619                      303
db$DiseaseName[db$AssociatedGenes == "TP53"]
##  [1] "Basal cell carcinoma, susceptibility to, 7"
##  [2] "Carcinoma of pancreas"                     
##  [3] "Familial cancer of breast"                 
##  [4] "Li-Fraumeni syndrome"                      
##  [5] "Li-Fraumeni syndrome"                      
##  [6] "Choroid plexus papilloma"                  
##  [7] "Choroid plexus papilloma"                  
##  [8] "Osteosarcoma"                              
##  [9] "Bone marrow failure syndrome 5"            
## [10] "Nasopharyngeal carcinoma"                  
## [11] "Li-Fraumeni syndrome 1"                    
## [12] "Hepatocellular carcinoma"                  
## [13] "Colorectal cancer"                         
## [14] "Adrenocortical carcinoma, hereditary"      
## [15] "Glioma susceptibility 1"
db$DiseaseName[which(db$AssociatedGenes == "TP53")]
##  [1] "Basal cell carcinoma, susceptibility to, 7"
##  [2] "Carcinoma of pancreas"                     
##  [3] "Familial cancer of breast"                 
##  [4] "Li-Fraumeni syndrome"                      
##  [5] "Li-Fraumeni syndrome"                      
##  [6] "Choroid plexus papilloma"                  
##  [7] "Choroid plexus papilloma"                  
##  [8] "Osteosarcoma"                              
##  [9] "Bone marrow failure syndrome 5"            
## [10] "Nasopharyngeal carcinoma"                  
## [11] "Li-Fraumeni syndrome 1"                    
## [12] "Hepatocellular carcinoma"                  
## [13] "Colorectal cancer"                         
## [14] "Adrenocortical carcinoma, hereditary"      
## [15] "Glioma susceptibility 1"
unique(db$DiseaseName[db$AssociatedGenes == "TP53"])
##  [1] "Basal cell carcinoma, susceptibility to, 7"
##  [2] "Carcinoma of pancreas"                     
##  [3] "Familial cancer of breast"                 
##  [4] "Li-Fraumeni syndrome"                      
##  [5] "Choroid plexus papilloma"                  
##  [6] "Osteosarcoma"                              
##  [7] "Bone marrow failure syndrome 5"            
##  [8] "Nasopharyngeal carcinoma"                  
##  [9] "Li-Fraumeni syndrome 1"                    
## [10] "Hepatocellular carcinoma"                  
## [11] "Colorectal cancer"                         
## [12] "Adrenocortical carcinoma, hereditary"      
## [13] "Glioma susceptibility 1"
these_have_tp53 <- db$AssociatedGenes == "TP53"
unique(db$DiseaseName[these_have_tp53])
##  [1] "Basal cell carcinoma, susceptibility to, 7"
##  [2] "Carcinoma of pancreas"                     
##  [3] "Familial cancer of breast"                 
##  [4] "Li-Fraumeni syndrome"                      
##  [5] "Choroid plexus papilloma"                  
##  [6] "Osteosarcoma"                              
##  [7] "Bone marrow failure syndrome 5"            
##  [8] "Nasopharyngeal carcinoma"                  
##  [9] "Li-Fraumeni syndrome 1"                    
## [10] "Hepatocellular carcinoma"                  
## [11] "Colorectal cancer"                         
## [12] "Adrenocortical carcinoma, hereditary"      
## [13] "Glioma susceptibility 1"
wc -l gene_condition_source_id.txt
## 11062 gene_condition_source_id.txt

Review: data types vs. data structures

  • Key programming concept
  • Primitive types
    • numeric, logical, string
  • Types
    • vector, matrix, array, list, data frame
  • Data structure
    • A structure created using types
    • e.g. vector, data frame, tree
All same data type Mixed data types
1d Vector List
2d Matrix Data frame
nd Array

One more variable type: Factors

  • Stores categorical variables
  • A variable which can be one of a limited number of values
  • Create a factor using the factor() function
x = c(1,2,3,3,2,1,2,3,1)
x
## [1] 1 2 3 3 2 1 2 3 1
f = factor(x)
f
## [1] 1 2 3 3 2 1 2 3 1
## Levels: 1 2 3
#A factor has levels: these are the possible values
levels(f)
## [1] "1" "2" "3"

Factors and data frames

name = c("patient1", "patient2", "patient3")
age = c(46, 49, 50)
smoker = c(T, F, T)

#By default data.frame() creates factors for string variable columns
patientRecords = data.frame(name, age, smoker)
patientRecords
str(patientRecords)   #Note: "Factor w/ 3 levels"
## 'data.frame':    3 obs. of  3 variables:
##  $ name  : chr  "patient1" "patient2" "patient3"
##  $ age   : num  46 49 50
##  $ smoker: logi  TRUE FALSE TRUE
#Tell data.frame() not to create factors
patientRecords2 = data.frame(name, age, smoker, stringsAsFactors=F)
str(patientRecords2)
## 'data.frame':    3 obs. of  3 variables:
##  $ name  : chr  "patient1" "patient2" "patient3"
##  $ age   : num  46 49 50
##  $ smoker: logi  TRUE FALSE TRUE

Exercise: Working with Factors (10 minutes)

Convert the following protein sequence (just the amino acid sequence part) to a factor and use it to count how many amino acids of each type are in a protein

Hint: ?summary

sp|P04637|P53_HUMAN Cellular tumor antigen p53 OS=Homo sapiens MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

Exercise solution (5 minutes)

prot_string <- "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD"
nchar(prot_string)
## [1] 393
# prot_string <- sub("\\n","",prot_string)
prot <- strsplit(prot_string,split="")
length(prot)
## [1] 1
str(prot)
## List of 1
##  $ : chr [1:393] "M" "E" "E" "P" ...
prot <- prot[[1]]
str(prot)
##  chr [1:393] "M" "E" "E" "P" "Q" "S" "D" "P" "S" "V" "E" "P" "P" "L" "S" ...
testing <- strsplit(c("this","class","rocks"),"")
str(testing)
## List of 3
##  $ : chr [1:4] "t" "h" "i" "s"
##  $ : chr [1:5] "c" "l" "a" "s" ...
##  $ : chr [1:5] "r" "o" "c" "k" ...
table(prot)
## prot
##  A  C  D  E  F  G  H  I  K  L  M  N  P  Q  R  S  T  V  W  Y 
## 24 10 20 30 11 23 12  8 20 32 12 14 45 15 26 38 22 18  4  9
summary(prot)
##    Length     Class      Mode 
##       393 character character
prot <- as.factor(prot)
str(prot)
##  Factor w/ 20 levels "A","C","D","E",..: 11 4 4 13 14 16 3 13 16 18 ...
table(prot)
## prot
##  A  C  D  E  F  G  H  I  K  L  M  N  P  Q  R  S  T  V  W  Y 
## 24 10 20 30 11 23 12  8 20 32 12 14 45 15 26 38 22 18  4  9
summary(prot)
##  A  C  D  E  F  G  H  I  K  L  M  N  P  Q  R  S  T  V  W  Y 
## 24 10 20 30 11 23 12  8 20 32 12 14 45 15 26 38 22 18  4  9

Files

Key Programming Concept

Reminder: Redirect output and input

Redirect output to a file

$ grep cancer genes2phenotype.txt > myfile.txt

Will print all lines in “genes2phenotype.txt” that contain “cancer” to the file “myfile.txt” (no output to screen). Use “cat” or “more” to view the contents of myfile.txt

What is a file

A file is a set of data stored somewhere and given a name Every file is a series of binary numbers, usually thought of as a series of bytes *A bit is the smallest unit of data a computer can store: either 0 or 1 (a switch is off or on)

Binary (base 2) Number (base 10)
00 0
01 1
10 2
11 3

A byte is 8 bits and can represent up to 2^8 = 256 values As a more compact representation, computers often display binary values in hexadecimal (base 16). Two hexademical digits = 16^2 = 256 values = 8 binary digits

Hexademical (base 16) Number (base 10)
0 0
1 1
2 2
… …
a 10
… …
f 15

So a file can be thought of as a series of pairs of hexadecimal digits UNIX commands are available to view a file as this raw data e.g.

xxd fileOutput.png | head
## 00000000: 8950 4e47 0d0a 1a0a 0000 000d 4948 4452  .PNG........IHDR
## 00000010: 0000 0307 0000 01da 0806 0000 0055 6735  .............Ug5
## 00000020: 9700 0000 0173 5247 4200 aece 1ce9 0000  .....sRGB.......
## 00000030: 0009 7048 5973 0000 1712 0000 1712 0167  ..pHYs.........g
## 00000040: 9fd2 5200 0040 0049 4441 5478 01ec dd07  ..R..@.IDATx....
## 00000050: b874 5579 f6f1 ad22 a002 4aef f0d2 11e9  .tUy..."..J.....
## 00000060: bd57 c150 c4cf 8a8a 251a a351 5193 8f8b  .W.P....%..QQ...
## 00000070: 688c 628d 5882 31d1 98cf 2488 05a2 4004  h.b.X.1...$...@.
## 00000080: 312a 5de9 1d01 69d2 5fa4 f75e 6ddf fb5b  1*]...i._..^m..[
## 00000090: ba4e e61d e69c 33a7 cf99 b9d7 75ed b3e7  .N....3.....u...
  • This is how a PNG image file is stored on disk (the first 80 bytes): #{bash} 00000000: 8950 4e47 0d0a 1a0a 0000 000d 4948 4452 .PNG........IHDR 00000010: 0000 0307 0000 01da 0806 0000 0055 6735 .............Ug5 00000020: 9700 0000 0173 5247 4200 aece 1ce9 0000 .....sRGB....... 00000030: 0009 7048 5973 0000 1712 0000 1712 0167 ..pHYs.........g 00000040: 9fd2 5200 0040 0049 4441 5478 01ec dd07 ..R..@.IDATx.... #

  • The data in a file can be represented in different ways

  • A common way for text to be represented is in ASCII format

  • Numbers 32 to 126 define printable characters e.g. decimal 32 = space

  • Note that there are two different ways to say “next line”

    • number 10 (hex: 0A) = line feed
    • number 13 (hex: 0D) = carriage return
  • Different operating systems unfortunately use different ways (encodings) to represent “next line”.

  • macOS and UNIX use line feed (hex: 0A). Windows uses a combination of carriage return and line feed (hex: 0D 0A).

  • See more details at https://en.wikipedia.org/wiki/Newline

  • Exercise on your own: Use the UNIX xxd command to view a text file and see how the text (and a line break) is represented as ASCII codes.

  • Why did we learn this? There are different types of hidden or unexpected characters possible in a file and you sometimes need to be aware of them when reading and writing files. If in doubt, examine your file in a text editor, or use head on the command line to view a file and make sure it is the format you are expecting.

  • Tip: on UNIX, you can use the file command to see what type of format a file is.

file mydata.txt
file fileOutput.png
## mydata.txt: ASCII text, with CRLF line terminators
## fileOutput.png: PNG image data, 775 x 474, 8-bit/color RGBA, non-interlaced
  • A common type of file in data science is a table (like a spreadsheet worksheet, or data frame), usually represented as a tab-delimited format. This type of file should be explored in a spreadsheet application before you try to read it in R

Working with files in R

#Remember to set the working directory
getwd()
## [1] "D:/Dropbox/TA_MMG1002H_2020/p4b"
#setwd("path")
#setwd("C:/MyDoc")

list.files()
##  [1] "_site.yml"                    "aaoutput.txt"                
##  [3] "ASCII-Table-wide.svg"         "assignment2_sequence.txt"    
##  [5] "column.txt"                   "docs"                        
##  [7] "fileOutput.png"               "gene_condition_source_id.txt"
##  [9] "GSE13924_RAW.tar"             "GSM350771.CEL.gz"            
## [11] "GSM350772.CEL.gz"             "index.html"                  
## [13] "index.md"                     "mod_me.html"                 
## [15] "mod_me.Rmd"                   "mod_me_files"                
## [17] "mod1.html"                    "mod1.Rmd"                    
## [19] "mod2.html"                    "mod2.Rmd"                    
## [21] "mod2.utf8.md"                 "mod2_files"                  
## [23] "mod3.html"                    "mod3.Rmd"                    
## [25] "mod4.Rmd"                     "mod5.Rmd"                    
## [27] "mod5exercise.fa"              "mod6.Rmd"                    
## [29] "mod7.knit.md"                 "mod7.Rmd"                    
## [31] "mod8.Rmd"                     "mod9.Rmd"                    
## [33] "mod9.utf8.md"                 "mydata.txt"                  
## [35] "p4b.Rproj"                    "README.md"                   
## [37] "recipebookTOC.png"            "row.txt"                     
## [39] "RStudioOptions.png"           "site_libs"                   
## [41] "sum.R"                        "test.txt"                    
## [43] "thisisapdf.pdf"               "thisisnotapdf.png"           
## [45] "yeast_properties.txt"
dir()  #this does the same thing
##  [1] "_site.yml"                    "aaoutput.txt"                
##  [3] "ASCII-Table-wide.svg"         "assignment2_sequence.txt"    
##  [5] "column.txt"                   "docs"                        
##  [7] "fileOutput.png"               "gene_condition_source_id.txt"
##  [9] "GSE13924_RAW.tar"             "GSM350771.CEL.gz"            
## [11] "GSM350772.CEL.gz"             "index.html"                  
## [13] "index.md"                     "mod_me.html"                 
## [15] "mod_me.Rmd"                   "mod_me_files"                
## [17] "mod1.html"                    "mod1.Rmd"                    
## [19] "mod2.html"                    "mod2.Rmd"                    
## [21] "mod2.utf8.md"                 "mod2_files"                  
## [23] "mod3.html"                    "mod3.Rmd"                    
## [25] "mod4.Rmd"                     "mod5.Rmd"                    
## [27] "mod5exercise.fa"              "mod6.Rmd"                    
## [29] "mod7.knit.md"                 "mod7.Rmd"                    
## [31] "mod8.Rmd"                     "mod9.Rmd"                    
## [33] "mod9.utf8.md"                 "mydata.txt"                  
## [35] "p4b.Rproj"                    "README.md"                   
## [37] "recipebookTOC.png"            "row.txt"                     
## [39] "RStudioOptions.png"           "site_libs"                   
## [41] "sum.R"                        "test.txt"                    
## [43] "thisisapdf.pdf"               "thisisnotapdf.png"           
## [45] "yeast_properties.txt"
#file.create(), file.rename() and many more

Data import and export

  • We usually work with files using read and write functions
write.table(patientRecords, "mydata.txt", sep="\t")

mydata = read.table("mydata.txt")
mydata

Exercise: Practice reading and writing files (10 minutes)

Read the “gene_condition_source_id.txt” file dataFrame = read.table(“gene_condition_source_id.txt”, header=TRUE, sep=“, fill=TRUE, comment.char =”“, quote =”")

count the number of columns and rows

write one column to a file “column.txt” and write one row to a file “row.txt”

examine original and resulting files in a spreadsheet - did you get the correct answer?

Exercise solution (5 minutes)

db <- read.table("gene_condition_source_id.txt",
                 header=T,sep="\t",quote="",comment.char="",
                 colClasses="character")
nrow(db)
## [1] 11061
ncol(db)
## [1] 9
identical(db[,"AssociatedGenes"],db$AssociatedGenes)
## [1] TRUE
write.table(db[,2],file="column.txt",quote=F,row.names=F,col.names=F)

write.table(db[1,],file="row.txt",sep="\t",
            quote=F,row.names=F,col.names=F)
head column.txt
wc -l column.txt
## A2M
## A2M
## A2ML1
## A4GALT
## AAAS
## AAGAB
## AARS1
## AARS1
## AARS2
## AARS2
## 11061 column.txt
less row.txt
awk -F '\t' 'END {print NF}' gene_condition_source_id.txt
awk -F '\t' '{print NF}' row.txt
## 2    A2M     C0002395    Alzheimer disease   MONDO   MONDO:0004975   104300  Feb 19 2020
## 9
## 9

Conditional statements

Key Programming Concept

Conditional statements

  • if condition is true, then run these commands
  • The “if statement”
if (test_expression) {
  
  
  statement
}

#Note: indenting is optional, but is good style because indented code is easier to read
  • Example
x = 5
if (x > 0) {
  print("Positive number")
}
## [1] "Positive number"
x = -5
if (x > 0) {
  print("Positive number")
}

if..else statement

if (test_expression) {
  statement1 
} else {
  statement2
}
  • Example
x = 5
if (x > 0) {
  print("Positive number")
} else {
  print("Negative number")
}
## [1] "Positive number"

if..else if statement

if (test_expression1) {
    statement1
} else if (test_expression2) {
    statement2
} else {
    statement3
}
  • Example
x = 5
if (x > 0) {
  print("Positive number")
} else if (x==0) {
  print("Zero")
} else {
  print("Negative number")
}
## [1] "Positive number"
x = 0
if (x > 0) {
  print("Positive number")
} else if (x==0) {
  print("Zero")
} else {
  print("Negative number")
}
## [1] "Zero"
x = -4
if (x > 0) {
  print("Positive number")
} else if (x==0) {
  print("Zero")
} else {
  print("Negative number")
}
## [1] "Negative number"

Comparison operators

<           less than
<=          less than or equal to
>           greater than
>=          greater than or equal to
==          exactly equal to
!=          not equal to
!x          not x
x | y       x OR y
x & y       x AND y e.g. if((x != 3) & (x > 0))

Note: identical() function is more generic than == and it can work on any R object