{"id":2232,"date":"2018-04-02T15:04:08","date_gmt":"2018-04-02T14:04:08","guid":{"rendered":"http:\/\/pcool.dyndns.org:8080\/statsbook\/?page_id=2232"},"modified":"2025-06-29T17:30:29","modified_gmt":"2025-06-29T16:30:29","slug":"idh2-gene","status":"publish","type":"page","link":"https:\/\/pcool.dyndns.org\/index.php\/idh2-gene\/","title":{"rendered":"IDH2 gene"},"content":{"rendered":"\n

Please make sure the necessary packages (seqinr1<\/a><\/sup> and Biostrings2<\/a><\/sup>) are installed<\/a> as described to allow analysis. <\/p>\n\n\n\n

The Biostrings<\/a> and pwalign<\/a> packages are part of Bioconductor<\/a> and installation is a little different:<\/p>\n\n\n\n

install.packages('BiocManager')\nBiocManager::install('Biostrings')<\/mark><\/em>\nBiocManager::install('pwalign') # this package is also needed<\/mark>\nlibrary(Biostrings)<\/mark><\/em><\/code><\/pre>\n\n\n\n
If you used reshape2 or dplyr, you may need to restart R if you get a can’t unload package error message.<\/p>\n\n\n\n
Furthermore, the ggplot23<\/a><\/sup> and reshape24<\/a><\/sup> libraries should be loaded.<\/p>\n\n\n\n
In this section, the genetic sequence of the isocitrate dehydrogenase (IDH) 2 gene in humans is compared to that of orangutans. Mutations in the gene have been associated with Ollier’s disease and Maffucci syndrome in humans.<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
IDH is an enzyme that catalyses the oxidative decarboxylation of isocitrate. In humans, IDH has three forms: IDH1, IDH2 and IDH3. IDH3 catalyses the reaction with nicotinamide adenine dinucleotide (NAD+) as cofactor within the citric acid cycle, whilst IDH1 and IDH2 catalyse the reaction outside the citric acid cycle (with nicotinamide adenine dinucleotide phosphate (NADP+) as cofactor). Consequently, IDH1 and IDH2 are NADP+ dependent and IDH3 is NAD+ dependent.<\/p>\n\n\n\n
IDH2 is a mitochondrial enzyme,<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
encoded on chromosome 15 (in humans and orangutans):<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
The genetic code for human (homo sapiens<\/em>) IDH2 can be downloaded here<\/a> and the genetic code for an orangutan (pongo abelii<\/em>) IDH2 here<\/a>. Select all the text (Ctrl-A), copy (Ctrl-C) and paste (Ctrl-V) it into a text editor. Save the files as idh2_human.fasta <\/strong>and idh2_orangutan.fasta<\/strong> respectively (without<\/strong><\/em> a txt extension) in the working directory (folder) of R.<\/p>\n\n\n\n
Load the files into R and examen the structure (str):<\/p>\n\n\n\n
library(seqinr)<\/mark><\/em>\nhuman <- read.fasta(file = \"idh2_human.fasta\")<\/em><\/span>\norangutan <- read.fasta(file = \"idh2_orangutan.fasta\")<\/em><\/span>\nstr(human)<\/em><\/span>\nList of 1<\/em><\/span>\n $ NG_023302.1:4923-23499:Class 'SeqFastadna' atomic [1:18577] t c c c ...<\/em><\/span>\n .. ..- attr(*, \"name\")= chr \"NG_023302.1:4923-23499\"<\/em><\/span>\n .. ..- attr(*, \"Annot\")= chr \">NG_023302.1:4923-23499 Homo sapiens isocitrate dehydrogenase (NADP(+)) 2, mitochondrial (IDH2), RefSeqGene (LR\"| __truncated__<\/em><\/span>\nstr(orangutan)<\/em><\/span>\nList of 1<\/em><\/span>\n $ NC_036918.1:c70426964-70408539:Class 'SeqFastadna' atomic [1:18426] g c a a ...<\/em><\/span>\n .. ..- attr(*, \"name\")= chr \"NC_036918.1:c70426964-70408539\"<\/em><\/span>\n .. ..- attr(*, \"Annot\")= chr \">NC_036918.1:c70426964-70408539 Pongo abelii isolate Susie chromosome 15, Susie_PABv2, whole genome shotgun sequence\"<\/em><\/span><\/code><\/pre>\n\n\n\n
The sequences can be viewed by:<\/p>\n\n\n\n
human[1]<\/mark>\n$`NG_023302.1:4923-23499`\n    [1] \"t\" \"c\" \"c\" \"c\" \"c\" \"g\" \"g\" \"c\" \"a\" \"a\" \"g\" \"g\" \"c\" \"c\" \"c\" \"a\" \"a\" \"t\" \"g\" \"g\" \"g\" \"g\" \"c\" \"g\" \"g\" \"c\" \"a\" \"g\" \"g\" \"c\" \"c\" \"c\"\n.....\n.....\n[18561] \"a\" \"a\" \"a\" \"a\" \"g\" \"c\" \"t\" \"c\" \"t\" \"t\" \"c\" \"a\" \"c\" \"a\" \"a\" \"a\" \"a\"\nattr(,\"name\")\n[1] \"NG_023302.1:4923-23499\"\nattr(,\"Annot\")\n[1] \">NG_023302.1:4923-23499 Homo sapiens isocitrate dehydrogenase (NADP(+)) 2, mitochondrial (IDH2), RefSeqGene (LRG_611) on chromosome 15\"\nattr(,\"class\")\n[1] \"SeqFastadna\"\n\n<\/mark>orangutan[1]<\/mark>\n$`NC_036918.1:c70426964-70408539`\n    [1] \"g\" \"c\" \"a\" \"a\" \"g\" \"g\" \"c\" \"c\" \"c\" \"a\" \"a\" \"t\" \"g\" \"g\" \"g\" \"g\" \"c\" \"g\" \"g\" \"c\" \"g\" \"g\" \"g\" \"c\" \"c\" \"c\" \"g\" \"g\" \"c\" \"a\" \"g\" \"c\"\n.....\n.....\n[18401] \"t\" \"a\" \"g\" \"c\" \"t\" \"a\" \"c\" \"t\" \"a\" \"a\" \"a\" \"a\" \"a\" \"g\" \"c\" \"t\" \"c\" \"t\" \"t\" \"c\" \"a\" \"c\" \"a\" \"a\" \"a\" \"a\"\nattr(,\"name\")\n[1] \"NC_036918.1:c70426964-70408539\"\nattr(,\"Annot\")\n[1] \">NC_036918.1:c70426964-70408539 Pongo abelii isolate Susie chromosome 15, Susie_PABv2, whole genome shotgun sequence\"\nattr(,\"class\")\n[1] \"SeqFastadna\"<\/mark><\/em><\/code><\/pre>\n\n\n\n
The sequences are not of equal length:<\/p>\n\n\n\n
length(human[[1]])<\/em><\/span>\n[1] 18577<\/em><\/span>\nlength(orangutan[[1]])<\/em><\/span>\n[1] 18426<\/em><\/span><\/code><\/pre>\n\n\n\n
G-C bonds are more stable as they contain 3 hydrogen bonds whilst the A-T combination only has two. It is straight forward to create a table and calculate the G-C content:<\/p>\n\n\n\n
human_table <- table(human[[1]])<\/em><\/span>\nhuman_table<\/em><\/span>\na c g t <\/em><\/span>\n4244 4645 5144 4544 <\/em><\/span>\nhuman_table[1]<\/em><\/span>\n a <\/em><\/span>\n4244 <\/em><\/span>\nhuman_table[2]<\/em><\/span>\n c <\/em><\/span>\n4645 <\/em><\/span>\nhuman_table[3]<\/em><\/span>\n g <\/em><\/span>\n5144 <\/em><\/span>\nhuman_table[4]<\/em><\/span>\n t <\/em><\/span>\n4544 <\/em><\/span>\n# calculate the g-c content:<\/em><\/span>\ngc_content_human <- (human_table[3] + human_table[2]) \/ sum(human_table)<\/em><\/span>\ngc_content_human<\/em><\/span>\n g <\/em><\/span>\n0.5269419 <\/em><\/span>\n# or calculate GC content with the build in function:<\/em><\/span>\nGC(human[[1]])<\/em><\/span>\n[1] 0.5269419<\/em><\/span>\n\norangutan_table <- table(orangutan[[1]])<\/em><\/span>\norangutan_table<\/em><\/span>\na c g t <\/em><\/span>\n4181 4591 5152 4502<\/span> <\/em><\/span>\norangutan_table[1]<\/em><\/span>\n a <\/em><\/span>\n4181<\/span> <\/em><\/span>\norangutan_table[2]<\/em><\/span>\n c <\/em><\/span>\n4591<\/span> <\/em><\/span>\norangutan_table[3]<\/em><\/span>\n g <\/em><\/span>\n5152<\/span> <\/em><\/span>\norangutan_table[4]<\/em><\/span>\n t <\/em><\/span>\n4502<\/span> <\/em><\/span>\n# calculate the g-c content:<\/em><\/span>\ngc_content_orangutan <- (orangutan_table[3] + orangutan_table[2]) \/ sum(orangutan_table)<\/em><\/span>\ngc_content_orangutan<\/em><\/span>\n g <\/em><\/span>\n0.5287637<\/span> <\/em><\/span>\n# or calculate GC content with the build in function:<\/em><\/em><\/span>\nGC(orangutan[[1]])<\/em><\/span>\n[1] 0.5287637<\/em><\/span><\/code><\/pre>\n\n\n\n
There is also a function to count individual nucleotides, or combinations of nucleotides:<\/p>\n\n\n\n
count(human[[1]], 1) # for single nucleotides\n<\/mark>a c g t \n4244 4645 5144 4544 <\/mark>\ncount(human[[1]], 2) # for pairs of nucleotides<\/mark><\/em>\naa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt <\/em><\/span>\n1100 834 1479 830 1317 1486 393 1449 1152 1286 1690 1016 675 1039 1582 1248<\/span> <\/em>\ncount(human[[1]], 3) # for triplets of nucleotides (codons)<\/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>\n438 162 295 204 272 249 70 243 341 348 512 278 157 177 257 239 219 291 538 269 437 464 123 462 <\/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 68 117 127 81 161 367 585 336 265 231 452 204 341 440 126 379 403 488 489 310 178 208 393 237 <\/em><\/span>\ntaa tac tag tat tca tcc tcg tct tga tgc tgg tgt tta ttc ttg ttt <\/em><\/span>\n178 150 194 153 267 333 74 365 340 333 562 347 179 286 347 436<\/span> <\/em>\n# the output is a table object. So the extract is easy:<\/em><\/span>\ncount(human[[1]], 1)[3] # using reference by place<\/em><\/span>\n g <\/em><\/span>\n5144<\/span> <\/em>\ncount(human[[1]], 1)[\"g\"] # using reference by name<\/em><\/span>\n g <\/em><\/span>\n5144<\/span> <\/em>\n# similarly:<\/em><\/span>\ncount(human[[1]], 3)[\"att\"]<\/em><\/span>\natt <\/em><\/span>\n239<\/em> <\/span>\n\ncount(orangutan[[1]], 1)\n<\/mark>a c g t\n 4181 4591 5152 4502<\/mark>\ncount(orangutan[[1]], 2)<\/mark><\/em>\naa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt<\/span><\/em>\n 1076 834 1463 807 1284 1444 413 1450 1157 1277 1701 1017 664 1036 1574 1228<\/span><\/em>\ncount(orangutan[[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<\/span><\/em>\n 422 176 280 197 267 244 70 253 344 348 497 274 144 179 246 238 216 278 533 257 413 443 134 454<\/span><\/em>\n cga cgc cgg cgt cta ctc ctg ctt gaa gac gag gat gca gcc gcg gct gga ggc ggg ggt gta gtc gtg gtt<\/span><\/em>\n 67 127 140 79 169 376 591 314 263 230 453 211 339 430 129 379 402 477 509 313 175 206 408 228<\/span><\/em>\n taa tac tag tat tca tcc tcg tct tga tgc tgg tgt tta ttc ttg ttt<\/span><\/em>\n 175 150 197 142 265 327 80 364 344 324 555 351 176 275 329 448<\/span><\/em>\n# the output is a table object. So the extract is easy:<\/em><\/span>\ncount(orangutan[[1]], 1)[3] # reference by place<\/em><\/span>\n g<\/span><\/em>\n 5152<\/span><\/em>\ncount(orangutan[[1]], 1)[\"g\"]  # reference by name<\/em><\/span>\n g<\/span><\/em>\n 5152<\/span><\/em>\n# similarly:<\/em><\/span>\ncount(orangutan[[1]], 3)[\"att\"]<\/em><\/span>\n att<\/span><\/em>\n 238<\/span><\/em><\/code><\/pre>\n\n\n\n
So, the length of the sequences is similar (18577<\/em> in the human sequence and 18426 in the orangutan sequence), the combination att occurs 239 times in humans and 238 times in orangutans and the G-C content is also similar (52.7% in humans and 52.9% in orangutans).<\/p>\n\n\n\n
To show the sequence lengths and G-C content in a bar chart (after reshaping with reshape25<\/a><\/sup>:<\/p>\n\n\n\n
# show sequence length and GC content in each sequence (human and orangutan):<\/em><\/span>\nseq_lengths_human <- lapply(human, length)<\/em><\/span>\nseq_lengths_human_gc <- lapply(human, GC)<\/em><\/span>\nseq_lengths_orangutan <- lapply(orangutan, length)<\/em><\/span>\nseq_lengths_orangutan_gc <- lapply(orangutan, GC)<\/em><\/span>\n\nlengths_human_plot <- reshape2::melt(seq_lengths_human)<\/em><\/span>\ngcs_human_plot <- reshape2::melt(seq_lengths_human_gc)<\/em><\/span>\nlengths_orangutan_plot <- reshape2::melt(seq_lengths_orangutan)<\/em><\/span>\ngcs_orangutan_plot <- reshape2::melt(seq_lengths_orangutan_gc)<\/em><\/span>\n\nplot_data <- data.frame(human_length = lengths_human_plot,human_gc = gcs_human_plot,orangutan_length = lengths_orangutan_plot,orangutan_gc = gcs_orangutan_plot)<\/em><\/span>\nplot_data<\/em><\/span>\n human_length.value human_length.L1 human_gc.value human_gc.L1<\/em><\/span>\n1 18577 NG_023302.1:4923-23499 0.5269419 NG_023302.1:4923-23499<\/em><\/span>\n orangutan_length.value orangutan_length.L1 orangutan_gc.value<\/em><\/span>\n1 18426 NC_036918.1:c70426964-70408539 0.5287637<\/em><\/span>\n orangutan_gc.L1<\/em><\/span>\n1 NC_036918.1:c70426964-70408539<\/em><\/span>\n\ndev.new()<\/em><\/span>\nplot_data %>%<\/em><\/span>\nggplot() +<\/em><\/span>\ngeom_bar(aes(x = human_length.value, y = human_gc.value), stat = \"identity\", colour = \"black\", fill = \"blue\", width = 10) +<\/em><\/span>\ngeom_bar(aes(x = orangutan_length.value, y = orangutan_gc.value), stat = \"identity\", colour = \"black\", fill = \"orange\", width = 10) +<\/em><\/span>\nscale_x_continuous(name = \"Sequence Length\", limits = c(18400, 18600), breaks = c(18400, 18500, 18600, 18700)) +<\/em><\/span>\nscale_y_continuous(name = \"G-C proportion\", limits = c(0.0, 0.6), breaks = c(0.0, 0.2, 0.4, 0.6)) +<\/em><\/span>\ntheme_bw()<\/em><\/span>\n<\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
A sliding window plot is useful to compare the G-C content in different parts of the sequence. Here, we compare the G-C sequence in blocks of 1000 nucleotides and plot them with ggplot26<\/a><\/sup>:<\/p>\n\n\n\n
# sliding window plot to show gc content per 1000 base pairs for sequences:<\/em><\/span>\nhuman_starts <- seq(1, length(human[[1]]) - 1000, by = 1000)<\/em><\/span>\nhuman_n <- length(human_starts) # how many chuncks?<\/em><\/span>\nhuman_chunkGCs <- numeric(human_n) # null vector that has the same length as n<\/em><\/span>\nfor (i in 1:human_n) {<\/em><\/span>\nhuman_chunk <- human[[1]][human_starts[i]:(human_starts[i] + 999)]<\/em><\/span>\nhuman_chunkGC <- GC(human_chunk)<\/em><\/span>\nhuman_chunkGCs[i] <- human_chunkGC<\/em><\/span>\n}<\/em><\/span>\nhuman_chunkGCs<\/em><\/span>\n [1] 0.730 0.521 0.495 0.534 0.470 0.427 0.456 0.452 0.507 0.497 0.526 0.555 0.584 0.525 0.550<\/em><\/span>\n[16] 0.514 0.488 0.633<\/em><\/span>\n# and for orangutans:<\/em><\/span>\norangutan_starts <- seq(1, length(orangutan[[1]]) - 1000, by = 1000)<\/em><\/span>\norangutan_n <- length(orangutan_starts) # how many chuncks?<\/em><\/span>\norangutan_chunkGCs <- numeric(orangutan_n) # null vector that has the same length as n<\/em><\/span>\nfor (i in 1:orangutan_n) {<\/em><\/span>\norangutan_chunk <- orangutan[[1]][orangutan_starts[i]:(orangutan_starts[i] + 999)]<\/em><\/span>\norangutan_chunkGC <- GC(orangutan_chunk)<\/em><\/span>\norangutan_chunkGCs[i] <- orangutan_chunkGC<\/em><\/span>\n}<\/em><\/span>\norangutan_chunkGCs<\/em><\/span>\n [1] 0.719 0.523 0.499 0.528 0.463 0.435 0.464 0.451 0.498 0.510 0.536 0.556 0.585 0.534 0.553<\/em><\/span>\n[16] 0.528 0.494 0.638<\/em><\/span>\n\ndev.new()<\/em><\/span>\nggplot() + <\/em><\/span>\ngeom_line(aes(x = human_starts, y = human_chunkGCs), colour = \"blue\") +<\/em><\/span>\ngeom_line(aes(x = orangutan_starts, y = orangutan_chunkGCs), colour = \"orange\") +<\/em><\/span>\nscale_x_continuous(name = \"Sequence Position\", limits = c(0, 16000), breaks = c(5000, 10000, 15000)) +<\/em><\/span>\nscale_y_continuous(name = \"GC proportion\", limits = c(0.0, 0.6), breaks = c(0.0, 0.2, 0.4, 0.6)) +<\/em><\/span>\ntheme_bw()<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
To show the amino acids the nucleotides code for:<\/p>\n\n\n\n
getTrans(human[[1]]) # this will give the standard one letter abbreviation<\/em><\/span>\n [1] \"S\" \"P\" \"A\" \"R\" \"P\" \"N\" \"G\" \"A\" \"A\" \"G\" \"P\" \"A\" \"A\" \"P\" \"P\" \"R\" \"W\" \"C\" \"P\" \"R\" \"G\" \"Q\" \"R\"<\/em><\/span>\n- - - -<\/em><\/span>\n[6188] \"K\" \"A\" \"L\" \"H\" \"K\"<\/em><\/span>\naaa(getTrans(human[[1]])) # will give the three letter abbreviation<\/em><\/span>\n [1] \"Ser\" \"Pro\" \"Ala\" \"Arg\" \"Pro\" \"Asn\" \"Gly\" \"Ala\" \"Ala\" \"Gly\" \"Pro\" \"Ala\" \"Ala\" \"Pro\" \"Pro\"<\/em><\/span>\n- - - -<\/em><\/span>\n[6181] \"Gly\" \"Ala\" \"Ile\" \"Phe\" \"Ser\" \"Tyr\" \"Stp\" \"Lys\" \"Ala\" \"Leu\" \"His\" \"Lys\"<\/em><\/span>\n\ngetTrans(orangutan[[1]]) # this will give the standard one letter abbreviation<\/em><\/span>\n [1] \"A\" \"R\" \"P\" \"N\" \"G\" \"A\" \"A\" \"G\" \"P\" \"A\" \"A\" \"P\" \"P\" \"R\" \"W\" \"C\" \"P\" \"R\" \"G\" \"P\" \"R\" \"P\" \"P\"<\/em><\/span>\n- - - - <\/em><\/span>\n[6142] \"K\"<\/em><\/span>\naaa(getTrans(orangutan[[1]])) # will give the three letter abbreviation<\/em><\/span>\n [1] \"Ala\" \"Arg\" \"Pro\" \"Asn\" \"Gly\" \"Ala\" \"Ala\" \"Gly\" \"Pro\" \"Ala\" \"Ala\" \"Pro\" \"Pro\" \"Arg\" \"Trp\"<\/em><\/span>\n- - - - - <\/em><\/span>\n[6136] \"Thr\" \"Lys\" \"Lys\" \"Leu\" \"Phe\" \"Thr\" \"Lys\"<\/em><\/span><\/code><\/pre>\n\n\n\n
Traditionally, dot plots are used to compare DNA sequences. The seqinr package7<\/a><\/sup> does include a function to create dot plots. The number of nucleotides (wsize), step (wstep for overlap) and number of matches (nmatch) can be set. To compare chunks of 9 nucleotides, without overlap,  in the DNA sequences of the IDH2 gene in humans and orangutans:<\/p>\n\n\n\n
dotPlot(human[[1]], orangutan[[1]], wsize = 9, wstep = 9,nmatch = 9, col = c(\"white\", \"black\"), xlab =deparse(substitute(human[[1]])), ylab =deparse(substitute(orangutan[[1]])))<\/em><\/span><\/code><\/pre>\n\n\n\n
Depending on the computer, it can take some time for the comparison to be completed:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
The plot produced by the standard plot is a little difficult to assess , but this will be addressed later. If the sequences were the same, there should be a diagonal line from the origin of the plot to the upper right corner (similar as seen in the dot plot here<\/a>).  There seems little resemblance. Why is this? The sequences are not of the same length and not aligned. To compare the sequences, it is necessary to align them first.<\/p>\n\n\n\n
To compare sequences, first define a scoring matrix<\/strong><\/em>. Here a score of +2 is given for a match and a score of -1 for a mismatch:<\/p>\n\n\n\n
sigma <- pwalign::nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = TRUE)<\/em><\/span>\nsigma # Print out the matrix<\/em><\/span>\n   A  C  G  T\nA  2 -1 -1 -1\nC -1  2 -1 -1\nG -1 -1  2 -1\nT -1 -1 -1  2<\/mark><\/em><\/code><\/pre>\n\n\n\n
If you used reshape2 or dplyr, you may need to restart R if you get a can’t unload package error message.<\/p>\n\n\n\n
It is also required to convert the sequences to a string:<\/p>\n\n\n\n
human_string <- getSequence(human[[1]], as.string = TRUE)<\/span><\/em>\norangutan_string <- getSequence(orangutan[[1]], as.string = TRUE) <\/span><\/em>\nhuman_string <- toupper(human_string[[1]]) # convert to upper case<\/span><\/em>\norangutan_string <- toupper(orangutan_string[[1]]) # convert to upper case<\/span><\/em><\/code><\/pre>\n\n\n\n
Now align the strings using the scoring matrix and setting gaps:<\/p>\n\n\n\n
comparison <- pwalign::pairwiseAlignment(human_string,orangutan_string, substitutionMatrix = sigma,gapOpening = -2,gapExtension = -8,scoreOnly = FALSE)\ncomparison<\/mark>\nGlobal PairwiseAlignmentsSingleSubject (1 of 1)\npattern: TCCCCGGCAAGGCCCAATGGGGCGGCAGGCCCGGCAGCCCCGCCCCGGTGGTGCCCGCGCGGCCAGCGCCCGCCAGGCC...ATGTTTTGCATACTGTAATTTATATTGCCCTTGGAACACATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA\nsubject: ------GCAAGGCCCAATGGGGCGGCGGGCCCGGCAGCCCCGCCCCGGTGGTGTCCGCGCGGCCCGCGCCCGCCAGGCC...ATGTTTTGCATACTGTAATTTATATTGCCCTTGGAACACATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA\nscore: 32478 <\/mark><\/em><\/code><\/pre>\n\n\n\n
Considering there are 18577 base pairs and a score of 2 is given for a match, a score of 32478 is very convincing for similarity!<\/p>\n\n\n\n
A function to review the aligned sequences is not part of the standard packages. However, Avril Coghlan8<\/a><\/sup> has published a function for this on her website<\/a>. The function can be can be seen on this page<\/a> and can be downloaded from here<\/a>. To use the function, copy and paste the function into the R console and execute it. The function should now be available to R as “printPairwiseAlignment”. To show the alignment using the function:<\/p>\n\n\n\n
printPairwiseAlignment(alignment=comparison) \nLoading required package: pwalign\n\nAttaching package: \u2018pwalign\u2019\n\nThe following objects are masked from \u2018package:Biostrings\u2019:\n\n    aligned, alignedPattern, alignedSubject, compareStrings, deletion, errorSubstitutionMatrices, indel, insertion, mismatchSummary, mismatchTable,\n    nedit, nindel, nucleotideSubstitutionMatrix, pairwiseAlignment, PairwiseAlignments, PairwiseAlignmentsSingleSubject, pattern, pid,\n    qualitySubstitutionMatrices, stringDist, unaligned, writePairwiseAlignments<\/mark>\n\n[1] \"GCAAGGCCCAATGGGGCGGCAGGCCCGGCAGCCCCGCCCCGGTGGTGCCCGCGCGGCCAG 60\"\n[1] \"GCAAGGCCCAATGGGGCGGCGGGCCCGGCAGCCCCGCCCCGGTGGTGTCCGCGCGGCCCG 60\"\n[1] \" \"<\/mark><\/em>\n.....\n.....\n18531\"\n[1] \"TGTTTTACCTCAGCCAGTCAGTATGTTTTGCATACTGTAATTTATATTGCCCTTGGAACA 18386\"\n[1] \" \"\n[1] \"CATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA 18591\"\n[1] \"CATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA 18446\"\n[1] \" \"<\/mark><\/em><\/code><\/pre>\n\n\n\n
This looks promising. To create a new dot plot with the aligned sequences:<\/p>\n\n\n\n
comparison@pattern # this is the human\n<\/em><\/mark>[7] GCAAGGCCCAATGGGGCGGCAGGCCCGGCAGCCCCGCCCCGGTGGTGCCCGCGCGGCCAGCGCCCGCCAGGCCCAGCGT...TATGTTTTGCATACTGTAATTTATATTGCCCTTGGAACACATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA <\/mark><\/em>\ncomparison@subject # this is the orangutan\n<\/em><\/mark>[1] GCAAGGCCCAATGGGGCGGCGGGCCCGGCAGCCCCGCCCCGGTGGTGTCCGCGCGGCCCGCGCCCGCCAGGCCCAGCGT...TATGTTTTGCATACTGTAATTTATATTGCCCTTGGAACACATGGTGCCATATTTAGCTACTAAAAAGCTCTTCACAAAA <\/mark><\/em>\n\nhuman_seq <- as.character(comparison@pattern)\norangutan_seq <- as.character(comparison@subject)\nhuman_seq <- getSequence(human_seq)\norangutan_seq <- getSequence(orangutan_seq)\ndotPlot(human_seq, orangutan_seq, wsize=9,wstep=9,nmatch=9,col=c('white', 'black'))<\/em><\/mark><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
This looks much better and the DNA sequences look very similar. However, the dot plot function included in the seqinr9<\/a><\/sup> package is not very versatile. It would be better to plot the comparison with ggplot210<\/a><\/sup> .<\/p>\n\n\n\n
The main advantage of open source software is that the code of 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) rather than a plot. 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 and execute it.<\/p>\n\n\n\n
The new function (dot_plot_data) should now be available to R (it can also be viewed and downloaded from the download page<\/a>). To access the function comparing nine nucleotides and matching all nine nucleotides takes the same format as the dotPlot function from the seqinr package11<\/a><\/sup>:<\/p>\n\n\n\n
new_plot <- dot_plot_data(human_seq, orangutan_seq, wsize = 9, wstep = 9, nmatch = 9)<\/span><\/em><\/code><\/pre>\n\n\n\n
To show the structure (str) and class of the returned object:<\/p>\n\n\n\n
str(new_plot)<\/em><\/span>\n logi [1:2071, 1:2071] TRUE FALSE FALSE FALSE FALSE FALSE ...<\/em><\/span>\n - attr(*, \"dimnames\")=List of 2<\/em><\/span>\n ..$ : chr [1:2071] \"GCAAGGCCC\" \"AATGGGGCG\" \"GCAGGCCCG\" \"GCAGCCCCG\" ...<\/em><\/span>\n ..$ : chr [1:2071] \"GCAAGGCCC\" \"AATGGGGCG\" \"GCGGGCCCG\" \"GCAGCCCCG\" ...<\/em><\/span>\nclass(new_plot)<\/em><\/span>\n[1] \"matrix\" \"array\"<\/em><\/span><\/code><\/pre>\n\n\n\n
The data are now available as a matrix called new_plot. The row names are sequence 1 and the column names sequence 2. Unfortunately, ggplot2 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)<\/em><\/span>\nsequence_y <- colnames(new_plot)<\/em><\/span><\/code><\/pre>\n\n\n\n
Secondly, remove the row and column names from the matrix (the names are not unique and otherwise ggplot2 will group by them):<\/p>\n\n\n\n
rownames(new_plot) <- NULL<\/span>\ncolnames(new_plot) <- NULL<\/span><\/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
library(reshape2<\/mark>)<\/mark><\/em>\nnew_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 <- new_plot_2 %>%<\/span><\/em>\nggplot(aes(x = Var1, y = Var2, fill = value)) +<\/span><\/em>\ngeom_raster() +<\/span><\/em>\nscale_fill_manual(values = c(\"lightgray\", \"blue\"), name = \"Pair\", <\/span><\/em>\n  labels = c(\"Different\", \"Same\")) +<\/span><\/em>\nscale_x_continuous(\"Human\") +<\/span><\/em>\nscale_y_continuous(\"Orangutan\") + <\/span><\/em>\ntheme_bw() <\/span><\/em>\n\ndev.new()<\/span><\/em>\ndot_plot<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
If the sequences would be indicated on the axes, the axes would become too crowded. However, the plot can be windowed (here between 10 and 20):<\/p>\n\n\n\n
dot_plot_2 <- new_plot_2 %>%<\/span><\/em>\nggplot(aes(x = Var1, y = Var2, fill = value)) +<\/span><\/em>\ngeom_raster() +<\/span><\/em>\nscale_fill_manual(values = c(\"lightgray\", \"blue\"), name = \"Pair\",<\/span><\/em>\n  labels = c(\"Different\", \"Same\")) +<\/span><\/em>\nscale_x_continuous(\"Human\", breaks = 1: length(sequence_x), <\/span><\/em>\n  labels = sequence_x, limit = c(10,20)<\/strong>) +<\/span><\/em>\nscale_y_continuous(\"Orangutan\", breaks = 1: length(sequence_y), <\/span><\/em>\n  labels = sequence_y, limit = c(10,20)<\/strong>) + <\/span><\/em>\ntheme_bw() <\/span><\/em>\n\ndev.new()<\/span><\/em>\ndot_plot_2<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n","protected":false},"excerpt":{"rendered":"
Please make sure the necessary packages (seqinr and Biostrings) are installed as described to allow analysis. The Biostrings and pwalign packages are part of Bioconductor and installation is a little different: If you used reshape2 or dplyr, you may need to restart R if you get a can’t unload package error message. Furthermore, the ggplot2 and reshape2 […]<\/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-2232","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2232","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=2232"}],"version-history":[{"count":8,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2232\/revisions"}],"predecessor-version":[{"id":4566,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/2232\/revisions\/4566"}],"wp:attachment":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/media?parent=2232"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}