Bioinformatics Final Project

Nupur Joshi (nkj219) - Fall 2019

1. Repeat questions 5 and 7 from Part 1 of the term project, except using the command-line version of hmmer.

#First, I changed my directory to make sure I am pulling from and saving to the right folder. Then...

phmmer Human_PUSL1.fasta uniprot_sprot.fasta

/var/folders/pz/rk934mt528g4zlchvtw8tty80000gn/T/Rtmpjzmj9M/chunk-code1186e391cb817.: line 3: phmmer: command not found

For Q5 of the midterm, I used HMMER in R and BLAST to find the best matches to my assigned gene. With R, I got a bacterial match with an E-value of 3.6e-24. With BLAST, my top match was a human one with an E-value of 0. With command-line, my top match is a human one with an E-value of 8.5e-203 and my second match is a mouse one with an E-value of 4.1e-153. Quite the contrast because, with BLAST, my mouse results didn’t even show up in the first 100. Please see accompanying figure 1.

phmmer Mouse_PUSL1.fasta uniprot_sprot.fasta

/var/folders/pz/rk934mt528g4zlchvtw8tty80000gn/T/Rtmpjzmj9M/chunk-code1186e207bb0fa.: line 1: phmmer: command not found

Unfortunately, on the midterm, I did not properly complete the HMMER and BLAST steps for the mouse gene so lost marks as a result. To complete the comparison between command-line and BLAST, I conducted a BLAST search for it for the final. I have found that the top, not only 2, but 5 results are mouse hits with E-values of 0. With command-line, however, the top hit is a mouse one with an E-value of 2.9e-195 and the second hit is a human one with an E-value of 4.4e-153. Please see accompanying figure 2.

For Q7 of the midterm, I used R to build the MSA and profile HMM for the given fly protein. Here, I used command-line. The portions done in clustalw are included as screenshots. Please refer to accompanying figures 3-12.

hmmbuild fly_protein.hmm fly_protein_new
hmmsearch fly_protein.hmm uniprot_sprot.fasta

/var/folders/pz/rk934mt528g4zlchvtw8tty80000gn/T/Rtmpjzmj9M/chunk-code1186e2fcc133a.: line 1: hmmbuild: command not found
/var/folders/pz/rk934mt528g4zlchvtw8tty80000gn/T/Rtmpjzmj9M/chunk-code1186e2fcc133a.: line 2: hmmsearch: command not found

The most significant human hit is ZMAT3_HUMAN (Q9HA38, Zinc finger matrin-type protein 3 O) with an E-value of 4.6e-12. Please see accompanying figure 13.

#Adding all libraries for the R section of the project.
library("HMM")
library(seqinr)
library("ape")
library(Biostrings)
library(phangorn)
library("bio3d")
library("msa")

2. The HMM model for 5’ splice site detection we looked at in class was a “toy” HMM model to use for splice site detection that would not be expected to perform well: it only used the information that the first base of an intron is likely to be a G, and the emission probabilities were not empirically derived. It is known that the full motif for U2 (major class) introns in pre-mRNA has the following 5’ splice site consensus sequence MAG|GTRAGT (using an ambiguity code, where | is the exon-intron boundary), and so including additional upstream and downstream base states could also improve the model (). Improve the model given in class to use the information (i.e. add states for) for two upstream and two downstream bases, and use the following empirically derived exon and intron emission probabilities. Exon base frequencies: A=0.20,C=0.30,G=0.30,T=0.20. Intron base frequencies: A=0.15,C=0.35,G=0.35,T=0.15 Motif base 1: A=0.997,C=0.001,G=0.001,T=0.001, Motif base 2: A=0.001,C=0.001,G=0.997,T=0.001, Motif base 3: A=0.001,C=0.001,G=0.997,T=0.001, Motif base 4: A=0.001,C=0.001,G=0.001,T=0.997, Use the hmm() R package to define an improved 5’ exon-intron junction detector using this information. Then use the UCSC genome browser to extract the human hg38 DNA sequence around the final exon-intron splice site for your assigned human gene (50 bp upstream and downstream); and also extract the following genomic region chr2:85,539,313-85,539,468 (on the + strand). Now try to detect the exon-intron boundary using your improved HMM and compare with the previous HMM given in class for both sequences. Compare with the gencode gene annotation in the UCSC browser- was the correct splice site location detected in the two sequences?

