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 show to save data from R. After this exercise you will know how to:
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.
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."
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)
which(nchar(shelley.filt2) == max(nchar(shelley.filt2)))
shelley.filt2[301]
## [1] 304
## [1] "even"
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 with 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.8554244
One 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.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 the write.table
function.
ed <- read.table("https://raw.githubusercontent.com/NBISweden/workshop-r/master/data/lab_loadingdata/example.data", sep=":")
head(ed)
str(ed)
## V1 V2 V3 V4
## 1 bZIP29_1 bZIP29_2 bZIP29_3 WT_1
## 2 <NA> 14.1671426761817 12.9727029171751 14.8181737869517
## 3 1.41723379939617 1.20145379585993 <NA> <NA>
## 4 41.8629060744716 41.1998530830302 42.6965659118674 37.528635786519
## 5 102.422396502516 98.3189689612045 92.0684061403395 104.210417827802
## 6 <NA> 4.00484598619978 3.58953430232514 3.70454344673793
## V5 V6
## 1 WT_2 WT_3
## 2 <NA> <NA>
## 3 1.096828754225 <NA>
## 4 34.8492408728761 48.6456870599298
## 5 97.4183357161658 86.065446336799
## 6 2.64236018063295 6.27033804098888
## 'data.frame': 18946 obs. of 6 variables:
## $ V1: chr "bZIP29_1" NA "1.41723379939617" "41.8629060744716" ...
## $ V2: chr "bZIP29_2" "14.1671426761817" "1.20145379585993" "41.1998530830302" ...
## $ V3: chr "bZIP29_3" "12.9727029171751" NA "42.6965659118674" ...
## $ V4: chr "WT_1" "14.8181737869517" NA "37.528635786519" ...
## $ V5: chr "WT_2" NA "1.096828754225" "34.8492408728761" ...
## $ V6: chr "WT_3" NA NA "48.6456870599298" ...
Encode all NA values as “missing”, at export.
write.table(x=ed, na="missing", file="example_mis.data")
df.test <- read.table("example_mis.data", na.strings="missing")
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).
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;" ...
levels(d.gtf$Chromosome)
## NULL
feature == 'exon'
)
d.gtf.exons <- d.gtf[(d.gtf$Feature == 'exon'),]
aggregate(d.gtf.exons$Feature, by=list(d.gtf.exons$Chromosome), summary)
## 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
d.gtf.gene <- d.gtf[d.gtf$Feature == "gene",]
mean(abs(d.gtf.gene$Start - d.gtf.gene$End))
## [1] 5753.282
sum(d.gtf.gene$Strand == "+") / length(d.gtf.gene$Strand)
## [1] 0.5016231
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
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"