# Lab 7 for the University of Tulsa's CS-6643 Bioinformatics Course
# PDB
# Professor: Dr. McKinney, Fall 2022
# Noah L. Schrick - 1492657

## Set Working Directory to file directory - RStudio approach
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#### Part A: Obtaining PDB - no supporting R Code

#### Part B: Visualize the 3D structure
## Install Rpdb and load the pdb
if (!require("Rpdb")) install.packages("Rpdb")
library(Rpdb)
x<-read.pdb("1TGH.pdb")
natom(x)
visualize(x,type="l")

## Visualize the B and C chains
B_chain_pdb <- subset(x$atoms, x$atoms$chainid=="B") 
C_chain_pdb <- subset(x$atoms, x$atoms$chainid=="C") 
# remove water:
C_chain_pdb <- subset(C_chain_pdb,C_chain_pdb$resname!="HOH") 

# visualize chains B and C
BC_chains_pdb <- subset(x$atoms, x$atoms$chainid=="B" | 				                             x$atoms$chainid=="C")
color.vec <- c(rep("red",natom(B_chain_pdb)),rep("green",natom(C_chain_pdb)))

visualize(BC_chains_pdb,col=color.vec)
addResLab(BC_chains_pdb)

## Visualize B-C and A Chains
A_chain_pdb <- subset(x$atoms, x$atoms$chainid=="A") 
# remove water
A_chain_pdb <- subset(A_chain_pdb, A_chain_pdb$resname!="HOH")  

# visualize complex complex
BCA_chains_pdb <- subset(x$atoms, x$atoms$chainid=="B" |  
                           x$atoms$chainid=="C" | x$atoms$chainid=="A")
BCA.color.vec <- c(rep("red",natom(B_chain_pdb)),rep("green",natom(C_chain_pdb)),rep("blue",natom(A_chain_pdb)))

visualize(BCA_chains_pdb,col=BCA.color.vec)

#### Part C: Primary structure and DNA Palindromes
# get coordinates of C1' atoms of the C-chain DNA molecule
C_chain_pdb$resname
C_chain_resids<-unique(C_chain_pdb$resid)
C_chain_C1prime <- subset(C_chain_pdb, C_chain_pdb$elename=="C1'")

# get chain C DNA sequence
C_chain_sequence_messy <- C_chain_C1prime$resname 
C_chain_sequence <- paste(sapply(C_chain_sequence_messy,function(x) {unlist(strsplit(x,""))[2]}),collapse = "")

## Find palindromes
if (!require("BiocManager")) install.packages("BiocManager")
library(BiocManager)
if (!require("Biostrings")) BiocManager::install("Biostrings")
library(snpStats)

C_chain_DNAString <- DNAString(C_chain_sequence)
dna.pals <- findPalindromes(C_chain_DNAString, min.armlength=3, 
                            max.looplength=5, max.mismatch = 0)