#My partner for this question is Zaina Atieh.
#My New Junction Detector
TheUnitedStates = c("Exon", "toA","toG1","toG2", "toT","Intron")
StarsAndStripes = c("A","C","G","T")
transProbs = matrix(c('E_toE'=0.9,'E_toA'=0.1, 'E_toG1'= 0, 'E_toG2'= 0, 'E_toT'= 0,'E_toI'=0, 'A_toE'=0,'A_toA'=0, 'A_toG1'=1.0, 'A_toG2'=0, 'A_toT'=0, 'A_toI'=0, 'G1_toE'=0,'G1_toA'=0,'G1_toG1'=0,'G1_toG2'=1.0,'G1_toT'=0,'G1_toI'=0,'G2_toE'=0,'G2_toA'=0,'G2_toG1'=0,'G2_toG2'=0,'G2_toT'=1.0,'G2_toI'=0,'T_toE'=0,'T_toA'=0,'T_toG1'=0,'T_toG2'=0,'T_toT'=0,'T_toI'=1.0, 'I_toE'=0, 'I_toA'=0,'I_toG1'=0,'I_toG2'=0,'I_toT'=0, 'I_toI'=1.0), c(length(TheUnitedStates), length(TheUnitedStates)), byrow = TRUE)
transProbs

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  0.9  0.1    0    0    0    0
[2,]  0.0  0.0    1    0    0    0
[3,]  0.0  0.0    0    1    0    0
[4,]  0.0  0.0    0    0    1    0
[5,]  0.0  0.0    0    0    0    1
[6,]  0.0  0.0    0    0    0    1

emissionProbs = matrix(c('A'=0.20,'C'=0.30,'G'=0.30,'T'=0.20, 'A'=0.997,'C'=0.001,'G'=0.001,'T'=0.001,'A'=0.001,'C'=0.001,'G'=0.997,'T'=0.001,'A'=0.001,'C'=0.001,'G'=0.997,'T'=0.001,'A'=0.001,'C'=0.001,'G'=0.001,'T'=0.997,'A'=0.15,'C'=0.35,'G'=0.35,'T'=0.15), c(length(TheUnitedStates), length(StarsAndStripes)), byrow = TRUE)
emissionProbs

      [,1]  [,2]  [,3]  [,4]
[1,] 0.200 0.300 0.300 0.200
[2,] 0.997 0.001 0.001 0.001
[3,] 0.001 0.001 0.997 0.001
[4,] 0.001 0.001 0.997 0.001
[5,] 0.001 0.001 0.001 0.997
[6,] 0.150 0.350 0.350 0.150

hmm <- initHMM(TheUnitedStates, StarsAndStripes, startProbs = c(1,0,0,0,0,0), transProbs = transProbs, emissionProbs = emissionProbs)
#Class Junction Detector
ParkerStates = c("Exon", "5site", "Intron")
ParkerSymbols = c("A","C","G","T")
ParkertransProbs = matrix(c('EE'=0.9,'E5'=0.1,'EI'=0, '5E'=0, '55'=0, '5I'=1.0, 'IE'=0, 'I5'=0, 'II'=1.0), c(length(ParkerStates), length(ParkerStates)), byrow = TRUE)
ParkertransProbs

     [,1] [,2] [,3]
[1,]  0.9  0.1    0
[2,]  0.0  0.0    1
[3,]  0.0  0.0    1

ParkeremissionProbs = matrix(c('A'=0.25,'C'=0.25,'G'=0.25,'T'=0.25, 'A'=0.05,'C'=0.0,'G'=0.95,'T'=0.0, 'A'=0.4,'C'=0.1,'G'=0.1,'T'=0.4), c(length(ParkerStates), length(ParkerSymbols)), byrow = TRUE)
ParkeremissionProbs

     [,1] [,2] [,3] [,4]
