Loading data into R
R Programming Foundation for Life Scientists
1 Introduction
Up until now we have mostly created the object we worked with on the fly from within R. The most common use-case is however to read in different data sets that are stored as files, either somewhere on a server or locally on your computer. In this exercise we will test some common ways to import data in R and also how to save data from R. After this exercise you will know how to:
- Read data from txt files and save the information as a vector, data frame or a list.
- Identify missing data and correctly encode this at import
- Check that imported objects are imported correctly
- Read data from online resource
- Write data to a file
2 scan()
The function scan() can be used both to read data from files and directly from keyboard. The function is very flexible and have many different settings that allow to read data in different formats. To read and store a set of words that you type on your keyboard try the following code that will prompt your for input. After each word press enter and R will prompt you for new input. After the last word have been typed press enter twice to get back to your R prompt and have your character vector named words available in R your session.
words <- scan(what=character())We will read in this book chapter. Read the manual for scan and read the text file named book_chapter.txt into R, first as vector and then as a list, with each word in the chapter saved as a entry in the vector or as a single vector in a list.
shelley.vec <- scan(file="https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/book_chapter.txt", what=character())
str(shelley.vec)
shelley.list <- scan(file="https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/book_chapter.txt", what=list(character()))
class(shelley.list)## chr [1:420] "My" "present" "situation" "was" "one" "in" "which" "all" ...
## [1] "list"Check that your newly created objects contain the correct information and have been saved as you have intended eg. each entry of the vector or the list should contain a single word. Once your convinced that you have a sound word vector and list.
- Identify the longest word in your vector.
sort(nchar(shelley.vec), decreasing=TRUE)
which(nchar(shelley.vec) == max(nchar(shelley.vec)))
shelley.vec[381]## [1] 690 12 11 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
## [19] 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
## [37] 9 9 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
## [55] 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 7 7 7
## [73] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6
## [91] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## [109] 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
## [127] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [145] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [163] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [181] 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 4 4
## [199] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [217] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [235] 4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3 3
## [253] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [271] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [289] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [307] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [325] 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [343] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [361] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [379] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [397] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [415] 1 1 1 1 1 1
## [1] 381
## [1] "By the sacred earth on which I kneel, by the shades that wander near me, by the deep and eternal grief that I feel, I swear; and by thee, O Night, and the spirits that preside over thee, to pursue the daemon who caused this misery, until he or I shall perish in mortal conflict. For this purpose I will preserve my life; to execute this dear revenge will I again behold the sun and tread the green herbage of earth, which otherwise should vanish from my eyes forever. And I call on you, spirits of the dead, and on you, wandering ministers of vengeance, to aid and conduct me in my work. Let the cursed and hellish monster drink deep of agony; let him feel the despair that now torments me."- Go back and fix the way you read in the text to make sure that you get a vector with all words in chapter as individual entries also filter any non-letter characters and now identify the longest word.
shelley.vec2 <- scan(file="https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/book_chapter.txt", what=character(), sep=' ', quote=NULL)
shelley.filt2 <- gsub(pattern='[^[:alnum:] ]', replacement="", x=shelley.vec2)
longest <- which(nchar(shelley.filt2) == max(nchar(shelley.filt2)))
shelley.filt2[longest]## [1] "uninterested"3 read.table()
This is the by far most common way to get data into R. As the function creates a data frame at import it will only work for data set that fits those criteria, meaning that the data needs to have a set of columns of equal length that are separated with a common string eg. tab, comma, semicolon etc.
In this code block we first import the data from normalized.txt and accept the defaults for all other arguments in the function. With this settings R will read it as a tab delimited file and will use the first row of the data as colnames (header) and the first column as rownames.
expr.At <- read.table("https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/normalized.txt")
head(expr.At)## bZIP29_1 bZIP29_2 bZIP29_3 WT_1 WT_2 WT_3
## AT1G01020.1 13.572739 14.167143 12.972703 14.8181738 15.904017 11.3270623
## AT1G01030.1 1.417234 1.201454 1.385434 0.8590246 1.096829 0.8596431
## AT1G01040.2 41.862906 41.199853 42.696566 37.5286358 34.849241 48.6456871
## AT1G01050.1 102.422397 98.318969 92.068406 104.2104178 97.418336 86.0654463
## AT1G01060.1 3.216031 4.004846 3.589534 3.7045434 2.642360 6.2703380
## AT1G01070.1 8.230858 17.871625 14.924906 7.9996663 8.824486 13.8554244One does however not have to have all data as a file an the local disk, instead one can read data from online resources. The following command will read in a file from a web server.
url <- 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data'
abalone <- read.table(url, header=FALSE , sep=',')
head(abalone)## V1 V2 V3 V4 V5 V6 V7 V8 V9
## 1 M 0.455 0.365 0.095 0.5140 0.2245 0.1010 0.150 15
## 2 M 0.350 0.265 0.090 0.2255 0.0995 0.0485 0.070 7
## 3 F 0.530 0.420 0.135 0.6770 0.2565 0.1415 0.210 9
## 4 M 0.440 0.365 0.125 0.5160 0.2155 0.1140 0.155 10
## 5 I 0.330 0.255 0.080 0.2050 0.0895 0.0395 0.055 7
## 6 I 0.425 0.300 0.095 0.3515 0.1410 0.0775 0.120 8- Read this example data to R using the
read.table()function. This files consist of gene expression values. Once you have the object in R validate that it looks okay and export it using thewrite.tablefunction.
ed <- read.table("https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/example.data", sep=":", header = T)
head(ed)
str(ed)## bZIP29_1 bZIP29_2 bZIP29_3 WT_1 WT_2 WT_3
## 1 NA 14.167143 12.972703 14.818174 NA NA
## 2 1.417234 1.201454 NA NA 1.096829 NA
## 3 41.862906 41.199853 42.696566 37.528636 34.849241 48.645687
## 4 102.422397 98.318969 92.068406 104.210418 97.418336 86.065446
## 5 NA 4.004846 3.589534 3.704543 2.642360 6.270338
## 6 8.230858 17.871625 14.924906 7.999666 NA NA
## 'data.frame': 18945 obs. of 6 variables:
## $ bZIP29_1: num NA 1.42 41.86 102.42 NA ...
## $ bZIP29_2: num 14.2 1.2 41.2 98.3 4 ...
## $ bZIP29_3: num 12.97 NA 42.7 92.07 3.59 ...
## $ WT_1 : num 14.8 NA 37.5 104.2 3.7 ...
## $ WT_2 : num NA 1.1 34.85 97.42 2.64 ...
## $ WT_3 : num NA NA 48.65 86.07 6.27 ...Encode all NA values as “missing”, at export.
write.table(x=ed, na="missing", file="example_write.txt")- Read in the file you just created and double-check that you have the same data as earlier.
df.test <- read.table("example_write.txt", na.strings="missing")- Analysing genome annotation in R using read.table
For this exercise we will load a GTF file into R and calculate some basic summary statistics from the file. In the first part we will use basic manipulations of data frames to extract the information. In the second part you get a try out a library designed to work with annotation data, that stores the information in a more complex format, that allow for easy manipulation and calculation of summaries from genome annotation files.
For those not familiar with the GTF format it is a file format containing annotation information for a genome. It does not contain the actual DNA sequence of the organism, but instead refers to positions along the genome.
A valid GTF file should contain the following tab delimited fields (taken from the ensembl home page).
- seqname - name of the chromosome or scaffold; chromosome names can be given with or without the ‘chr’ prefix.
- source - name of the program that generated this feature, or the data source (database or project name)
- feature - feature type name, e.g. gene, transcript, exon, CDS, start_codon, end_codon
- start - Start position of the feature, with sequence numbering starting at 1.
- end - End position of the feature, with sequence numbering starting at 1.
- score - A floating point value.
- strand - defined as + (forward) or - (reverse).
- frame - One of ‘0’, ‘1’ or ‘2’. ‘0’ indicates that the first base of the feature is the first base of a codon, ‘1’ that the second base is the first base of a codon, and so on..
- attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|
| 1 | transcribed_unprocessed_pseudogene | gene | 11869 | 14409 | . | + | . | gene_id; “ENSG00000223972”; |
| 1 | processed_transcript | transcript | 11869 | 14409 | . | + | . | gene_id; “ENSG00000223972”; |
The last column can contain a large number of attributes that are semicolon-separated.
As these files for many organisms are large we will in this exercise use the latest version of Drosophila melanogaster genome annotation available at ftp://ftp.ensembl.org/pub/release-86/gtf/drosophila_melanogaster that is small enough for analysis even on a laptop.
Download the file named Drosophila_melanogaster.BDGP6.86.gtf.gz to your computer. Unzip this file and keep track of where your store the file.
With this done read this file into R using the function read.table() and add meaningful column names to the table.
d.gtf <- read.table("Drosophila_melanogaster.BDGP6.86.gtf", header=FALSE, comment.char="#", sep="\t")
colnames(d.gtf) <- c("Chromosome","Source","Feature", "Start","End","Score","Strand","Frame","Attribute")Prior to any analysis you should make sure that your attempt to read in the file has worked as expected. This can for example be done by having a look at the dimension of the stored object and making sure that it has the structure you expect.
dim(d.gtf)
str(d.gtf)## [1] 538684 9
## 'data.frame': 538684 obs. of 9 variables:
## $ Chromosome: chr "3R" "3R" "3R" "3R" ...
## $ Source : chr "FlyBase" "FlyBase" "FlyBase" "FlyBase" ...
## $ Feature : chr "gene" "transcript" "exon" "gene" ...
## $ Start : int 722370 722370 722370 835381 835381 835381 835381 835381 869486 869486 ...
## $ End : int 722621 722621 722621 2503907 2503907 835491 835491 835383 869548 869548 ...
## $ Score : chr "." "." "." "." ...
## $ Strand : chr "-" "-" "-" "+" ...
## $ Frame : chr "." "." "." "." ...
## $ Attribute : chr "gene_id FBgn0085804; gene_name CR41571; gene_source FlyBase; gene_biotype pseudogene;" "gene_id FBgn0085804; transcript_id FBtr0114258; gene_name CR41571; gene_source FlyBase; gene_biotype pseudogene"| __truncated__ "gene_id FBgn0085804; transcript_id FBtr0114258; exon_number 1; gene_name CR41571; gene_source FlyBase; gene_bio"| __truncated__ "gene_id FBgn0267431; gene_name CG45784; gene_source FlyBase; gene_biotype protein_coding;" ...- How many chromosome names can be found in the annotation file?
length(levels(as.factor(d.gtf$Chromosome)))## [1] 57- How many exons is there in total and per chromosome? (hint: first extract lines that have
feature == 'exon')
d.gtf.exons <- d.gtf[(d.gtf$Feature == 'exon'),]
nrow(d.gtf.exons)
aggregate(d.gtf.exons$Feature, by=list(d.gtf.exons$Chromosome), summary)## [1] 186414
## Group.1 x.Length x.Class x.Mode
## 1 211000022278158 1 character character
## 2 211000022278279 1 character character
## 3 211000022278282 1 character character
## 4 211000022278298 1 character character
## 5 211000022278307 1 character character
## 6 211000022278309 1 character character
## 7 211000022278436 1 character character
## 8 211000022278449 2 character character
## 9 211000022278498 1 character character
## 10 211000022278522 1 character character
## 11 211000022278603 1 character character
## 12 211000022278604 1 character character
## 13 211000022278664 1 character character
## 14 211000022278724 1 character character
## 15 211000022278750 1 character character
## 16 211000022278760 2 character character
## 17 211000022278875 1 character character
## 18 211000022278877 1 character character
## 19 211000022278878 1 character character
## 20 211000022278879 1 character character
## 21 211000022278880 1 character character
## 22 211000022278985 1 character character
## 23 211000022279055 1 character character
## 24 211000022279108 1 character character
## 25 211000022279132 1 character character
## 26 211000022279134 1 character character
## 27 211000022279165 1 character character
## 28 211000022279188 3 character character
## 29 211000022279222 1 character character
## 30 211000022279264 1 character character
## 31 211000022279342 1 character character
## 32 211000022279392 1 character character
## 33 211000022279446 1 character character
## 34 211000022279528 1 character character
## 35 211000022279529 1 character character
## 36 211000022279531 1 character character
## 37 211000022279555 1 character character
## 38 211000022279681 1 character character
## 39 211000022279708 1 character character
## 40 211000022280133 1 character character
## 41 211000022280328 4 character character
## 42 211000022280341 1 character character
## 43 211000022280347 1 character character
## 44 211000022280481 1 character character
## 45 211000022280494 2 character character
## 46 211000022280645 1 character character
## 47 211000022280703 1 character character
## 48 2L 32747 character character
## 49 2R 38551 character character
## 50 3L 34347 character character
## 51 3R 43159 character character
## 52 4 3165 character character
## 53 dmel_mitochondrion_genome 37 character character
## 54 rDNA 21 character character
## 55 Unmapped_Scaffold_8 14 character character
## 56 X 34136 character character
## 57 Y 182 character character
by(data=d.gtf$Feature, d.gtf[,"Chromosome"], summary)## d.gtf[, "Chromosome"]: 211000022278158
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278279
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278282
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278298
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278307
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278309
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278436
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278449
## Length Class Mode
## 4 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278498
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278522
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278603
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278604
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278664
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278724
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278750
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278760
## Length Class Mode
## 7 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278875
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278877
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278878
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278879
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278880
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022278985
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279055
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279108
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279132
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279134
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279165
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279188
## Length Class Mode
## 12 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279222
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279264
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279342
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279392
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279446
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279528
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279529
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279531
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279555
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279681
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022279708
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280133
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280328
## Length Class Mode
## 14 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280341
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280347
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280481
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280494
## Length Class Mode
## 6 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280645
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 211000022280703
## Length Class Mode
## 3 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 2L
## Length Class Mode
## 96740 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 2R
## Length Class Mode
## 110108 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 3L
## Length Class Mode
## 100745 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 3R
## Length Class Mode
## 125133 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: 4
## Length Class Mode
## 7844 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: dmel_mitochondrion_genome
## Length Class Mode
## 146 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: rDNA
## Length Class Mode
## 59 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: Unmapped_Scaffold_8
## Length Class Mode
## 43 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: X
## Length Class Mode
## 97193 character character
## ------------------------------------------------------------
## d.gtf[, "Chromosome"]: Y
## Length Class Mode
## 504 character character- Filter the data frame to only retain gene annotations
d.gtf.gene <- d.gtf[d.gtf$Feature == "gene",]- What is the average gene length of in the Drosophila genome?
mean(abs(d.gtf.gene$Start - d.gtf.gene$End))## [1] 5753.282- What fraction of the genes are encoded on the plus strand of the genome.
sum(d.gtf.gene$Strand == "+") / length(d.gtf.gene$Strand)## [1] 0.5016231- What is the median and mean length of the exons found on chromosome 3R in the data set?
d.gtf3R <- d.gtf[d.gtf$Chromosome == "3R",]
exon.position <- d.gtf3R[d.gtf3R$Feature == "exon",c("Start", "End")]
median(abs(exon.position$Start - exon.position$End))
mean(abs(exon.position$Start - exon.position$End))## [1] 251
## [1] 468.6693- Do the same calculations for the chromosomes 2L, 2R, 3L, 4, X and Y using a for loop.
chr <- c("2L","2R","3L","4","X","Y")
for (i in chr) {
d.gtf.tmp <- d.gtf[d.gtf$Chromosome == i,]
exon.position <- d.gtf.tmp[d.gtf.tmp$Feature == "exon", c("Start", "End")]
exon.med <- median(abs(exon.position$Start - exon.position$End))
exon.mean <- mean(abs(exon.position$Start - exon.position$End))
txt <- sprintf("The median and mean exon length for %s is %g and %g, respectively", i, exon.med, exon.mean)
print(txt)
}## [1] "The median and mean exon length for 2L is 279 and 502.617, respectively"
## [1] "The median and mean exon length for 2R is 225 and 437.187, respectively"
## [1] "The median and mean exon length for 3L is 255 and 502.116, respectively"
## [1] "The median and mean exon length for 4 is 198 and 429.573, respectively"
## [1] "The median and mean exon length for X is 257 and 526.301, respectively"
## [1] "The median and mean exon length for Y is 410.5 and 657.055, respectively"