{"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 (seqinr and Biostrings) 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 ggplot2 and reshape2 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 reshape2:<\/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 ggplot2:<\/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 package 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 Coghlan 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 seqinr package is not very versatile. It would be better to plot the comparison with ggplot2 .<\/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 package:<\/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 package 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}]}}