Introduction

This code compiles summary information about the gene GYPA (Glycophorin A). This protein is that this gene codes for is associated with diseases such as Malaria and Hepatitis A.

It also generates alignments and a phylogenetic tree to indicate the evolutionary relationship between the human version of the gene and its homologs in other species.

Resources / References

Key information use to make this script can be found here: - Refseq Gene: https://www.ncbi.nlm.nih.gov/gene/2993 - Refseq Homologene: https://www.ncbi.nlm.nih.gov/homologene/48076

Other resources consulted includes - Neanderthal genome: http://neandertal.ensemblgenomes.org/index.html

Other interesting resources and online tools include: - Sub-cellular locations prediction: https://wolfpsort.hgc.jp/

Preparation

Load necessary packages:

Download and load drawProteins from Bioconductor

library(BiocManager)

## Bioconductor version '3.13' is out-of-date; the current release version '3.14'
##   is available with R version '4.1'; see https://bioconductor.org/install

#install("drawProteins")

library(drawProteins)

Load other packages

# github packages
library(compbio4all)
library(ggmsa)

## Registered S3 methods overwritten by 'ggalt':
##   method                  from   
##   grid.draw.absoluteGrob  ggplot2
##   grobHeight.absoluteGrob ggplot2
##   grobWidth.absoluteGrob  ggplot2
##   grobX.absoluteGrob      ggplot2
##   grobY.absoluteGrob      ggplot2

# CRAN packages
library(rentrez)
library(seqinr)
library(ape)

## 
## Attaching package: 'ape'

## The following objects are masked from 'package:seqinr':
## 
##     as.alignment, consensus

library(pander)


library(ggplot2)

# Bioconductor packages
## msa
### The msa package is having problems on some platforms
### You can skip the msa steps if necessary.  The msa output
### is used to make a distance matrix and then phylogenetics trees,
### but I provide code to build the matrix by hand so
### you can proceed even if msa doesn't work for you.

## Biostrings
library(Biostrings)

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min

## Loading required package: S4Vectors

## Loading required package: stats4

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## Loading required package: XVector

## Loading required package: GenomeInfoDb

## 
## Attaching package: 'Biostrings'

## The following object is masked from 'package:ape':
## 
##     complement

## The following object is masked from 'package:seqinr':
## 
##     translate

## The following object is masked from 'package:base':
## 
##     strsplit

library(msa)

## 
## Attaching package: 'msa'

## The following object is masked from 'package:BiocManager':
## 
##     version

library(drawProteins)



library(HGNChelper)

Accession numbers

Accession numbers were obtained from RefSeq, Refseq Homlogene, UniProt and PDB. UniProt accession numbers can be found by searching for the gene name. PDB accessions can be found by searching with a UniProt accession or a gene name, though many proteins are not in PDB.

A protein BLAST search (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome) was carried out to find more species, since the HomoloGene search only yielded 2 species. The gene not appears to be primate-only.

OPTIONAL: Use the function to confirm the validity of your gene name and any aliases

# this is optional
HGNChelper::checkGeneSymbols(x = c("GYPA"))

## Maps last updated on: Thu Oct 24 12:31:05 2019

##      x Approved Suggested.Symbol
## 1 GYPA     TRUE             GYPA

Accession number table

Not available: - Drosophila Does not occur: - Outside of primates

#                       RefSeq           Uniprot     PDB    sci name                common name          gene name
GYPA_table_vector<-c("NP_002090.4",     "P02724", "1AFO",  "Homo sapiens" ,           "Human",               "GYPA",
                     "XP_001145917.1",  "NA",     "NA",    "Pan troglodytes" ,        "Chimpanzee",          "GYPA",
                     "XP_003815879.2",  "NA",     "NA",    "Pan paniscus",            "Bonobo",              "GYPA",
                     "XP_021794517.1",  "NA",     "NA",    "Papio anubis ",           "Baboon",              "GYPA",
                     "XP_011842541.1",  "NA",     "NA",    "Mandrillus leucophaeus",  "Drill",               "GYPA",
                     "XP_009238620.2",  "NA",     "NA",    "Pongo abelii",            "Sumatran orangutan",  "GYPA",
                     "XP_014994775.2",  "NA",     "NA",    "Macaca mulatta",          "Rhesus macaque",      "GYPA",
                     "XP_037848439.1",  "NA",     "NA",    "Chlorocebus sabaeus",     "Green monkey",        "GYPA",
                     "XP_032027955.1",  "NA",     "NA",    "Hylobates moloch",        "Silvery gibbon",      "GYPA",
                     "XP_025241177.1",  "NA",     "NA",    "Theropithecus gelada",    "Gelada",              "GYPA")

GYPA_matrix <- matrix( GYPA_table_vector, ncol = 6, byrow = TRUE)
GYPA_df <- data.frame( GYPA_matrix )

colnames( GYPA_df ) <- c("ncbi.protein.accession", "UniProt.id", "PDB", "species", "common.name",
                          "gene.name")

