Assignment: Your assignment is to use your notes from class - along with help from classmates, UTAs, and me - to turn this script into a fleshed-out description of what is going on.
This is a substantial project - we’ll work on it in steps over the rest of the unit.
We are currently focused on the overall process and will cover the details over the rest of this unit.
Your first assignment is to get this script to run from top to bottom by adding all of the missing R commands. Once you have done that, you can knit it into an HTML file and upload it to RPubs. (Note - you’ll need to add the YAML header!)
Your second assignment, which will be posted later, is to answer all the TODO and other prompts to add information. You can start on this, but you don’t have to do this on your first time through the code.
Delete all the prompts like TODO() as you compete them. Use RStudio’s search function to see if you’ve missed any - there are a LOT!
Add YAML header!!! Give it a title
title: “MSA Walkthrough assignment” author: “Thulasi Varatharajan” date: “9/26/21” output: html_document
A complete bioinformatics workflow in R
By: Nathan L. Brouwer
“Worked example: Building a phylogeny in R”
Introduction
Describe how phylogeneies can be used in biology (readings will be assigned)
Phylogenies can be used to find relationships between different groups of organisms. Phylogenies can show where groups of organisms started to relate in the evolutionary lines, and also where related species started to diverge in evolutionary traits.
Vocab
Make a list of at least 10 vocab terms that are important (don’t have to define)
phylogenetic trees multiple sequence alignment pairwise sequence alignment accession number fasta file distance matrix shroom reproducible work flow nucleotide amino acid
Key functions
Make a list of at least 5 key functions Put in the format of package::function
Biostrings::PairWiseAlignment() stringr::str_split() rentrez::entrez_fetch() Biostrings::pid() BioStrings::AAStringSet()
Software Preliminaires
Add the necessary calls to library() to load call packages Indicate which packages cam from Bioconducotr, CRAN, and GitHub
Load packages into memory
# github packages
library(compbio4all)
# CRAN packages
library(rentrez)
library(seqinr)
library(ape)
# Bioconductor packages
library(msa)
library(Biostrings)Aligning macromolecular sequences
TODO: Fill in the XXXXXs and write a 1-2 sentence of what is going on here.
Add the package that is where entrez_fetch is from using :: notation
Pairwise alignment is happening here- it compares two sequences to find the best alignment between both.
align.h3.vs.m3a <- Biostrings::pairwiseAlignment (hShroom3, mShroom3a)
# Human shroom 3 (H. sapiens)
hShroom3 <- entrez_fetch(db = "protein",
id = "NP_065910",
rettype = "fasta")TODO:explain what cat() is doing
Takes contents of hShroom3 and converts it to characters and then prints it out- also called concatenate and print.
cat(hShroom3)## >NP_065910.3 protein Shroom3 [Homo sapiens]
## MMRTTEDFHKPSATLNSNTATKGRYIYLEAFLEGGAPWGFTLKGGLEHGEPLIISKVEEGGKADTLSSKL
## QAGDEVVHINEVTLSSSRKEAVSLVKGSYKTLRLVVRRDVCTDPGHADTGASNFVSPEHLTSGPQHRKAA
## WSGGVKLRLKHRRSEPAGRPHSWHTTKSGEKQPDASMMQISQGMIGPPWHQSYHSSSSTSDLSNYDHAYL
## RRSPDQCSSQGSMESLEPSGAYPPCHLSPAKSTGSIDQLSHFHNKRDSAYSSFSTSSSILEYPHPGISGR
## ERSGSMDNTSARGGLLEGMRQADIRYVKTVYDTRRGVSAEYEVNSSALLLQGREARASANGQGYDKWSNI
## PRGKGVPPPSWSQQCPSSLETATDNLPPKVGAPLPPARSDSYAAFRHRERPSSWSSLDQKRLCRPQANSL
## GSLKSPFIEEQLHTVLEKSPENSPPVKPKHNYTQKAQPGQPLLPTSIYPVPSLEPHFAQVPQPSVSSNGM
## LYPALAKESGYIAPQGACNKMATIDENGNQNGSGRPGFAFCQPLEHDLLSPVEKKPEATAKYVPSKVHFC
## SVPENEEDASLKRHLTPPQGNSPHSNERKSTHSNKPSSHPHSLKCPQAQAWQAGEDKRSSRLSEPWEGDF
## QEDHNANLWRRLEREGLGQSLSGNFGKTKSAFSSLQNIPESLRRHSSLELGRGTQEGYPGGRPTCAVNTK
## AEDPGRKAAPDLGSHLDRQVSYPRPEGRTGASASFNSTDPSPEEPPAPSHPHTSSLGRRGPGPGSASALQ
## GFQYGKPHCSVLEKVSKFEQREQGSQRPSVGGSGFGHNYRPHRTVSTSSTSGNDFEETKAHIRFSESAEP
## LGNGEQHFKNGELKLEEASRQPCGQQLSGGASDSGRGPQRPDARLLRSQSTFQLSSEPEREPEWRDRPGS
## PESPLLDAPFSRAYRNSIKDAQSRVLGATSFRRRDLELGAPVASRSWRPRPSSAHVGLRSPEASASASPH
## TPRERHSVTPAEGDLARPVPPAARRGARRRLTPEQKKRSYSEPEKMNEVGIVEEAEPAPLGPQRNGMRFP
## ESSVADRRRLFERDGKACSTLSLSGPELKQFQQSALADYIQRKTGKRPTSAAGCSLQEPGPLRERAQSAY
## LQPGPAALEGSGLASASSLSSLREPSLQPRREATLLPATVAETQQAPRDRSSSFAGGRRLGERRRGDLLS
## GANGGTRGTQRGDETPREPSSWGARAGKSMSAEDLLERSDVLAGPVHVRSRSSPATADKRQDVLLGQDSG
## FGLVKDPCYLAGPGSRSLSCSERGQEEMLPLFHHLTPRWGGSGCKAIGDSSVPSECPGTLDHQRQASRTP
## CPRPPLAGTQGLVTDTRAAPLTPIGTPLPSAIPSGYCSQDGQTGRQPLPPYTPAMMHRSNGHTLTQPPGP
## RGCEGDGPEHGVEEGTRKRVSLPQWPPPSRAKWAHAAREDSLPEESSAPDFANLKHYQKQQSLPSLCSTS
## DPDTPLGAPSTPGRISLRISESVLRDSPPPHEDYEDEVFVRDPHPKATSSPTFEPLPPPPPPPPSQETPV
## YSMDDFPPPPPHTVCEAQLDSEDPEGPRPSFNKLSKVTIARERHMPGAAHVVGSQTLASRLQTSIKGSEA
## ESTPPSFMSVHAQLAGSLGGQPAPIQTQSLSHDPVSGTQGLEKKVSPDPQKSSEDIRTEALAKEIVHQDK
## SLADILDPDSRLKTTMDLMEGLFPRDVNLLKENSVKRKAIQRTVSSSGCEGKRNEDKEAVSMLVNCPAYY
## SVSAPKAELLNKIKEMPAEVNEEEEQADVNEKKAELIGSLTHKLETLQEAKGSLLTDIKLNNALGEEVEA
## LISELCKPNEFDKYRMFIGDLDKVVNLLLSLSGRLARVENVLSGLGEDASNEERSSLYEKRKILAGQHED
## ARELKENLDRRERVVLGILANYLSEEQLQDYQHFVKMKSTLLIEQRKLDDKIKLGQEQVKCLLESLPSDF
## IPKAGALALPPNLTSEPIPAGGCTFSGIFPTLTSPLTODO: explain what this code chunk is doing
This code chunk is assigning each fasta data set that entrez_fetch gets and assigns it to each appropriate variable.
# Mouse shroom 3a (M. musculus)
mShroom3a <- entrez_fetch(db = "protein",
id = "AAF13269",
rettype = "fasta")
# Human shroom 2 (H. sapiens)
hShroom2 <- entrez_fetch(db = "protein",
id = "CAA58534",
rettype = "fasta")
# Sea-urchin shroom
sShroom <- entrez_fetch(db = "protein",
id = "XP_783573",
rettype = "fasta")TODO: Explain what this code chunk is doing
This code chunk is taking every data set and finding the number of amino acids(characters) in each data set with the nchar function.
nchar(hShroom3)## [1] 2070nchar(mShroom3a)## [1] 2083nchar(sShroom)## [1] 1758nchar(hShroom2)## [1] 1673Prepping macromolecular sequences
TODO: Explain what this function does
This function is cleaning the fasta files that are being obtained by rentrez_fetch
fasta_cleaner## function (fasta_object, parse = TRUE)
## {
## fasta_object <- sub("^(>)(.*?)(\\n)(.*)(\\n\\n)", "\\4",
## fasta_object)
## fasta_object <- gsub("\n", "", fasta_object)
## if (parse == TRUE) {
## fasta_object <- stringr::str_split(fasta_object, pattern = "",
## simplify = FALSE)
## }
## return(fasta_object[[1]])
## }
## <bytecode: 0x00000000230cc458>
## <environment: namespace:compbio4all>TODO: explain how to add the function to your R session even if you can’t download compbio4all
You can create a new variable called fasta_cleaner and then assign fasta_cleaner to the necessary objects to the variable so that it does essentially the same thing
fasta_cleaner <- function(fasta_object, parse = TRUE){
fasta_object <- sub("^(>)(.*?)(\\n)(.*)(\\n\\n)","\\4",fasta_object)
fasta_object <- gsub("\n", "", fasta_object)
if(parse == TRUE){
fasta_object <- stringr::str_split(fasta_object,
pattern = "",
simplify = FALSE)
}
return(fasta_object[[1]])
}TODO: briefly explain what this code chunk is doing
This code chunk takes the fasta files and cleans them for the alignments to be done after
hShroom3 <- fasta_cleaner(hShroom3, parse = F)
mShroom3a <- fasta_cleaner(mShroom3a, parse = F)
hShroom2 <- fasta_cleaner(hShroom2, parse = F)
sShroom <- fasta_cleaner(sShroom, parse = F)hShroom3## [1] "MMRTTEDFHKPSATLNSNTATKGRYIYLEAFLEGGAPWGFTLKGGLEHGEPLIISKVEEGGKADTLSSKLQAGDEVVHINEVTLSSSRKEAVSLVKGSYKTLRLVVRRDVCTDPGHADTGASNFVSPEHLTSGPQHRKAAWSGGVKLRLKHRRSEPAGRPHSWHTTKSGEKQPDASMMQISQGMIGPPWHQSYHSSSSTSDLSNYDHAYLRRSPDQCSSQGSMESLEPSGAYPPCHLSPAKSTGSIDQLSHFHNKRDSAYSSFSTSSSILEYPHPGISGRERSGSMDNTSARGGLLEGMRQADIRYVKTVYDTRRGVSAEYEVNSSALLLQGREARASANGQGYDKWSNIPRGKGVPPPSWSQQCPSSLETATDNLPPKVGAPLPPARSDSYAAFRHRERPSSWSSLDQKRLCRPQANSLGSLKSPFIEEQLHTVLEKSPENSPPVKPKHNYTQKAQPGQPLLPTSIYPVPSLEPHFAQVPQPSVSSNGMLYPALAKESGYIAPQGACNKMATIDENGNQNGSGRPGFAFCQPLEHDLLSPVEKKPEATAKYVPSKVHFCSVPENEEDASLKRHLTPPQGNSPHSNERKSTHSNKPSSHPHSLKCPQAQAWQAGEDKRSSRLSEPWEGDFQEDHNANLWRRLEREGLGQSLSGNFGKTKSAFSSLQNIPESLRRHSSLELGRGTQEGYPGGRPTCAVNTKAEDPGRKAAPDLGSHLDRQVSYPRPEGRTGASASFNSTDPSPEEPPAPSHPHTSSLGRRGPGPGSASALQGFQYGKPHCSVLEKVSKFEQREQGSQRPSVGGSGFGHNYRPHRTVSTSSTSGNDFEETKAHIRFSESAEPLGNGEQHFKNGELKLEEASRQPCGQQLSGGASDSGRGPQRPDARLLRSQSTFQLSSEPEREPEWRDRPGSPESPLLDAPFSRAYRNSIKDAQSRVLGATSFRRRDLELGAPVASRSWRPRPSSAHVGLRSPEASASASPHTPRERHSVTPAEGDLARPVPPAARRGARRRLTPEQKKRSYSEPEKMNEVGIVEEAEPAPLGPQRNGMRFPESSVADRRRLFERDGKACSTLSLSGPELKQFQQSALADYIQRKTGKRPTSAAGCSLQEPGPLRERAQSAYLQPGPAALEGSGLASASSLSSLREPSLQPRREATLLPATVAETQQAPRDRSSSFAGGRRLGERRRGDLLSGANGGTRGTQRGDETPREPSSWGARAGKSMSAEDLLERSDVLAGPVHVRSRSSPATADKRQDVLLGQDSGFGLVKDPCYLAGPGSRSLSCSERGQEEMLPLFHHLTPRWGGSGCKAIGDSSVPSECPGTLDHQRQASRTPCPRPPLAGTQGLVTDTRAAPLTPIGTPLPSAIPSGYCSQDGQTGRQPLPPYTPAMMHRSNGHTLTQPPGPRGCEGDGPEHGVEEGTRKRVSLPQWPPPSRAKWAHAAREDSLPEESSAPDFANLKHYQKQQSLPSLCSTSDPDTPLGAPSTPGRISLRISESVLRDSPPPHEDYEDEVFVRDPHPKATSSPTFEPLPPPPPPPPSQETPVYSMDDFPPPPPHTVCEAQLDSEDPEGPRPSFNKLSKVTIARERHMPGAAHVVGSQTLASRLQTSIKGSEAESTPPSFMSVHAQLAGSLGGQPAPIQTQSLSHDPVSGTQGLEKKVSPDPQKSSEDIRTEALAKEIVHQDKSLADILDPDSRLKTTMDLMEGLFPRDVNLLKENSVKRKAIQRTVSSSGCEGKRNEDKEAVSMLVNCPAYYSVSAPKAELLNKIKEMPAEVNEEEEQADVNEKKAELIGSLTHKLETLQEAKGSLLTDIKLNNALGEEVEALISELCKPNEFDKYRMFIGDLDKVVNLLLSLSGRLARVENVLSGLGEDASNEERSSLYEKRKILAGQHEDARELKENLDRRERVVLGILANYLSEEQLQDYQHFVKMKSTLLIEQRKLDDKIKLGQEQVKCLLESLPSDFIPKAGALALPPNLTSEPIPAGGCTFSGIFPTLTSPL"Aligning sequences
TODO: give this a title. Explain what code below is doing
The code below is lining up the sequences of human Shroom3 and mice.
# Pairwise Alignment
align.h3.vs.m3a <- Biostrings::pairwiseAlignment(
hShroom3,
mShroom3a)TODO: In 1-2 sentence explain what this object shows
This object shows the aligned human and mice DNA and points out the most similar regions and the number of base pairs in common
align.h3.vs.m3a## Global PairwiseAlignmentsSingleSubject (1 of 1)
## pattern: MMRTTEDFHKPSATLN-SNTATKGRYIYLEAFLE...KAGALALPPNLTSEPIPAGGCTFSGIFPTLTSPL
## subject: MK-TPENLEEPSATPNPSRTPTE-RFVYLEALLE...KAGAISLPPALTGHATPGGTSVFGGVFPTLTSPL
## score: 2189.934TODO: explain what this is showing
This function is showing percent identity in the aligned sequence between aligned human shroom 3 and mice DNA
# add necessary function
Biostrings::pid(align.h3.vs.m3a)## [1] 70.56511TODO: briefly explain what is going on here versus the previous code chunk
In the previous code chunk, the human shroom and mice shroom sequences were getting aligned. However, in this chunk, the human shroom and another human shroom (human shroom 2) are getting aligned and compared(human hemoglobin and myoglobin sequences)
align.h3.vs.h2 <- Biostrings::pairwiseAlignment(
hShroom3,
hShroom2)TODO: explain what is going on here and compare and contrast with previous ouput
This chunk is giving the similarity score value between human myoglobin and human hemoglobin compared to the mice and human hemoglobin, this has more of a similarity
score(align.h3.vs.h2)## [1] -5673.853TODO: briefly explian the difference between the output of score() and pid() (can be very brief - we’ll get into the details later)
score() is a generalized similarity index for genomic sequences, while pid calculates sequence identity for specifically pairwise sequence alignment.
Biostrings::pid(align.h3.vs.h2)## [1] 33.83277The shroom family of genes
TODO: briefly explain why I have this whole table here
This whole table is a series of the best sequenced shroom sequences and their accession numbers, and they are being added to the database of already existing shroom sequences.
shroom_table <- c("CAA78718" , "X. laevis Apx" , "xShroom1",
"NP_597713" , "H. sapiens APXL2" , "hShroom1",
"CAA58534" , "H. sapiens APXL", "hShroom2",
"ABD19518" , "M. musculus Apxl" , "mShroom2",
"AAF13269" , "M. musculus ShroomL" , "mShroom3a",
"AAF13270" , "M. musculus ShroomS" , "mShroom3b",
"NP_065910", "H. sapiens Shroom" , "hShroom3",
"ABD59319" , "X. laevis Shroom-like", "xShroom3",
"NP_065768", "H. sapiens KIAA1202" , "hShroom4a",
"AAK95579" , "H. sapiens SHAP-A" , "hShroom4b",
#"DQ435686" , "M. musculus KIAA1202" , "mShroom4",
"ABA81834" , "D. melanogaster Shroom", "dmShroom",
"EAA12598" , "A. gambiae Shroom", "agShroom",
"XP_392427" , "A. mellifera Shroom" , "amShroom",
"XP_783573" , "S. purpuratus Shroom" , "spShroom") #sea urchinTODO: write a short sentence explaining what this next code chunk will do, then annotate each line with what was done.
This code chunk is taking the really long vector from the concatenation before and making it into a data table, which is a little more organized to look at.
# convert to matrix
shroom_table_matrix <- matrix(shroom_table,
byrow = T,
nrow = 14)
# convert to data frame
shroom_table <- data.frame(shroom_table_matrix,
stringsAsFactors = F)
# Three columns
names(shroom_table) <- c("accession", "name.orig","name.new")
# Create simplified species names
shroom_table$spp <- "Homo"
shroom_table$spp[grep("laevis",shroom_table$name.orig)] <- "Xenopus"
shroom_table$spp[grep("musculus",shroom_table$name.orig)] <- "Mus"
shroom_table$spp[grep("melanogaster",shroom_table$name.orig)] <- "Drosophila"
shroom_table$spp[grep("gambiae",shroom_table$name.orig)] <- "mosquito"
shroom_table$spp[grep("mellifera",shroom_table$name.orig)] <- "bee"
shroom_table$spp[grep("purpuratus",shroom_table$name.orig)] <- "sea urchin"TODO: in a brief sentence explain what this is doing
This is just accessing what is inside the vector shroom_table.
shroom_table## accession name.orig name.new spp
## 1 CAA78718 X. laevis Apx xShroom1 Xenopus
## 2 NP_597713 H. sapiens APXL2 hShroom1 Homo
## 3 CAA58534 H. sapiens APXL hShroom2 Homo
## 4 ABD19518 M. musculus Apxl mShroom2 Mus
## 5 AAF13269 M. musculus ShroomL mShroom3a Mus
## 6 AAF13270 M. musculus ShroomS mShroom3b Mus
## 7 NP_065910 H. sapiens Shroom hShroom3 Homo
## 8 ABD59319 X. laevis Shroom-like xShroom3 Xenopus
## 9 NP_065768 H. sapiens KIAA1202 hShroom4a Homo
## 10 AAK95579 H. sapiens SHAP-A hShroom4b Homo
## 11 ABA81834 D. melanogaster Shroom dmShroom Drosophila
## 12 EAA12598 A. gambiae Shroom agShroom mosquito
## 13 XP_392427 A. mellifera Shroom amShroom bee
## 14 XP_783573 S. purpuratus Shroom spShroom sea urchinAligning multiple sequences
TODO: in a brief sentence explain what the $ allows us to do
The $ allows us to access a specific column in the table; in this case, the column with the accession numbers can be accessed.
shroom_table$accession## [1] "CAA78718" "NP_597713" "CAA58534" "ABD19518" "AAF13269" "AAF13270"
## [7] "NP_065910" "ABD59319" "NP_065768" "AAK95579" "ABA81834" "EAA12598"
## [13] "XP_392427" "XP_783573"TODO: briefly explain what this chunk is doing and add the correct function
Takes all 14 sequences from the shroom_table vector with the imported package rentrez
# add necessary function
shrooms <-rentrez::entrez_fetch(db = "protein",
id = shroom_table$accession,
rettype = "fasta")TODO: in a very brief sentence explain what this is doing.
This concatenates the entire vector of characters (amino acids) inside shrooms and prints it out like it is in a text editor
cat(shrooms)TODO: in a brief sentence explain what this is doing and if/how its different from the previous code chunks
This chunk downloads the data of the 14 sequences in a different way. Instead of obtaining the sequence data from an imported package in rentrez, this uses a function written by the user called entrez_fetch_list, which is what makes it different from the previous chunk
shrooms_list <- entrez_fetch_list(db = "protein",
id = shroom_table$accession,
rettype = "fasta")TODO: briefly explain what I am doing this
This chunk is displaying the length of the object shrooms_list
length(shrooms_list)## [1] 14TODO: briefly explain what I am doing this. We will get into the details of for() loops in R later in the semester.
This section of code is cycling through every element (fasta sequence) in shrooms_list and cleaning it up.
for(i in 1:length(shrooms_list)){
shrooms_list[[i]] <- fasta_cleaner(shrooms_list[[i]], parse = F)
}TODO: summarize what is going on in this code chunk, then annotate each line of code with what its doing
In this code chunk, the variable shrooms_vector is being initialized to the length of shrooms_list and then is being filled with shrooms_list’s elements.
# Initializing shrooms_vector to the length of shrooms_list
shrooms_vector <- rep(NA, length(shrooms_list))
# Copying shrooms_list into shrooms_vector
for(i in 1:length(shrooms_vector)){
shrooms_vector[i] <- shrooms_list[[i]]
}
# Copying names of objects in shrooms_list into shrooms_vector
names(shrooms_vector) <- names(shrooms_list)TODO: explain what this is doing then add the necessary function.
This chunk converts the vector into a string set
# add necessary function
shrooms_vector_ss <-Biostrings::AAStringSet(shrooms_vector)MSA
TODO: briefly summarize what this section of the document will do.
Readings will be assigned to explain what MSAs are.
This section of the document will be doing a type of alignment called multiple sequence alignment, where three or more sequences of similar length are aligned; in this case, the shroom sequences are aligned together.
Building an XXXXXXXX (MSA)
TODO: briefly explain what this chunk does, then add the necessary function.
This chunk is building a multiple sequence alignment using ClustalW to align homologous nucleotide protein sequences.
# add necessary function
library(msa)
shrooms_align <- msa(shrooms_vector_ss,
method = "ClustalW")## use default substitution matrixViewing an MSA
TODO: briefly summarize what this section will do.
This section will allow the user to view the aligned sequences from msa. #### Viewing an MSA in R
TODO: Briefly summarize what output is shown below
The output that is shown here is the result of the multiple sequence alignment that happened on the shrooms_vector_ss object with all the fasta files in it.
shrooms_align## CLUSTAL 2.1
##
## Call:
## msa(shrooms_vector_ss, method = "ClustalW")
##
## MsaAAMultipleAlignment with 14 rows and 2252 columns
## aln names
## [1] -------------------------...------------------------- NP_065768
## [2] -------------------------...------------------------- AAK95579
## [3] -------------------------...SVFGGVFPTLTSPL----------- AAF13269
## [4] -------------------------...SVFGGVFPTLTSPL----------- AAF13270
## [5] -------------------------...CTFSGIFPTLTSPL----------- NP_065910
## [6] -------------------------...NKS--LPPPLTSSL----------- ABD59319
## [7] -------------------------...------------------------- CAA58534
## [8] -------------------------...------------------------- ABD19518
## [9] -------------------------...LT----------------------- NP_597713
## [10] -------------------------...------------------------- CAA78718
## [11] -------------------------...------------------------- EAA12598
## [12] -------------------------...------------------------- ABA81834
## [13] MTELQPSPPGYRVQDEAPGPPSCPP...------------------------- XP_392427
## [14] -------------------------...AATSSSSNGIGGPEQLNSNATSSYC XP_783573
## Con -------------------------...------------------------- ConsensusTODO: briefly explain what is being done in this chunk. This is tricky (and annoying) so do your best
The first chunk is taking the shrooms_align class and assigning it to be a multiple sequence alignment
The second chunk is converting the shrooms_align_seqinr from a sequence alignment file to different formats that can be used by other packages.
# WHAT IS THE LINE BELOW DOING? (its tricky - do your best)
class(shrooms_align) <- "AAMultipleAlignment"
# WHAT IS THE LINE BELOW DOING? This is simpler
shrooms_align_seqinr <- msaConvert(shrooms_align, type = "seqinr::alignment")TODO: what is the output this produces
It prints out the sequences and their corresponding amino acids; they’re getting closer and closer to multiple sequence alignment.
print_msa(alignment = shrooms_align_seqinr,
chunksize = 60)Displaying an MSA XXXXXXXX
TODO: explain this output and how its differnet from the prevoius
This output is showing a part of the alignment that starts 2000 amino acids in and goes for a length 100 amino acids. The alignment is completely done at this point and the colors show up where the amino acids align up the best.
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 ggplot2ggmsa::ggmsa(shrooms_align, # shrooms_align, NOT shrooms_align_seqinr
start = 2000,
end = 2100) 
Saving an MSA as PDF
TODO: explain what this command is doing. Add the package the function is coming from using :: notation This may not work for everyone. If its not working you can comment it out.
This command is taking the aligned shrooms sequences and making plots out of them and turning them into a pdf.
msaPrettyPrint(shrooms_align, # alignment
file = "shroom_msa.pdf", # file name
y=c(2000, 2100), # range
askForOverwrite=FALSE)TODO: explain what this command is doing This command is getting the filepath of the current working directory of the rmd file.
getwd()## [1] "C:/Users/thula/OneDrive/Documents/R"