Updated 2021/04/26
The Book of R: chapter 4.3 (factors), 8 (files)
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
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
| All same data type | Mixed data types | |
|---|---|---|
| 1d | Vector | List |
| 2d | Matrix | Data frame |
| nd | Array |
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"
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
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
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
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
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
ASCII table from https://simple.wikipedia.org/wiki/ASCII
Numbers 32 to 126 define printable characters e.g. decimal 32 = space
Note that there are two different ways to say “next line”
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
#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
write.table(patientRecords, "mydata.txt", sep="\t")
mydata = read.table("mydata.txt")
mydata
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?
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
if (test_expression) {
statement
}
#Note: indenting is optional, but is good style because indented code is easier to read
x = 5
if (x > 0) {
print("Positive number")
}
## [1] "Positive number"
x = -5
if (x > 0) {
print("Positive number")
}
if (test_expression) {
statement1
} else {
statement2
}
x = 5
if (x > 0) {
print("Positive number")
} else {
print("Negative number")
}
## [1] "Positive number"
if (test_expression1) {
statement1
} else if (test_expression2) {
statement2
} else {
statement3
}
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"
< 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