The finished table

pander::pander( GYPA_df )

Table continues below
ncbi.protein.accession	UniProt.id	PDB	species
NP_002090.4	P02724	1AFO	Homo sapiens
XP_001145917.1	NA	NA	Pan troglodytes
XP_003815879.2	NA	NA	Pan paniscus
XP_021794517.1	NA	NA	Papio anubis
XP_011842541.1	NA	NA	Mandrillus leucophaeus
XP_009238620.2	NA	NA	Pongo abelii
XP_014994775.2	NA	NA	Macaca mulatta
XP_037848439.1	NA	NA	Chlorocebus sabaeus
XP_032027955.1	NA	NA	Hylobates moloch
XP_025241177.1	NA	NA	Theropithecus gelada

common.name	gene.name
Human	GYPA
Chimpanzee	GYPA
Bonobo	GYPA
Baboon	GYPA
Drill	GYPA
Sumatran orangutan	GYPA
Rhesus macaque	GYPA
Green monkey	GYPA
Silvery gibbon	GYPA
Gelada	GYPA

Data Preparation

Download Sequences

All sequences were downloaded using a wrapper compbio4all::entrez_fetch_list() which uses rentrez::entrez_fetch() to access NCBI databases.

# download FASTA sequences
GYPA_list <- compbio4all::entrez_fetch_list( db = "protein",
                                              id = GYPA_df$ncbi.protein.accession,
                                              rettype = "fasta"
                                              )

Number of FASTA files obtained

length( GYPA_list )

## [1] 10

The first entry

GYPA_list[[1]]

## [1] ">NP_002090.4 glycophorin-A isoform 1 precursor [Homo sapiens]\nMYGKIIFVLLLSEIVSISALSTTEVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTV\nYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSV\nEIENPETSDQ\n\n"

# output should be the FASTA sequence with header information and newlines still included

Initial Data Cleaning

Remove FASTA header

for(i in 1:length(GYPA_list)){
  GYPA_list[[i]] <- compbio4all::fasta_cleaner(GYPA_list[[i]], parse = F)
}

Specific additional cleaning steps will be as needed for particular analyses

General Protein Information

Protein Diagram

For code see https://rpubs.com/lowbrowR/drawProtein

P02724_json <- drawProteins::get_features("P02724")

## [1] "Download has worked"

my_prot_df <- drawProteins::feature_to_dataframe(P02724_json)
is(my_prot_df)

## [1] "data.frame"       "list"             "oldClass"         "vector"          
## [5] "list_OR_List"     "vector_OR_Vector" "vector_OR_factor"

my_canvas <- draw_canvas(my_prot_df)  
my_canvas <- draw_chains(my_canvas, my_prot_df, 
                         label_size = 2.5)
my_canvas <- draw_domains(my_canvas, my_prot_df)
my_canvas

Draw dotplot

Prepare Data

GYPA_list[[1]]

## [1] "MYGKIIFVLLLSEIVSISALSTTEVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ"

GYPA_human_vector <- unlist(strsplit( GYPA_list[[1]], "" ))
seqinr::dotPlot( GYPA_human_vector, GYPA_human_vector )

TODO:

Create a 2x2 panel to show different values for dotPlots settings
Create a single large plot with the best version of the plot

par(mfrow = c(2,2), 
    mar = c(0,0,2,1))

# plot 1: Defaults
seqinr::dotPlot(GYPA_human_vector, GYPA_human_vector, 
        wsize = 1, 
        nmatch = 1, 
        main = "size=1, num match=1")

# plot 2 size = 10, nmatch = 10
seqinr::dotPlot(GYPA_human_vector, GYPA_human_vector, 
        wsize = 10, 
        nmatch = 1, 
        main = "size = 10, nmatch = 10")

# plot 3: size = 10, nmatch = 5
seqinr::dotPlot(GYPA_human_vector, GYPA_human_vector, 
        wsize = 10, 
        nmatch = 5, 
        main = "size = 10, nmatch = 5")

# plot 4: size = 20, nmatch = 5
seqinr::dotPlot(GYPA_human_vector, GYPA_human_vector, 
        wsize = 20,
        nmatch = 5,
        main = "size = 20, nmatch = 5")

par(mfrow = c(1,1), 
    mar = c(4,4,4,4))

seqinr::dotPlot(GYPA_human_vector, GYPA_human_vector, 
        wsize = 20,
        nmatch = 5,
        main = "GYPA human dot plot")

Protein properties compiled from databases

TODO: Create table

Below are links to relevant information. 1. Pfam: http://pfam.xfam.org/protein/P02724; mostly “Glycophorin A”, transmembrane region from 92 to 114 2. DisProt: NA 3. RepeatDB: NA 4. PDB secondary structural location: NA

Protein feature prediction

Predict protein fold

NOTE: My protein does NOT contain “U”.

First, I need the data from Chou and Zhang (1994) Table 5. Code to build this table is available at https://rpubs.com/lowbrowR/843543