[1,] 0.25 0.25 0.25 0.25
[2,] 0.05 0.00 0.95 0.00
[3,] 0.40 0.10 0.10 0.40

Parkerhmm <- initHMM(ParkerStates, ParkerSymbols, startProbs = c(1,0,0), transProbs = ParkertransProbs, emissionProbs = ParkeremissionProbs)
#PUSL1 Human
BaldEagle <- s2c('AGCCCCAGCCCACGGCTTATTCCTCAAGTCAGTGCTGTACGGGAACCTCGgtaagaaaaacaggcacgagaagctcctgtcatgtgcccagtgactactg')
BaldEagle <- toupper(BaldEagle)
vitBaldEagle <- viterbi(hmm, BaldEagle)
vitBaldEagle

  [1] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [13] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [25] "Exon"   "toA"    "toG1"   "toG2"   "toT"    "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [37] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [49] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [61] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [73] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [85] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [97] "Intron" "Intron" "Intron" "Intron"

ParkervitBaldEagle <- viterbi(Parkerhmm, BaldEagle)
ParkervitBaldEagle

  [1] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [17] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [33] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [49] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [65] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [81] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [97] "Exon" "Exon" "Exon" "Exon"

#chr2:85,539,313-85,539,468
NorthAmericanBison <- s2c('ACGAGGCGTTCATCGAGGAGGGCACATTCCTTTTCACCTCAGAGTCGGTCGGGGAAGGCCACCCAGGTGAGGGGACGGCCTGAAGCGAAGCGTGGGGCGGGGCAGAAGGCAGCGCCAAGGTCCGGCTGGCTGCGGCCGGCCGGTGGTGGGGCCCGC')
vitNorthAmericanBison <- viterbi(hmm, NorthAmericanBison)
vitNorthAmericanBison

  [1] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [13] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [25] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [37] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [49] "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"   "Exon"  
 [61] "Exon"   "Exon"   "Exon"   "Exon"   "toA"    "toG1"   "toG2"   "toT"    "Intron" "Intron" "Intron" "Intron"
 [73] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [85] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
 [97] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
[109] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
[121] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
[133] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"
[145] "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron" "Intron"

ParkervitNorthAmericanBison <- viterbi(Parkerhmm, NorthAmericanBison)
ParkervitNorthAmericanBison

  [1] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [17] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [33] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [49] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [65] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [81] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
 [97] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
[113] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
[129] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"
[145] "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon" "Exon"

As expected, the improved junction detector did return a splice site while the previous junction detector we used in class did not. This occurred in both the sequence for the PUSL1 Human gene assigned to me as well as the sequence at chr2:85,539,313-85,539,468. In the latter, the improved detector picked up on the exact splice site found in the UCSC genome browser. However, in the former, with my own gene, the detector picked up something about 25 bases upstream of the actual splice site seen in the browser.

3. In this question you will examine a simple toy phylogenetic tree, and associated MSA, similar to the one used as a maximum likelihood example in lectures. In the file toy_MSA.fasta is defined a (very simple) multiple sequence alignment (MSA) of a region of DNA of length 1 across four primate species. In the file toy_tree.tre is defined a potential phylogenetic tree in newick format for which we want to determine its fit to the MSA. Read in the tree (using read.tree() in R) and also read in the MSA. Display the unrooted tree. Assume that in this simple (unrealistic) example all tree branches have the same length t, and assume base composition is equal for all bases. Compute the log likelihood of the tree by hand, using the Jukes-Cantor evolutionary model and for parameter values 𝛼t = 0.09, 0.103, and 0.12. Which value of 𝛼t shows the maximum likelihood? Show your working. Then compute the log likelihood of the tree defined in toy_tree.tre using the R phangorn package and compare your results.

ChristmasTree <- read.tree("toy_tree.tre")
SantasToys = readDNAMultipleAlignment("toy_MSA.fasta", format="FASTA")
SantasToys <- as.phyDat(SantasToys)
plot(ChristmasTree, type="unrooted")

