Accessing data from Google Docs
Nathan Brouwer
11/12/2021
The goal of this exercise is to make you familiar with how to download data from Google Sheets and to briefly review some key concepts R functions and coding concepts.
We’ll do the following things
download a list of RefSeq accessions from a Google sheet remove the NAs using na.omit() select out all but one isoform using duplicated())m
Packages
## Google sheets download package
# comment this out when you are done
# install.packages("googlesheets4")
library(googlesheets4)
# comp bio packages
library(seqinr)
library(rentrez)
library(compbio4all)
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##
## Attaching package: 'IRanges'## The following object is masked from 'package:grDevices':
##
## windows## Loading required package: XVector## Loading required package: GenomeInfoDb##
## Attaching package: 'Biostrings'## The following object is masked from 'package:seqinr':
##
## translate## The following object is masked from 'package:base':
##
## strsplitDownload data
web address (URL) for the spreadsheet with the data
spreadsheet_sp <- "https://docs.google.com/spreadsheets/d/1spC_ZA3_cVuvU3e_Jfcj2nEIfzp-vaP7SA5f-qwQ1pg/edit?usp=sharing" make sure we tall the package we aren’t interested in checking user access credentials / authorization.
# be sure to run this!
googlesheets4::gs4_deauth() # <====== MUST RUN THISThird, we download our data.
NOTE!: sometimes Google Sheets or the function gets cranky and throws this error:
“Error in curl::curl_fetch_memory(url, handle = handle) : Error in the HTTP2 framing layer”
If that happens, just re-run the code.
# I include this again in case you missed is the first time : )
googlesheets4::gs4_deauth()
# download
## NOTE: if you get an error, just run the code again
refseq_column <- read_sheet(ss = spreadsheet_sp, # the url
sheet = "RefSeq_prot", # the name of the worksheet
range = "selenoprot!H1:H364",
col_names = TRUE,
na = "", # fill in empty spaces "" w/NA
trim_ws = TRUE)## v Reading from "human_gene_table".## v Range ''selenoprot'!H1:H364'.## NOTE: if you get an error, just run the code again
# for reasons we won't get into I'm going to do this
protein_refseq <- refseq_column$RefSeq_prota snapshot of the results
protein_refseq[1:10]## [1] "NP_000783.2" "NP_998758.1" "NP_001034804.1" "NP_001034805.1"
## [5] "NP_001311245.1" NA NA "NP_054644.1"
## [9] "NP_001353425.1" "NP_000784.3"downloaded another sheet of the same url
# download
## NOTE: if you get an error, just run the code again
gene_name_column <- read_sheet(ss = spreadsheet_sp, # the url
sheet = "gene", # the name of the worksheet
range = "selenoprot!A1:A364",
col_names = TRUE,
na = "", # fill in empty spaces "" w/NA
trim_ws = TRUE)## v Reading from "human_gene_table".## v Range ''selenoprot'!A1:A364'.## NOTE: if you get an error, just run the code again
# for reasons we won't get into I'm going to do this
gene <- gene_name_column$geneCheck data
We downloaded a vector of character data:
is(protein_refseq)## [1] "character" "vector"
## [3] "data.frameRowLabels" "SuperClassMethod"
## [5] "character_OR_connection" "character_OR_NULL"
## [7] "atomic" "EnumerationValue"
## [9] "vector_OR_Vector" "vector_OR_factor"class(protein_refseq)## [1] "character"length(protein_refseq)## [1] 363protein_refseq[1:10]## [1] "NP_000783.2" "NP_998758.1" "NP_001034804.1" "NP_001034805.1"
## [5] "NP_001311245.1" NA NA "NP_054644.1"
## [9] "NP_001353425.1" "NP_000784.3"check for the presence of NAs using is.na
is.na(protein_refseq)## [1] FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE TRUE
## [13] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE
## [37] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## [73] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE
## [253] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
## [289] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
## [337] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [349] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSEWe can count them up a couple different ways. First, with table(), were “false” means “not an NA”
table(is.na(protein_refseq))##
## FALSE TRUE
## 334 29another way to count NAs
# make a vector TRUE/FALSE
temp <- is.na(protein_refseq)
# this subsets by the TRUE statements
protein_refseq[temp]## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [26] NA NA NA NAtemp2 <- protein_refseq[temp]
# count them up
length(temp2)## [1] 29Make a dataframe
make a dataframe
seleno_df <- data.frame(gene = gene,
protein_refseq = protein_refseq)check df
summary(seleno_df)## gene protein_refseq
## Length:363 Length:363
## Class :character Class :character
## Mode :character Mode :characterhead(seleno_df)## gene protein_refseq
## 1 DIO1 NP_000783.2
## 2 DIO1 NP_998758.1
## 3 DIO1 NP_001034804.1
## 4 DIO1 NP_001034805.1
## 5 DIO1 NP_001311245.1
## 6 DIO1 <NA>Remove NAs
remove NAs
# omit NAs
seleno_df_noNA <- na.omit(seleno_df)
# check length- should be shorter
dim(seleno_df)## [1] 363 2dim(seleno_df_noNA)## [1] 334 2Select just 1 isoform
The same gene can appear multiple times because multiple isoforms are listed.
head(seleno_df_noNA)## gene protein_refseq
## 1 DIO1 NP_000783.2
## 2 DIO1 NP_998758.1
## 3 DIO1 NP_001034804.1
## 4 DIO1 NP_001034805.1
## 5 DIO1 NP_001311245.1
## 8 DIO2 NP_054644.1We can use the function unique() to select just one row for each gene.
genes_unique <- unique(seleno_df_noNA$gene)
length(genes_unique)## [1] 37genes_unique## [1] "DIO1" "DIO2" "DIO3" "GPX1" "GPX2" "GPX3"
## [7] "GPX4" "GPX6" "MSRB1" "SELENOF" "SELENOH" "SELENOI"
## [13] "SELENOK" "SELENOM" "SELENON" "SELENOO" "SELENOP" "SELENOS"
## [19] "SELENOT" "SELENOV" "SELENOW" "SEPHS2" "TXNRD1" "TXNRD2"
## [25] "TXNRD3" "SELENOP1" "SELENOP2" "SELENOU" "SELENOW1" "SELENOW2"
## [31] "SELENOE" "SELENOJ" "SELENOL" "SELENOO1" "SELENOO2" "SELENOT1"
## [37] "SELENOT2"unique() just gives us the unique elements. A related function, duplicated(), gives us the location of duplicated elements in the vector. FALSE means “not duplicated yet” or “first instance so far”.
i.dups <- duplicated(seleno_df_noNA$gene)We can remove the duplicates using a form of reverse indexing where the “!” means “not”. (You don’t need to know this for the exam)
seleno_df_noNA[!i.dups, ]## gene protein_refseq
## 1 DIO1 NP_000783.2
## 8 DIO2 NP_054644.1
## 14 DIO3 NP_001353.4
## 15 GPX1 NP_000572.2
## 20 GPX2 NP_002074.2
## 24 GPX3 NP_002075.2
## 26 GPX4 NP_002076.2
## 29 GPX6 NP_874360.1
## 30 MSRB1 NP_057416.1
## 31 SELENOF NP_004252.2
## 35 SELENOH NP_734467.1
## 37 SELENOI NP_277040.1
## 39 SELENOK NP_067060.2
## 40 SELENOM NP_536355.1
## 41 SELENON NP_996809.1
## 43 SELENOO NP_113642.1
## 44 SELENOP NP_005401.3
## 47 SELENOS NP_060915.2
## 49 SELENOT NP_057359.2
## 50 SELENOV NP_874363.1
## 53 SELENOW NP_003000.1
## 54 SEPHS2 NP_036380.2
## 55 TXNRD1 NP_877393.1
## 62 TXNRD2 NP_006431.2
## 69 TXNRD3 NP_443115.1
## 232 SELENOP1 NP_001026780.2
## 233 SELENOP2 NP_001335698.1
## 236 SELENOU NP_001180447.1
## 268 SELENOW1 NP_001291715.2
## 269 SELENOW2 NP_001341647.1
## 334 SELENOE NP_001182713.2
## 338 SELENOJ NP_001180398.1
## 340 SELENOL NP_001177311.1
## 343 SELENOO1 NP_001038336.2
## 344 SELENOO2 NP_001335014.1
## 348 SELENOT1 NP_840075.2
## 350 SELENOT2 NP_001091957.2Make a dataframe of non-duplicated genes
seleno_df_noDups <- seleno_df_noNA[!i.dups, ]
dim(seleno_df_noDups)## [1] 37 2Select 2 random sequences
Let’s select 2 random sequences to work with. We’ll use sample() to select a random index number to get
First, lets make a vector that contains a unique number for each row of data
indices <- 1:nrow(seleno_df_noDups)This would do the same thing
# with dim
indices <- 1:dim(seleno_df_noDups)[1]
# with length
indices <- 1:length(seleno_df_noDups$gene)or hard-coded
indices <- 1:37We can then use sample() to select 2 random numbers from this vector.
For x = we’ll use our vector of indices (1 to 37). For size we’ll use 2, since we want to pull out just 2 numbers. For replace we’ll use FALSE since we don’t want to be ale to select the same number twice.
i.random.genes <- sample(x = indices,
size = 2,
replace = FALSE)Hard coded this would be
i.random.genes <- sample(x = c(1:37),
size = 2,
replace = FALSE)This gives me 2 indices values.
i.random.genes## [1] 31 28I can now use these index values to pull out 2 rows of data
seleno_df_noNA[i.random.genes, ]## gene protein_refseq
## 41 SELENON NP_996809.1
## 37 SELENOI NP_277040.1Hard coded, this would be something like this for whichever genes happen to have been selected
seleno_df_noNA[c(8,7), ]## gene protein_refseq
## 10 DIO2 NP_000784.3
## 9 DIO2 NP_001353425.1Downloading genes
I will now download the two FASTA files.
rentrez::entrez_fetch(id = "NP_000784.3",
db = "protein",
rettype = "fasta")## [1] ">NP_000784.3 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"rentrez::entrez_fetch(id = "NP_001353425.1",
db = "protein",
rettype = "fasta")## [1] ">NP_001353425.1 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"save this into vectors
prot1 <- rentrez::entrez_fetch(id = "NP_000784.3",
db = "protein",
rettype = "fasta")
prot2 <- rentrez::entrez_fetch(id = "NP_001353425.1",
db = "protein",
rettype = "fasta")I can put them into a list like this
# make the list
seleno_thingy <- vector("list", 1)
# add the first fasta
seleno_thingy[[1]] <- prot1
# See the result
seleno_thingy## [[1]]
## [1] ">NP_000784.3 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"# add the first fasta
seleno_thingy[[2]] <- prot2
# see the result
seleno_thingy## [[1]]
## [1] ">NP_000784.3 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"
##
## [[2]]
## [1] ">NP_001353425.1 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"# WHAT DOES THIS DO?
names(seleno_thingy) <- c("prot1", "prot2")
#Output
seleno_thingy## $prot1
## [1] ">NP_000784.3 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"
##
## $prot2
## [1] ">NP_001353425.1 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"Elements of the list are accessed like this
seleno_thingy[[1]]## [1] ">NP_000784.3 type II iodothyronine deiodinase isoform b [Homo sapiens]\nMGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLD\nAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQL\nPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRV\nVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR\n\n"I’ll clean them with fasta_cleaner()
# first, make a copy of the list for storing the clean data
## I'm just going to copy over the old data
seleno_thingy_clean <- seleno_thingy
# HOW TO MAKE THIS MORE COMPACT?
for(i in 1:length(seleno_thingy_clean)){
clean_fasta_temp <- compbio4all::fasta_cleaner(seleno_thingy[[i]],
parse = T)
seleno_thingy_clean[[i]] <- clean_fasta_temp
}Now the data looks like this Each element of the list is a vector of character data
seleno_thingy_clean## $prot1
## [1] "M" "G" "I" "L" "S" "V" "D" "L" "L" "I" "T" "L" "Q" "I" "L" "P" "V" "F"
## [19] "F" "S" "N" "C" "L" "F" "L" "A" "L" "Y" "D" "S" "V" "I" "L" "L" "K" "H"
## [37] "V" "V" "L" "L" "L" "S" "R" "S" "K" "S" "T" "R" "G" "E" "W" "R" "R" "M"
## [55] "L" "T" "S" "E" "G" "L" "R" "C" "V" "W" "K" "S" "F" "L" "L" "D" "A" "Y"
## [73] "K" "Q" "V" "K" "L" "G" "E" "D" "A" "P" "N" "S" "S" "V" "V" "H" "V" "S"
## [91] "S" "T" "E" "G" "G" "D" "N" "S" "G" "N" "G" "T" "Q" "E" "K" "I" "A" "E"
## [109] "G" "A" "T" "C" "H" "L" "L" "D" "F" "A" "S" "P" "E" "R" "P" "L" "V" "V"
## [127] "N" "F" "G" "S" "A" "T" "U" "P" "P" "F" "T" "S" "Q" "L" "P" "A" "F" "R"
## [145] "K" "L" "V" "E" "E" "F" "S" "S" "V" "A" "D" "F" "L" "L" "V" "Y" "I" "D"
## [163] "E" "A" "H" "P" "S" "D" "G" "W" "A" "I" "P" "G" "D" "S" "S" "L" "S" "F"
## [181] "E" "V" "K" "K" "H" "Q" "N" "Q" "E" "D" "R" "C" "A" "A" "A" "Q" "Q" "L"
## [199] "L" "E" "R" "F" "S" "L" "P" "P" "Q" "C" "R" "V" "V" "A" "D" "R" "M" "D"
## [217] "N" "N" "A" "N" "I" "A" "Y" "G" "V" "A" "F" "E" "R" "V" "C" "I" "V" "Q"
## [235] "R" "Q" "K" "I" "A" "Y" "L" "G" "G" "K" "G" "P" "F" "S" "Y" "N" "L" "Q"
## [253] "E" "V" "R" "H" "W" "L" "E" "K" "N" "F" "S" "K" "R"
##
## $prot2
## [1] "M" "G" "I" "L" "S" "V" "D" "L" "L" "I" "T" "L" "Q" "I" "L" "P" "V" "F"
## [19] "F" "S" "N" "C" "L" "F" "L" "A" "L" "Y" "D" "S" "V" "I" "L" "L" "K" "H"
## [37] "V" "V" "L" "L" "L" "S" "R" "S" "K" "S" "T" "R" "G" "E" "W" "R" "R" "M"
## [55] "L" "T" "S" "E" "G" "L" "R" "C" "V" "W" "K" "S" "F" "L" "L" "D" "A" "Y"
## [73] "K" "Q" "V" "K" "L" "G" "E" "D" "A" "P" "N" "S" "S" "V" "V" "H" "V" "S"
## [91] "S" "T" "E" "G" "G" "D" "N" "S" "G" "N" "G" "T" "Q" "E" "K" "I" "A" "E"
## [109] "G" "A" "T" "C" "H" "L" "L" "D" "F" "A" "S" "P" "E" "R" "P" "L" "V" "V"
## [127] "N" "F" "G" "S" "A" "T" "U" "P" "P" "F" "T" "S" "Q" "L" "P" "A" "F" "R"
## [145] "K" "L" "V" "E" "E" "F" "S" "S" "V" "A" "D" "F" "L" "L" "V" "Y" "I" "D"
## [163] "E" "A" "H" "P" "S" "D" "G" "W" "A" "I" "P" "G" "D" "S" "S" "L" "S" "F"
## [181] "E" "V" "K" "K" "H" "Q" "N" "Q" "E" "D" "R" "C" "A" "A" "A" "Q" "Q" "L"
## [199] "L" "E" "R" "F" "S" "L" "P" "P" "Q" "C" "R" "V" "V" "A" "D" "R" "M" "D"
## [217] "N" "N" "A" "N" "I" "A" "Y" "G" "V" "A" "F" "E" "R" "V" "C" "I" "V" "Q"
## [235] "R" "Q" "K" "I" "A" "Y" "L" "G" "G" "K" "G" "P" "F" "S" "Y" "N" "L" "Q"
## [253] "E" "V" "R" "H" "W" "L" "E" "K" "N" "F" "S" "K" "R"class(seleno_thingy_clean[[1]])## [1] "character"is(seleno_thingy_clean[[1]])## [1] "character" "vector"
## [3] "data.frameRowLabels" "SuperClassMethod"
## [5] "character_OR_connection" "character_OR_NULL"
## [7] "atomic" "EnumerationValue"
## [9] "vector_OR_Vector" "vector_OR_factor"is.vector(seleno_thingy_clean[[1]])## [1] TRUEMake an dotplot
For old-times sake we can make a dotplot.
Now for a dotplot
We’ll take each vector in the list and make it into a separate object
prot1_vector <- seleno_thingy_clean[[1]]
prot2_vector <- seleno_thingy_clean[[2]]We can dotplot like this
seqinr::dotPlot(prot1_vector,
prot1_vector)
WHAT DID I DO DIFFERENTLY HERE?
seqinr::dotPlot(seleno_thingy_clean[[1]],
seleno_thingy_clean[[2]])Pairwise alignment
dotPlot likes things in a single vector, but pairwiseAlignment like a single string of characters, so as always we have to process the data.
We’ll use the paste() command to take everything in the vector and squish it together “WHAT DOES”" MEAN?"" means no seperator
prot1_str <- paste(seleno_thingy_clean[[1]],sep = "", collapse = "")
prot2_str <- paste(seleno_thingy_clean[[2]],sep = "", collapse = "")So now things look like this a vector of length 1 with the entire sequence in a single element of the vector.
prot1_str## [1] "MGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGLRCVWKSFLLDAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASPERPLVVNFGSATUPPFTSQLPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSFEVKKHQNQEDRCAAAQQLLERFSLPPQCRVVADRMDNNANIAYGVAFERVCIVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR"Protein alignments need a amino acid transition matrix, and we need to use data() to bring those up into active memory (VERY IMPORTANT STEP!)
data(BLOSUM50)The alignment
align_out <- Biostrings::pairwiseAlignment(pattern = prot1_str,
subject = prot2_str,
type = "global",
gapOpening = -9.5,
gapExtension = -0.5)The raw output
align_out## Global PairwiseAlignmentsSingleSubject (1 of 1)
## pattern: MGILSVDLLITLQILPVFFSNCLFLALYDSVILL...IVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR
## subject: MGILSVDLLITLQILPVFFSNCLFLALYDSVILL...IVQRQKIAYLGGKGPFSYNLQEVRHWLEKNFSKR
## score: 1159.127The full alignment looks like this. It should have LOTS of gaps because these are two random sequences.
compbio4all::print_pairwise_alignment(align_out)## [1] "MGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGL 60"
## [1] "MGILSVDLLITLQILPVFFSNCLFLALYDSVILLKHVVLLLSRSKSTRGEWRRMLTSEGL 60"
## [1] " "
## [1] "RCVWKSFLLDAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASP 120"
## [1] "RCVWKSFLLDAYKQVKLGEDAPNSSVVHVSSTEGGDNSGNGTQEKIAEGATCHLLDFASP 120"
## [1] " "
## [1] "ERPLVVNFGSATUPPFTSQLPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSF 180"
## [1] "ERPLVVNFGSATUPPFTSQLPAFRKLVEEFSSVADFLLVYIDEAHPSDGWAIPGDSSLSF 180"
## [1] " "
## [1] "EVKKHQNQEDRCAAAQQLLERFSLPPQCRVVADRMDNNANIAYGVAFERVCIVQRQKIAY 240"
## [1] "EVKKHQNQEDRCAAAQQLLERFSLPPQCRVVADRMDNNANIAYGVAFERVCIVQRQKIAY 240"
## [1] " "
## [1] "LGGKGPFSYNLQEVRHWLEKNFSKR 300"
## [1] "LGGKGPFSYNLQEVRHWLEKNFSKR 300"
## [1] " "These are two randomly chosen sequences, so the alignment should be pretty …WHAT?
The score is negative, but on its own that is hard to interpret. Even good alignments can be negative.
score(align_out)## [1] 1159.127pid gives us percent identity
pid(align_out)## [1] 100Of course, pid can be calculated several ways (WHY IS THIS AN ISSUE / POSSIBLE?)
pid(align_out,type = "PID1")## [1] 100pid(align_out,type = "PID2")## [1] 100pid(align_out,type = "PID3")## [1] 100pid(align_out,type = "PID4")## [1] 100