The table looks like this:

# enter once
aa.1.1 <- c("A","R","N","D","C","Q","E","G","H","I",
            "L","K","M","F","P","S","T","W","Y","V")
# alpha proteins
alpha <- c(285, 53, 97, 163, 22, 67, 134, 197, 111, 91, 
           221, 249, 48, 123, 82, 122, 119, 33, 63, 167)
# beta proteins
beta <- c(203, 67, 139, 121, 75, 122, 86, 297, 49, 120, 
          177, 115, 16, 85, 127, 341, 253, 44, 110, 229)
# alpha + beta
a.plus.b <- c(175, 78, 120, 111, 74, 74, 86, 171, 33, 93,
              110, 112, 25, 52, 71, 126, 117, 30, 108, 123)
# alpha/beta
a.div.b <- c(361, 146, 183, 244, 63, 114, 257, 377, 107, 239, 
             339, 321, 91, 158, 188, 327, 238, 72, 130, 378)

pander(data.frame(aa.1.1, alpha, beta, a.plus.b, a.div.b))

aa.1.1	alpha	beta	a.plus.b	a.div.b
A	285	203	175	361
R	53	67	78	146
N	97	139	120	183
D	163	121	111	244
C	22	75	74	63
Q	67	122	74	114
E	134	86	86	257
G	197	297	171	377
H	111	49	33	107
I	91	120	93	239
L	221	177	110	339
K	249	115	112	321
M	48	16	25	91
F	123	85	52	158
P	82	127	71	188
S	122	341	126	327
T	119	253	117	238
W	33	44	30	72
Y	63	110	108	130
V	167	229	123	378

Convert to frequencies Table 5 therefore becomes this

alpha.prop <- alpha/sum(alpha)
beta.prop <- beta/sum(beta)
a.plus.b.prop <- a.plus.b/sum(a.plus.b)
a.div.b <- a.div.b/sum(a.div.b)

aa.prop <- data.frame(alpha.prop,
                      beta.prop,
                      a.plus.b.prop,
                      a.div.b)
row.names(aa.prop) <- aa.1.1

pander::pander(aa.prop)

	alpha.prop	beta.prop	a.plus.b.prop	a.div.b
A	0.1165	0.07313	0.09264	0.08331
R	0.02166	0.02414	0.04129	0.03369
N	0.03964	0.05007	0.06353	0.04223
D	0.06661	0.04359	0.05876	0.05631
C	0.008991	0.02702	0.03917	0.01454
Q	0.02738	0.04395	0.03917	0.02631
E	0.05476	0.03098	0.04553	0.05931
G	0.08051	0.107	0.09052	0.08701
H	0.04536	0.01765	0.01747	0.02469
I	0.03719	0.04323	0.04923	0.05516
L	0.09031	0.06376	0.05823	0.07824
K	0.1018	0.04143	0.05929	0.07408
M	0.01962	0.005764	0.01323	0.021
F	0.05027	0.03062	0.02753	0.03646
P	0.03351	0.04575	0.03759	0.04339
S	0.04986	0.1228	0.0667	0.07547
T	0.04863	0.09114	0.06194	0.05493
W	0.01349	0.01585	0.01588	0.01662
Y	0.02575	0.03963	0.05717	0.03
V	0.06825	0.08249	0.06511	0.08724

Determine the number of each amino acid in my protein.

A Function to convert a table into a vector is helpful here because R is goofy about tables not being the same as vectors.

table_to_vector <- function(table_x){
  table_names <- attr(table_x, "dimnames")[[1]]
  table_vect <- as.vector(table_x)
  names(table_vect) <- table_names
  return(table_vect)
}

GYPA_human_table <- table(GYPA_human_vector)/length(GYPA_human_vector)
GYPA.human.aa.freq <- table_to_vector(GYPA_human_table)
GYPA.human.aa.freq

##          A          D          E          F          G          H          I 
## 0.04666667 0.04000000 0.08666667 0.02000000 0.04000000 0.03333333 0.10000000 
##          K          L          M          N          P          Q          R 
## 0.04000000 0.07333333 0.02000000 0.01333333 0.06666667 0.02000000 0.04000000 
##          S          T          V          Y 
## 0.14000000 0.10000000 0.08666667 0.03333333

Check for the presence of “U” (unknown aa.)

aa.names <- names(GYPA.human.aa.freq)
i.U <- which(aa.names == "U")
aa.names[i.U]

## character(0)

GYPA.human.aa.freq[i.U]

## named numeric(0)

GYPA does not have any occurrences of “U”

Add data on my focal protein to the amino acid frequency table.

# needed to do some work because my gene does not have C or W
new.GYPA.human.aa.freq <- append(unname(GYPA.human.aa.freq), c(0,0))
names(new.GYPA.human.aa.freq) <- append(names(GYPA.human.aa.freq), c("C", "W"))


new.GYPA.human.aa.freq