fit = pml(ChristmasTree, data=SantasToys, model="Jukes-Cantor")
fit

 loglikelihood: -4.778555 

unconstrained loglikelihood: 0 

Rate matrix:
  a c g t
a 0 1 1 1
c 1 0 1 1
g 1 1 0 1
t 1 1 1 0

Base frequencies:  
0.25 0.25 0.25 0.25 

logLik(fit)

'log Lik.' -4.778555 (df=5)

Reindeer<- function(a){
  Rudolph<-(1+3*exp(-4*a))/4
  Vixen<-(1-exp(-4*a))/4
  Prancer<-(7*Vixen^5+2*Rudolph^3*Vixen^2+Rudolph^4*Vixen+4*Rudolph^2*Vixen^3+2*Rudolph*Vixen^4)/4
  Blitzen<-log(Prancer)
  return(Blitzen)
}
Reindeer(0.09)

[1] -4.78537

Reindeer(0.103)

[1] -4.778555

Reindeer(0.12)

[1] -4.787469

Of the three values, 0.103 shows the maximum likelihood which is -4.778555. The log likelihood is also computed using phangorn above the manual calculations. The result of this is also -4.778555.

4. Using PFAM, extract the multiple sequence seed alignment for the most significant protein domain from your assigned human gene (if no known domain was found in your gene, use the piwi domain). Using phanghorn, compute the gene tree from this MSA using both neighbour- joining and maximum parsimony and maximum likelihood. Show the parsimony score and log likelihood of your maximum parsimony and maximum likelihood trees, respectively. Show your R code.

#Extract the multiple sequence seed alignment for the most significant protein domain from PFAM.
seedalign <- pfam("PseudoU_synth_1", alignment = "seed") 
PseudoU <- as.phyDat(seedalign$ali, type="AA")
#phangorn part
#neighbor-joining
dm <- dist.ml(PseudoU)
treeNJ <- NJ(dm)
plot(treeNJ, "unrooted", main="NJ")

#maximum parsimony
treePars <- optim.parsimony(treeNJ, PseudoU,method = "sankoff")

Final p-score 11170 after  89 nni operations 

plot(treePars, "unrooted", main="parsimony")

parsimony(treeNJ, PseudoU,method = "sankoff")

[1] 11337

The parsimony score is 11337.

#maximum likelihood
fitpars <- pml(treeNJ, data=PseudoU)

negative edges length changed to 0!

fitpars

 loglikelihood: -59152.72 

unconstrained loglikelihood: -1249.257 
Rate matrix: 

logLik(fitpars)

'log Lik.' -59152.72 (df=619)

The log likelihood is -59152.72.

plot(fitpars$tree)

5. The SARS epidemic of 2003 was a viral respiratory disease caused by a mutation in a coronavirus. Various samples of the virus were sequenced as the virus mutated and spread throughout the world. It is also known that the SARS coronavirus is distantly related to the Himalayan Palm Civet coronavirus and so this can be used as an outgroup in our phylogenetic analysis. The sequences of a key protein in each of these samples are provided in the fasta file “sars.fasta”. Perform a phylogenetic analysis of the spread of the disease from these samples to determine the likely site of origin of the disease. Show a NJ phylogenetic tree as your evidence and give your arguments for (i) the site of origin, (ii) the next major site SARS spread to. You should use the phangorn package in R. Note that the root function can be used to set the root of the tree e.g. root(treeNJ, outgroup = “Himalayan palm civet sars cov, complete genome”) .

SARS <- read.phyDat("sars.fasta", type="DNA", format = "fasta")
SARSdm <- dist.ml(SARS)
treeNJ <- NJ(SARSdm)
treeNJroot <- root(treeNJ, outgroup = "Himalayan palm civet sars cov, complete genome", resolve.root=TRUE)
plot(treeNJroot, type="phylogram")

Based on this, we can deduce that the city of origin is Guangzhou in Guangdong. Depending on the definition of the next major site, Zhongshan is the next major city but still in Guangdong and Hong Kong is the next major region after Guangdong as a whole.