{"id":2239,"date":"2018-04-02T15:07:24","date_gmt":"2018-04-02T14:07:24","guid":{"rendered":"http:\/\/pcool.dyndns.org:8080\/statsbook\/?page_id=2239"},"modified":"2025-06-29T09:35:19","modified_gmt":"2025-06-29T08:35:19","slug":"ladys-slipper-orchid","status":"publish","type":"page","link":"https:\/\/pcool.dyndns.org\/index.php\/ladys-slipper-orchid\/","title":{"rendered":"Lady's Slipper Orchid"},"content":{"rendered":"\n

Please make sure the necessary packages (seqinr1<\/a><\/sup> and Biostrings2<\/a><\/sup>) are installed as described<\/a> to allow analysis. Furthermore, the ggplot23<\/a><\/sup> and reshape24<\/a><\/sup> packages should be loaded.<\/p>\n\n\n\n

The genetic code for the lady’s slipper orchid (<\/a>Cypripedium calceolus)<\/i> can be downloaded<\/a>. Select all the text (Ctrl-A), copy (Ctrl-C) and paste (Ctrl-V) it into a text editor. Save the file as ls_orchid.fasta<\/strong> (without<\/strong><\/em> a txt extension!) in the working directory (folder) of R.<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

Load the required package, import the fasta file as a list (orchid) and review the structure (str):<\/p>\n\n\n\n

library(seqinr)\norchid <- read.fasta(file='~\/Desktop\/ls_orchid.fasta')\nstr(orchid)<\/mark>\nList of 94\n $ gi|2765658|emb|Z78533.1|CIZ78533: 'SeqFastadna' chr [1:740] \"c\" \"g\" \"t\" \"a\" ...\n  ..- attr(*, \"name\")= chr \"gi|2765658|emb|Z78533.1|CIZ78533\"\n  ..- attr(*, \"Annot\")= chr \">gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\"\n.....\n.....\nP.purpuratum 5.8S rRNA gene and ITS1 and ITS2 DNA\"\n $ gi|2765564|emb|Z78439.1|PBZ78439: 'SeqFastadna' chr [1:592] \"c\" \"a\" \"t\" \"t\" ...\n  ..- attr(*, \"name\")= chr \"gi|2765564|emb|Z78439.1|PBZ78439\"\n  ..- attr(*, \"Annot\")= chr \">gi|2765564|emb|Z78439.1|PBZ78439 P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA\"<\/mark><\/em><\/code><\/pre>\n\n\n\n
The information is stored in a list. Every sequence can be viewed separately. To view the first sequence:<\/p>\n\n\n\n
orchid[1]<\/mark>\n$`gi|2765658|emb|Z78533.1|CIZ78533`\n  [1] \"c\" \"g\" \"t\" \"a\" \"a\" \"c\" \"a\" \"a\" \"g\" \"g\" \"t\" \"t\" \"t\" \"c\" \"c\" \"g\" \"t\" \"a\" \"g\" \"g\" \"t\" \"g\" \"a\" \"a\" \"c\" \"c\" \"t\" \"g\" \"c\" \"g\" \"g\" \"a\" \"a\"\n.....\n.....\n[694] \"t\" \"t\" \"c\" \"g\" \"g\" \"a\" \"t\" \"g\" \"t\" \"g\" \"a\" \"c\" \"c\" \"c\" \"c\" \"a\" \"g\" \"g\" \"t\" \"c\" \"a\" \"g\" \"g\" \"c\" \"g\" \"g\" \"g\" \"g\" \"g\" \"c\" \"a\" \"c\" \"c\"\n[727] \"c\" \"g\" \"c\" \"t\" \"g\" \"a\" \"g\" \"t\" \"t\" \"t\" \"a\" \"c\" \"g\" \"c\"\nattr(,\"name\")\n[1] \"gi|2765658|emb|Z78533.1|CIZ78533\"\nattr(,\"Annot\")\n[1] \">gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA\"\nattr(,\"class\")\n[1] \"SeqFastadna\"<\/mark><\/em><\/code><\/pre>\n\n\n\n
and the last:<\/p>\n\n\n\n
orchid[94]<\/mark>\n$`gi|2765564|emb|Z78439.1|PBZ78439`\n  [1] \"c\" \"a\" \"t\" \"t\" \"g\" \"t\" \"t\" \"g\" \"a\" \"g\" \"a\" \"t\" \"c\" \"a\" \"c\" \"a\" \"t\" \"a\" \"a\" \"t\" \"a\" \"a\" \"t\" \"t\" \"g\" \"a\" \"t\" \"c\" \"g\" \"a\" \"g\" \"t\" \"t\"\n.....<\/mark><\/em>\n.....\n[529] \"g\" \"t\" \"g\" \"g\" \"a\" \"c\" \"g\" \"a\" \"a\" \"c\" \"t\" \"a\" \"t\" \"g\" \"c\" \"t\" \"a\" \"c\" \"a\" \"a\" \"c\" \"a\" \"a\" \"a\" \"a\" \"t\" \"t\" \"g\" \"t\" \"t\" \"g\" \"t\" \"g\"\n[562] \"c\" \"a\" \"a\" \"a\" \"t\" \"g\" \"c\" \"c\" \"c\" \"c\" \"g\" \"g\" \"t\" \"t\" \"g\" \"g\" \"c\" \"c\" \"g\" \"t\" \"t\" \"t\" \"a\" \"g\" \"t\" \"t\" \"g\" \"g\" \"g\" \"c\" \"c\"\nattr(,\"name\")\n[1] \"gi|2765564|emb|Z78439.1|PBZ78439\"\nattr(,\"Annot\")\n[1] \">gi|2765564|emb|Z78439.1|PBZ78439 P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA\"\nattr(,\"class\")\n[1] \"SeqFastadna\"<\/mark><\/em><\/code><\/pre>\n\n\n\n
To show the total number of sequences:<\/p>\n\n\n\n
> number_sequences <- length(orchid)\n> number_sequences\n<\/mark>[1] 94<\/mark><\/em><\/code><\/pre>\n\n\n\n
To show the length of the first and the last sequence:<\/p>\n\n\n\n
length(orchid[[1]])<\/em><\/span>\n[1] 740<\/em><\/span>\nlength(orchid[[number_sequences]])<\/em><\/span>\n[1] 592<\/em><\/span><\/code><\/pre>\n\n\n\n
To create a table of the nucleotides and calculate the GC (Guanine – Cytosine) content (more stable as three hydrogen bonds):<\/p>\n\n\n\n
table <- table(orchid[[1]])<\/span><\/em>\ntable<\/span><\/em>\na c g t <\/em><\/span>\n144 200 241 155<\/span> <\/span><\/em>\nstr(table)<\/span><\/em>\n 'table' int [1:4(1d)] 144 200 241 155<\/em><\/span>\n - attr(*, \"dimnames\")=List of 1<\/em><\/span>\n ..$ : chr [1:4] \"a\" \"c\" \"g\" \"t\"<\/em><\/span>\ntable[1]<\/span><\/em>\n a <\/em><\/span>\n144<\/span> <\/span><\/em>\ntable[2]<\/span><\/em>\n c <\/em><\/span>\n200<\/span> <\/span><\/em>\ntable[3]<\/span><\/em>\n g <\/em><\/span>\n241<\/span> <\/span><\/em>\ntable[4]<\/span><\/em>\n t <\/em><\/span>\n155<\/span> <\/span><\/em>\n# calculate the gc content:<\/span><\/em>\ngc_content <- (table[3] + table[2]) \/ sum(table)<\/span><\/em>\ngc_content<\/span><\/em>\n g <\/em><\/span>\n0.5959459 <\/em><\/span><\/code><\/pre>\n\n\n\n
Or, the easy way:<\/p>\n\n\n\n
GC(orchid[[1]])<\/em><\/span>\n[1] 0.5959459<\/em><\/span><\/code><\/pre>\n\n\n\n
So, almost 60%. It is also easy to count single nucleotides or combinations thereof:<\/p>\n\n\n\n
count(orchid[[1]], 1)<\/em><\/span>\na c g t <\/em><\/span>\n144 200 241 155<\/span> <\/em>\ncount(orchid[[1]], 2)<\/em><\/span>\naa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt <\/em><\/span>\n41 27 39 37 39 66 61 33 53 70 73 45 11 36 68 40<\/span> <\/em>\ncount(orchid[[1]], 3)<\/em><\/span>\naaa aac aag aat aca acc acg act aga agc agg agt ata atc atg att caa cac cag cat cca ccc ccg cct <\/em><\/span>\n 7 11 12 11 3 15 7 2 6 11 16 6 4 10 15 8 11 4 12 12 14 23 20 9 <\/em><\/span>\ncga cgc cgg cgt cta ctc ctg ctt gaa gac gag gat gca gcc gcg gct gga ggc ggg ggt gta gtc gtg gtt <\/em><\/span>\n 9 14 25 13 3 11 9 10 19 9 12 13 16 18 22 13 19 26 15 13 3 10 22 10 <\/em><\/span>\ntaa tac tag tat tca tcc tcg tct tga tgc tgg tgt tta ttc ttg ttt <\/em><\/span>\n 4 3 3 1 6 10 11 9 19 19 17 13 1 5 22 12<\/span> <\/em>\n# the output is a table object. So the extract is easy:<\/em><\/span>\ncount(orchid[[1]], 1)[3] <\/em><\/span>\n g <\/em><\/span>\n241 <\/em><\/span><\/code><\/pre>\n\n\n\n
Or better, refer by name:<\/p>\n\n\n\n
count(orchid[[1]], 1)[\"g\"]<\/em><\/span>\n g <\/em><\/span>\n241<\/span> <\/em>\n# similarly to count a specific triplets of nucleotides:<\/em><\/span>\ncount(orchid[[1]], 3)[\"att\"]<\/em><\/span>\natt <\/em><\/span>\n 8<\/em><\/span><\/code><\/pre>\n\n\n\n
Show the GC content for every sequence, reshape the data with reshape25<\/a><\/sup>  and plot with ggplot26<\/a><\/sup>:<\/p>\n\n\n\n
Show sequence length and GC content in each sequence, reshape and plot the data:<\/p>\n\n\n\n
library(reshape2)<\/mark><\/em>\nlibrary(ggplot2)<\/mark><\/em>\nseq_lengths <- lapply(orchid, length)<\/em><\/span>\nseq_lengths_gc <- lapply(orchid, GC)<\/em><\/span>\ngcs_plot <- melt(seq_lengths_gc)<\/span><\/em>\nlengths_plot <- melt(seq_lengths)<\/span><\/em>\nplot_data <- data.frame(length = lengths_plot, gc = gcs_plot)<\/span><\/em>\nggplot(plot_data, aes(x = length.value, y = gc.value)) +<\/span><\/em>\ngeom_point() +<\/span><\/em>\ntheme_bw()<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
So the majority of sequences have a GC content of around 50%. There is however a group of longer sequences that have a GC content of more than 55%.<\/p>\n\n\n\n
Sliding window plots are often used to show how the GC content varies along the sequence. Here, we create a sliding window plot of the longest (38th) sequence per 100 base pairs:<\/p>\n\n\n\n
# sliding window plot to show GC content per 100 base pairs for sequence 38:<\/em><\/span>\nstarts <- seq(1, length(orchid[[38]]) - 100, by = 100)<\/em><\/span>\nn <- length(starts) # how many chunks?<\/em><\/span>\nchunkGCs <- numeric(n) # null vector that has the same length as n<\/em><\/span>\nfor (i in 1:n) {<\/em><\/span>\nchunk <- orchid[[38]][starts[i]:(starts[i] + 99)]<\/em><\/span>\nchunkGC <- GC(chunk)<\/em><\/span>\nchunkGCs[i] <- chunkGC<\/em><\/span>\n}<\/em><\/span>\nchunkGCs<\/em><\/span>\n[1] 0.43 0.60 0.55 0.45 0.59 0.57 0.49<\/em><\/span>\n\nplot_data_2 <- data.frame(location = starts, GC = chunkGCs)<\/em><\/span>\ndev.new()<\/em><\/span>\nggplot(plot_data_2, aes(x = location, y = GC)) +<\/em><\/span>\ngeom_line() +<\/em><\/span>\ntheme_bw()<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
It is straight forward to get the single and triple letter amino acids coded for:<\/p>\n\n\n\n
getTrans(orchid[[38]]) # this will give the standard one letter abbreviation<\/em><\/span>\n [1] \"R\" \"N\" \"K\" \"V\" \"S\" \"V\" \"G\" \"E\" \"P\" \"P\" \"E\" \"G\" \"S\" \"L\" \"L\" \"R\" \"S\" \"H\" \"N\" \"N\" \"*\" \"S\" \"R\"<\/em><\/span>\n- - -<\/em><\/span>\n[254] \"T\" \"P\" \"G\" \"Q\" \"*\" \"A\" \"T\" \"R\" \"*\" \"V\"<\/em><\/span>\naaa(getTrans(orchid[[38]])) # will give the three letter abbreviation<\/em><\/span>\n [1] \"Arg\" \"Asn\" \"Lys\" \"Val\" \"Ser\" \"Val\" \"Gly\" \"Glu\" \"Pro\" \"Pro\" \"Glu\" \"Gly\" \"Ser\" \"Leu\" \"Leu\"<\/em><\/span>\n-  --<\/em><\/span>\n[256] \"Gly\" \"Gln\" \"Stp\" \"Ala\" \"Thr\" \"Arg\" \"Stp\" \"Val\"<\/em><\/span><\/code><\/pre>\n\n\n\n
Dot plots are traditionally used to compare DNA sequences. The seqinr7<\/a><\/sup> package contains a function that creates a dot plot (called dotPlot). The number of nucleotides (wsize), step (wstep for overlap) and number of matches (nmatch) can be set. To compare single nucleotides of the 38th orchid sequence with itself:<\/p>\n\n\n\n
dev.new()<\/em><\/span>\ndotPlot(orchid[[38]], orchid[[38]], wsize = 1, wstep = 1, nmatch = 1)<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
The diagonal line is uninterrupted, indicating the sequences are the same (obviously as the sequence was compared to itself).<\/p>\n\n\n\n
Or to compare triplets of nucleotides of the 38th orchid sequence with itself, specifying that all three nucleotides need to match (nmatch =3):<\/p>\n\n\n\n
dev.new()<\/em><\/span>\ndotPlot(orchid[[38]], orchid[[38]], wsize = 3, wstep = 3, nmatch = 3)<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
The dot plot function included in the package is not very versatile. It would be better to plot with ggplot28<\/a><\/sup>. The main advantage of open source software is that the code of the functions is available to the user. The code of the dot plot was changed so that it returns a matrix (with sequence 1 as rows and sequence 2 as columns). The row names are the first sequence and the column names the second sequence. Copy and paste this function<\/a> (called dot_plot_data) into the R console.<\/p>\n\n\n\n
The new function (dot_plot_data) should now be available to R (it can also be viewed and downloaded<\/a> from the download page). To access the function comparing three nucleotides and matching all three nucleotides takes the same format as the dotPlot function from the seqinr package9<\/a><\/sup> .<\/p>\n\n\n\n
new_plot <- dot_plot_data(orchid[[38]], orchid[[38]], wsize = 3, wstep = 3, nmatch = 3)<\/em><\/span>\nModification of dotPlot function in package seqinr.<\/span><\/em>\nThe function returns a matrix for plotting.<\/span><\/em>\nRow names contain Sequence1, Col names Sequence2.<\/span><\/em>\nUse the reshape2 package to convert and plot with ggplot2 and geom_raster.<\/span><\/em><\/code><\/pre>\n\n\n\n
To show the structure and class of the returned object:<\/p>\n\n\n\n
str(new_plot)<\/em><\/span>\n logi [1:263, 1:263] TRUE FALSE FALSE FALSE FALSE FALSE ...<\/em><\/span>\n - attr(*, \"dimnames\")=List of 2<\/em><\/span>\n ..$ : chr [1:263] \"cgt\" \"aac\" \"aag\" \"gtt\" ...<\/em><\/span>\n ..$ : chr [1:263] \"cgt\" \"aac\" \"aag\" \"gtt\" ...<\/em><\/span>\nclass(new_plot)<\/em><\/span>\n[[1] \"matrix\" \"array\"<\/em><\/span><\/code><\/pre>\n\n\n\n
The data is now available as a matrix called new_plot. The row names are sequence 1 and the column names sequence 2. Unfortunately, ggplot210<\/a><\/sup> can’t plot matrices directly. Consequently, the matrix need to be converted to a data frame suitable for plotting,<\/p>\n\n\n\n
First, extract the row names and column names as vectors:<\/p>\n\n\n\n
sequence_x <- rownames(new_plot)<\/span><\/em>\nsequence_y <- colnames(new_plot)<\/span><\/em><\/code><\/pre>\n\n\n\n
<\/p>\n\n\n\n
Secondly, remove the row and column names from the matrix (the names are not unique and otherwise ggplot211<\/a><\/sup> will group by them)<\/p>\n\n\n\n
rownames(new_plot) <- NULL<\/span><\/em>\ncolnames(new_plot) <- NULL<\/span><\/em><\/code><\/pre>\n\n\n\n
Thirdly, use the reshape2 package12<\/a><\/sup> to reshape the data to allow plotting (the row and column numbers are unique, representing the number in the sequence):<\/p>\n\n\n\n
new_plot_2 <- melt(new_plot)<\/span><\/em><\/code><\/pre>\n\n\n\n
The data frame new_plot_2 has 3 columns (variables in tidy format). Var1<\/strong> represents sequence 1 (x coordinate) and Var2<\/strong> sequence 2 (y coordinate). The value <\/strong>variable contains the value stored in the matrix. To plot the data:<\/p>\n\n\n\n
dot_plot <-\u00a0 ggplot(new_plot_2, <\/em><\/span>\n     aes(x = Var1, y = Var2, fill = value)) +<\/em><\/span>\ngeom_raster() +<\/em><\/span>\nscale_fill_manual(values = c(\"lightgray\", \"blue\"), name = \"Pair\", labels = c(\"Different\", \"Same\")) +<\/em><\/span>\nscale_x_continuous(\"Sequence 1\") +<\/em><\/span>\nscale_y_continuous(\"Sequence 2\") + <\/em><\/span>\ntheme_bw()<\/em><\/span> \ndot_plot<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
It is possible to show the nucleotides on the axes. However, they will become too crowded. It is possible to zoom in though:<\/p>\n\n\n\n
dot_plot_2 <- ggplot(new_plot_2, <\/em><\/span>\n     aes(x = Var1, y = Var2, fill = value)) +<\/em><\/span>\ngeom_raster() +<\/em><\/span>\nscale_fill_manual(values = c(\"lightgray\", \"blue\"), name = \"Pair\", labels = c(\"Different\", \"Same\")) +<\/em><\/span>\nscale_x_continuous(\"Sequence 1\", breaks = 1: length(sequence_x), labels = sequence_x, limit = c(10,20)) +<\/em><\/span>\nscale_y_continuous(\"Sequence 2\", breaks = 1: length(sequence_y), labels = sequence_y, limit = c(10,20)) + <\/em><\/span>\ntheme_bw() <\/em><\/span>\ndot_plot_2<\/em><\/span>\nWarning message:<\/em><\/span>\nRemoved 69048 rows containing missing values (geom_raster).<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Finally, to save fasta files is straight forward. To save the orchid sequences to a fasta file (orchid_2) in the working directory:<\/p>\n\n\n\n
write.fasta(names<\/span> = c(1<\/span>:94<\/span>), sequences<\/span> = orchid<\/span>, file.out<\/span> = \"orchid_2.fasta\"<\/span>)<\/span><\/em><\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"
Please make sure the necessary packages (seqinr and Biostrings) are installed as described to allow analysis. Furthermore, the ggplot2 and reshape2 packages should be loaded. The genetic code for the lady’s slipper orchid (Cypripedium calceolus) can be downloaded. Select all the text (Ctrl-A), copy (Ctrl-C) and paste (Ctrl-V) it into a text editor. Save the file as ls_orchid.fasta (without a […]<\/p>\n","protected":false},"author":2,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"inline_featured_image":false,"footnotes":""},"class_list":["post-2239","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2239","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/comments?post=2239"}],"version-history":[{"count":3,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2239\/revisions"}],"predecessor-version":[{"id":4552,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2239\/revisions\/4552"}],"wp:attachment":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/media?parent=2239"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}