##          A          D          E          F          G          H          I 
## 0.04666667 0.04000000 0.08666667 0.02000000 0.04000000 0.03333333 0.10000000 
##          K          L          M          N          P          Q          R 
## 0.04000000 0.07333333 0.02000000 0.01333333 0.06666667 0.02000000 0.04000000 
##          S          T          V          Y          C          W 
## 0.14000000 0.10000000 0.08666667 0.03333333 0.00000000 0.00000000

aa.prop$GYPA.human.aa.freq <- new.GYPA.human.aa.freq
pander::pander(aa.prop)

	alpha.prop	beta.prop	a.plus.b.prop	a.div.b	GYPA.human.aa.freq
A	0.1165	0.07313	0.09264	0.08331	0.04667
R	0.02166	0.02414	0.04129	0.03369	0.04
N	0.03964	0.05007	0.06353	0.04223	0.08667
D	0.06661	0.04359	0.05876	0.05631	0.02
C	0.008991	0.02702	0.03917	0.01454	0.04
Q	0.02738	0.04395	0.03917	0.02631	0.03333
E	0.05476	0.03098	0.04553	0.05931	0.1
G	0.08051	0.107	0.09052	0.08701	0.04
H	0.04536	0.01765	0.01747	0.02469	0.07333
I	0.03719	0.04323	0.04923	0.05516	0.02
L	0.09031	0.06376	0.05823	0.07824	0.01333
K	0.1018	0.04143	0.05929	0.07408	0.06667
M	0.01962	0.005764	0.01323	0.021	0.02
F	0.05027	0.03062	0.02753	0.03646	0.04
P	0.03351	0.04575	0.03759	0.04339	0.14
S	0.04986	0.1228	0.0667	0.07547	0.1
T	0.04863	0.09114	0.06194	0.05493	0.08667
W	0.01349	0.01585	0.01588	0.01662	0.03333
Y	0.02575	0.03963	0.05717	0.03	0
V	0.06825	0.08249	0.06511	0.08724	0

Functions to calculate similarities

Two custom functions are needed: one to calculate correlates between two columns of our table, and one to calculate correlation similarities.

# Correlation used in Chou and Zhange 1992.
chou_cor <- function(x,y){
  numerator <- sum(x*y)
denominator <- sqrt((sum(x^2))*(sum(y^2)))
result <- numerator/denominator
return(result)
}

# Cosine similarity used in Higgs and Attwood (2005). 
chou_cosine <- function(z.1, z.2){
  z.1.abs <- sqrt(sum(z.1^2))
  z.2.abs <- sqrt(sum(z.2^2))
  my.cosine <- sum(z.1*z.2)/(z.1.abs*z.2.abs)
  return(my.cosine)
}

Calculate correlation between each column

corr.alpha <- chou_cor(aa.prop[,5], aa.prop[,1])
corr.beta  <- chou_cor(aa.prop[,5], aa.prop[,2])
corr.apb   <- chou_cor(aa.prop[,5], aa.prop[,3])
corr.adb   <- chou_cor(aa.prop[,5], aa.prop[,4])

Calculate cosine similarity

cos.alpha <- chou_cosine(aa.prop[,5], aa.prop[,1])
cos.beta  <- chou_cosine(aa.prop[,5], aa.prop[,2])
cos.apb   <- chou_cosine(aa.prop[,5], aa.prop[,3])
cos.adb   <- chou_cosine(aa.prop[,5], aa.prop[,4])

Calculate distance. Note: we need to flip the dataframe on its side using a command called t()

aa.prop.flipped <- t(aa.prop)
round(aa.prop.flipped,2)

##                       A    R    N    D    C    Q    E    G    H    I    L    K
## alpha.prop         0.12 0.02 0.04 0.07 0.01 0.03 0.05 0.08 0.05 0.04 0.09 0.10
## beta.prop          0.07 0.02 0.05 0.04 0.03 0.04 0.03 0.11 0.02 0.04 0.06 0.04
## a.plus.b.prop      0.09 0.04 0.06 0.06 0.04 0.04 0.05 0.09 0.02 0.05 0.06 0.06
## a.div.b            0.08 0.03 0.04 0.06 0.01 0.03 0.06 0.09 0.02 0.06 0.08 0.07
## GYPA.human.aa.freq 0.05 0.04 0.09 0.02 0.04 0.03 0.10 0.04 0.07 0.02 0.01 0.07
##                       M    F    P    S    T    W    Y    V
## alpha.prop         0.02 0.05 0.03 0.05 0.05 0.01 0.03 0.07
## beta.prop          0.01 0.03 0.05 0.12 0.09 0.02 0.04 0.08
## a.plus.b.prop      0.01 0.03 0.04 0.07 0.06 0.02 0.06 0.07
## a.div.b            0.02 0.04 0.04 0.08 0.05 0.02 0.03 0.09
## GYPA.human.aa.freq 0.02 0.04 0.14 0.10 0.09 0.03 0.00 0.00

We can get distance matrix like this

dist(aa.prop.flipped, method = "euclidean")

##                    alpha.prop  beta.prop a.plus.b.prop    a.div.b
## beta.prop          0.13342098                                    
## a.plus.b.prop      0.09281824 0.08289406                         
## a.div.b            0.06699039 0.08659174    0.06175113           
## GYPA.human.aa.freq 0.20888670 0.19389741    0.18967626 0.19197448

Individual distances using dist()

dist.alpha <- dist((aa.prop.flipped[c(1,5),]),  method = "euclidean")
dist.beta  <- dist((aa.prop.flipped[c(2,5),]),  method = "euclidean")
dist.apb   <- dist((aa.prop.flipped[c(3,5),]),  method = "euclidean")
dist.adb  <- dist((aa.prop.flipped[c(4,5),]), method = "euclidean")

Compile the information. Rounding makes it easier to read

# fold types
fold.type <- c("alpha","beta","alpha plus beta", "alpha/beta")

# data
corr.sim <- round(c(corr.alpha,corr.beta,corr.apb,corr.adb),5)
cosine.sim <- round(c(cos.alpha,cos.beta,cos.apb,cos.adb),5)
Euclidean.dist <- round(c(dist.alpha,dist.beta,dist.apb,dist.adb),5)

# summary
sim.sum <- c("","","most.sim","")
dist.sum <- c("","","min.dist","")

df <- data.frame(fold.type,
           corr.sim ,
           cosine.sim ,
           Euclidean.dist ,
           sim.sum ,
           dist.sum )

Display output

pander::pander(df)

fold.type	corr.sim	cosine.sim	Euclidean.dist	sim.sum	dist.sum
alpha	0.6972	0.6972	0.2089
beta	0.742	0.742	0.1939
alpha plus beta	0.7412	0.7412	0.1897	most.sim	min.dist
alpha/beta	0.7375	0.7375	0.192

# fold.type
# corr.sim
# cosine.sim
# Euclidean.dist
# sim.sum
# dist.sum

Percent Identity Comparisons (PID)

Data Preparation

Convert all FASTA records intro entries in a single vector. FASTA entries are contained in a list produced at the beginning of the script. They were cleaned to remove the header and newline characters.

GYPA_list

## $NP_002090.4
## [1] "MYGKIIFVLLLSEIVSISALSTTEVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ"
## 
## $XP_001145917.1
## [1] "MYGKIIFVLLLSAIVSISASSTTEVAMHTSTSSVTKSYISSETSDKHKWDIYPATPRAHEVSEIYVTTVYPPEEENGEGVQLVHRFSEPEITLIIFGVMAGVIGTILLIYYSICRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ"
## 
## $XP_003815879.2
## [1] "MYGKIIFVLLLSAIVSISASSTTEVAMHTSTSSVTKSYISSETSDKQKWDTYPATPRAHEVSEIYVTTVYPPEEENGERGQLVHRFSEPEITLIIFGVMAGVIGTILLIYYSICRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ"
## 
## $XP_021794517.1
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSSQSNDKHTSDTYPTTPSAHEVSGFSGRTHYPPEEDNRERVQLVHEFSELVIALIIFGVMAGVIGTILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT"
## 
## $XP_011842541.1
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSLQPNDKHTSDTHPTTPSAHEVSEFSGRTRYPPEEDNRERVQLVHEFSELVITLIIFGVMAGVIGTILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT"
## 
## $XP_009238620.2
## [1] "MYEKIIFVLLLSEIVSIPASNTTGEVMHTSISSSVTKSYITPQTNDKHKQDTYPATPSAHEVSEISVITIHSPEEENGERGQLVHRFSEPVITLIVFGVMAGVIGTILLISYCIRRLRKQSPSDVQPLPSPDTDVPLSSVEIENPETIDQ"
## 
## $XP_014994775.2
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLGPESYVSSQSNDKHTSDTHPTTPSAHEVSEFSGRTHYPPEEDNRERVQLVHEFSELVIALIIFGVMAGVIGTILFISYCIRRLRKKNQSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEERT"
## 
## $XP_037848439.1
## [1] "MHTSISSLGPESYVSSQSNGERVQLVHEFSELVIALIIFGVMAGVIGTILFISYGICRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQLNLFTKPNEERT"
## 
## $XP_032027955.1
## [1] "MYEKIKFVLLLLDKNKWYTYPARSVNEVSEISVTTVYPPEEENGEWRQGQLVHLFSEPVITLIIFGVMAGVIGTILSISYCIRLLRKKSPSDVQPLPSPDTEVPLSSVEIENPETIDQ"
## 
## $XP_025241177.1
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSSQSNDKHTSDTYPTTPSAHEVSGFSGRTHYPPEEDNMIALIIFGVMAGVIGTILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT"

names( GYPA_list )

##  [1] "NP_002090.4"    "XP_001145917.1" "XP_003815879.2" "XP_021794517.1"
##  [5] "XP_011842541.1" "XP_009238620.2" "XP_014994775.2" "XP_037848439.1"
##  [9] "XP_032027955.1" "XP_025241177.1"

# 10 accession numbers

length( GYPA_list )

## [1] 10

Each entry is a full entry with no spaces or parsing, and no header

GYPA_list[1]

## $NP_002090.4
## [1] "MYGKIIFVLLLSEIVSISALSTTEVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAHEVSEISVRTVYPPEEETGERVQLAHHFSEPEITLIIFGVMAGVIGTILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ"

Make each entry of the list into a vector. There are several ways to do this.

GYPA_vector <- unlist( GYPA_list )

Name the vector

names( GYPA_list )

##  [1] "NP_002090.4"    "XP_001145917.1" "XP_003815879.2" "XP_021794517.1"
##  [5] "XP_011842541.1" "XP_009238620.2" "XP_014994775.2" "XP_037848439.1"
##  [9] "XP_032027955.1" "XP_025241177.1"

names( GYPA_vector ) <- names( GYPA_list )

PID table

Do pairwise alignments for humans, chimps and 2-other species.

GYPA_human      <- GYPA_vector["NP_002090.4"]
GYPA_chimp      <- GYPA_vector["XP_001145917.1"]
GYPA_bonobo     <- GYPA_vector["XP_003815879.2"]
GYPA_orangutan  <- GYPA_vector["XP_009238620.2"]

align.human.chimp      <- Biostrings::pairwiseAlignment(GYPA_human, GYPA_chimp)
align.human.bonobo     <- Biostrings::pairwiseAlignment(GYPA_human, GYPA_bonobo)
align.human.orangutan  <- Biostrings::pairwiseAlignment(GYPA_human, GYPA_orangutan)
align.chimp.bonobo     <- Biostrings::pairwiseAlignment(GYPA_chimp, GYPA_bonobo)
align.chimp.orangutan  <- Biostrings::pairwiseAlignment(GYPA_chimp, GYPA_orangutan)
align.bonobo.orangutan <- Biostrings::pairwiseAlignment(GYPA_bonobo, GYPA_orangutan)

Build matrix

pids <- c(1,                  NA,     NA,     NA,
          Biostrings::pid(align.human.chimp),          1,     NA,     NA,
          Biostrings::pid(align.human.bonobo), Biostrings::pid(align.chimp.bonobo),      1,     NA,
          Biostrings::pid(align.human.orangutan), Biostrings::pid(align.chimp.orangutan), Biostrings::pid(align.bonobo.orangutan), 1)

mat <- matrix(pids, nrow = 4, byrow = T)
row.names(mat) <- c("Homo","Chimp","Bonobo","Orangutan")   
colnames(mat) <- c("Homo","Chimp","Bonobo","Orangutan")   

pander::pander(mat)

	Homo	Chimp	Bonobo	Orangutan
Homo	1	NA	NA	NA
Chimp	88	1	NA	NA
Bonobo	88	97.32	1	NA
Orangutan	80.67	78.67	80	1

PID methods comparison

Compare different PID methods. I did this for Humans vs. chimps

PID1 <- Biostrings::pid(align.human.chimp, type="PID1")
PID2 <- Biostrings::pid(align.human.chimp, type="PID2")
PID3 <- Biostrings::pid(align.human.chimp, type="PID3")
PID4 <- Biostrings::pid(align.human.chimp, type="PID4")

method <- c("PID1", "PID2", "PID3", "PID4")

PID <- c( PID1, PID2, PID3, PID4 )

pid.comparison <- data.frame( method, PID )

pander::pander(pid.comparison)

method	PID
PID1	88
PID2	88.59
PID3	88.59
PID4	88.29

Multiple Sequence Alignment

MSA data preparation *my MSA package was not working so I just commented this out

For use with R bioinformatics tools we need to convert our named vector to a string set using Biostrings::AAStringSet(). Note the _ss tag at the end of the object we’re assigning the output to, which designates this as a string set.

GYPA_vector_ss <- Biostrings::AAStringSet( GYPA_vector )

Building Multiple Sequence Alignment (MSA)

GYPA_align <- msa(GYPA_vector_ss, method = "ClustalW")

## use default substitution matrix

Cleaning / Setting up an MSA

msa produces a species MSA object

class( GYPA_align )

## [1] "MsaAAMultipleAlignment"
## attr(,"package")
## [1] "msa"

is( GYPA_align )

## [1] "MsaAAMultipleAlignment" "AAMultipleAlignment"    "MsaMetaData"           
## [4] "MultipleAlignment"

Default output of MSA

GYPA_align

## CLUSTAL 2.1  
## 
## Call:
##    msa(GYPA_vector_ss, method = "ClustalW")
## 
## MsaAAMultipleAlignment with 10 rows and 176 columns
##      aln                                                   names
##  [1] MYGKIIFVLLLSEIVRISASSTTVP...SSVEIENPEETDQSNLFTKPNEERT XP_014994775.2
##  [2] -------------------------...SSVEIENPEETDQLNLFTKPNEERT XP_037848439.1
##  [3] MYGKIIFVLLLSEIVRISASSTTVP...SSVEIENPEETDQSNLFTKPNEEGT XP_021794517.1
##  [4] MYGKIIFVLLLSEIVRISASSTTVP...SSVEIENPEETDQSNLFTKPNEEGT XP_025241177.1
##  [5] MYGKIIFVLLLSEIVRISASSTTVP...SSVEIENPEETDQSNLFTKPNEEGT XP_011842541.1
##  [6] MYGKIIFVLLLSAIVSISASSTTEV...SSVEIENPETSDQ------------ XP_001145917.1
##  [7] MYGKIIFVLLLSAIVSISASSTTEV...SSVEIENPETSDQ------------ XP_003815879.2
##  [8] MYGKIIFVLLLSEIVSISALSTTEV...SSVEIENPETSDQ------------ NP_002090.4
##  [9] MYEKIIFVLLLSEIVSIPASNTTGE...SSVEIENPETIDQ------------ XP_009238620.2
## [10] MYEKIKFVLLL--------------...SSVEIENPETIDQ------------ XP_032027955.1
##  Con MYGKIIFVLLLSEIV?ISASSTT??...SSVEIENPE?TDQ-?????????-? Consensus

Change class of alignment

class(GYPA_align) <- "AAMultipleAlignment"

Convert to seqinr format

GYPA_align_seqinr <- msaConvert(GYPA_align, type = "seqinr::alignment")

OPTIONAL: show output with print_msa

compbio4all::print_msa(GYPA_align_seqinr)

## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLGPESYVSSQSNDKHTSDTHPTTPSAH 0"
## [1] "--------------------------MHTSISSLGPESYVSSQSNG-------------- 0"
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSSQSNDKHTSDTYPTTPSAH 0"
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSSQSNDKHTSDTYPTTPSAH 0"
## [1] "MYGKIIFVLLLSEIVRISASSTTVPATHTSTSSLVPESYVSLQPNDKHTSDTHPTTPSAH 0"
## [1] "MYGKIIFVLLLSAIVSISASSTTEVAMHTSTSS-VTKSYISSETSDKHKWDIYPATPRAH 0"
## [1] "MYGKIIFVLLLSAIVSISASSTTEVAMHTSTSS-VTKSYISSETSDKQKWDTYPATPRAH 0"
## [1] "MYGKIIFVLLLSEIVSISALSTTEVAMHTSTSSSVTKSYISSQTNDTHKRDTYAATPRAH 0"
## [1] "MYEKIIFVLLLSEIVSIPASNTTGEVMHTSISSSVTKSYITPQTNDKHKQDTYPATPSAH 0"
## [1] "MYEKIKFVLLL---------------------------------LDKNKWYTYPAR-SVN 0"
## [1] " "
## [1] "EVSEFSGRTHYPPEEDNRER--VQLVHEFSEL------------VIALIIFGVMAGVIGT 0"
## [1] "------------------ER--VQLVHEFSEL------------VIALIIFGVMAGVIGT 0"
## [1] "EVSGFS------------GR--THYPPEEDNRERVQLVHEFSELVIALIIFGVMAGVIGT 0"
## [1] "EVSGFS------------GR--THYPPEEDN-------------MIALIIFGVMAGVIGT 0"
## [1] "EVSEFSGRTRYPPEEDNRER--VQLVHEFSEL------------VITLIIFGVMAGVIGT 0"
## [1] "EVSEIYVTTVYPPEEENGEG--VQLVHRFSEP------------EITLIIFGVMAGVIGT 0"
## [1] "EVSEIYVTTVYPPEEENGER--GQLVHRFSEP------------EITLIIFGVMAGVIGT 0"
## [1] "EVSEISVRTVYPPEEETGER--VQLAHHFSEP------------EITLIIFGVMAGVIGT 0"
## [1] "EVSEISVITIHSPEEENGER--GQLVHRFSEP------------VITLIVFGVMAGVIGT 0"
## [1] "EVSEISVTTVYPPEEENGEWRQGQLVHLFSEP------------VITLIIFGVMAGVIGT 0"
## [1] " "
## [1] "ILFISYCIRRLRKKNQSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEERT 4"
## [1] "ILFISYGICRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQLNLFTKPNEERT 4"
## [1] "ILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT 4"
## [1] "ILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT 4"
## [1] "ILFISYCIRRLRKKSPSDVQPLPPPDAEVPLSSVEIENPEETDQSNLFTKPNEEGT 4"
## [1] "ILLIYYSICRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ------------ 4"
## [1] "ILLIYYSICRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ------------ 4"
## [1] "ILLISYGIRRLIKKSPSDVKPLPSPDTDVPLSSVEIENPETSDQ------------ 4"
## [1] "ILLISYCIRRLRKQSPSDVQPLPSPDTDVPLSSVEIENPETIDQ------------ 4"
## [1] "ILSISYCIRLLRKKSPSDVQPLPSPDTEVPLSSVEIENPETIDQ------------ 4"
## [1] " "

Finished MSA

Amino Acids 105-150 seem to be sort of conserved.

class(GYPA_align) <- "AAMultipleAlignment"
ggmsa::ggmsa(GYPA_align, start = 105, end = 150)

Distance Matrix

Make a distance matrix

GYPA_dist <- seqinr::dist.alignment(GYPA_align_seqinr, 
                                       matrix = "identity")

This produces a “dist” class object

is( GYPA_dist )

## [1] "dist"     "oldClass"

class( GYPA_dist )

## [1] "dist"

Round for display

GYPA_align_seqinr_rnd <- round(GYPA_dist, 3)
GYPA_align_seqinr_rnd

##                XP_014994775.2 XP_037848439.1 XP_021794517.1 XP_025241177.1
## XP_037848439.1          0.277                                             
## XP_021794517.1          0.327          0.416                              
## XP_025241177.1          0.328          0.418          0.082               
## XP_011842541.1          0.222          0.325          0.316          0.317
## XP_001145917.1          0.537          0.514          0.567          0.562
## XP_003815879.2          0.537          0.514          0.560          0.556
## NP_002090.4             0.516          0.489          0.538          0.534
## XP_009238620.2          0.535          0.500          0.545          0.547
## XP_032027955.1          0.517          0.450          0.555          0.557
##                XP_011842541.1 XP_001145917.1 XP_003815879.2 NP_002090.4
## XP_037848439.1                                                         
## XP_021794517.1                                                         
## XP_025241177.1                                                         
## XP_011842541.1                                                         
## XP_001145917.1          0.518                                          
## XP_003815879.2          0.518          0.164                           
## NP_002090.4             0.497          0.338          0.338            
## XP_009238620.2          0.510          0.456          0.441       0.440
## XP_032027955.1          0.491          0.455          0.435       0.473
##                XP_009238620.2
## XP_037848439.1               
## XP_021794517.1               
## XP_025241177.1               
## XP_011842541.1               
## XP_001145917.1               
## XP_003815879.2               
## NP_002090.4                  
## XP_009238620.2               
## XP_032027955.1          0.405

Phylognetic trees of sequences

Build a phylogenetic tree from distance matrix

tree <- nj(GYPA_align_seqinr_rnd)

Plotting phylogenetic trees

Plot the tree

plot.phylo(tree, main="GYPA Phylogenetic Tree", 
            use.edge.length = F)

mtext(text = "GYPA Phylogenetic Tree - rooted, no branch lengths")

Jaideep Dilip Final Portfolio: GYPA

Jaideep Dilip

12/16/2021

Introduction

Resources / References

Preparation

Accession numbers

Accession number table

Data Preparation

Download Sequences

Initial Data Cleaning

General Protein Information

Protein Diagram

Draw dotplot

Protein properties compiled from databases

Protein feature prediction

Predict protein fold

Functions to calculate similarities

Percent Identity Comparisons (PID)

Data Preparation

PID table

PID methods comparison

Multiple Sequence Alignment

MSA data preparation *my MSA package was not working so I just commented this out

Building Multiple Sequence Alignment (MSA)

Cleaning / Setting up an MSA

Finished MSA

Distance Matrix

Phylognetic trees of sequences

Plotting phylogenetic trees

aa.1.1	alpha	beta	a.plus.b	a.div.b
A	285	203	175	361
R	53	67	78	146
N	97	139	120	183
D	163	121	111	244
C	22	75	74	63
Q	67	122	74	114
E	134	86	86	257
G	197	297	171	377
H	111	49	33	107
I	91	120	93	239
L	221	177	110	339
K	249	115	112	321
M	48	16	25	91
F	123	85	52	158
P	82	127	71	188
S	122	341	126	327
T	119	253	117	238
W	33	44	30	72
Y	63	110	108	130
V	167	229	123	378

aa.1.1	alpha	beta	a.plus.b	a.div.b
A	285	203	175	361
R	53	67	78	146
N	97	139	120	183
D	163	121	111	244
C	22	75	74	63
Q	67	122	74	114
E	134	86	86	257
G	197	297	171	377
H	111	49	33	107
I	91	120	93	239
L	221	177	110	339
K	249	115	112	321
M	48	16	25	91
F	123	85	52	158
P	82	127	71	188
S	122	341	126	327
T	119	253	117	238
W	33	44	30	72
Y	63	110	108	130
V	167	229	123	378

aa.1.1	alpha	beta	a.plus.b	a.div.b
A	285	203	175	361
R	53	67	78	146
N	97	139	120	183
D	163	121	111	244
C	22	75	74	63
Q	67	122	74	114
E	134	86	86	257
G	197	297	171	377
H	111	49	33	107
I	91	120	93	239
L	221	177	110	339
K	249	115	112	321
M	48	16	25	91
F	123	85	52	158
P	82	127	71	188
S	122	341	126	327
T	119	253	117	238
W	33	44	30	72
Y	63	110	108	130
V	167	229	123	378