Step 1: Applying Packages Needed ——————————————-

Loading in ape package to use when importing FASTA, FASTQ, and VCF files

Step 2: Reading in FASTA File ———————————————-

Using read.FASTA will read in the file, and specifying the type of data to DNA, although DNA is default for this command. The head command will print a report describing the data type and how it is formatted. It will also display the labels within the list. The structure (str) command will give an output of the values within the list with their associated labels.

fasta <- read.FASTA('MT103168.fasta', type = 'DNA')
# Reads in FASTA File specifying DNA sequences
head(fasta)
## 1 DNA sequence in binary format stored in a list.
## 
## Sequence length: 1560 
## 
## Label:
## MT103168.1 Bifidobacterium longum strain BB536 cell division...
## 
## Base composition:
##     a     c     g     t 
## 0.156 0.319 0.289 0.236 
## (Total: 1.56 kb)
# Output that describes the fastq file including data type and format as well as labels
str(fasta)
## List of 1
##  $ MT103168.1 Bifidobacterium longum strain BB536 cell division protein FtsW (rodA) gene, complete cds: raw [1:1560] 88 18 48 88 ...
##  - attr(*, "class")= chr "DNAbin"
# Output that shows the values of each string within the list with associated labels

Step 3: Reading in fastq file ———————————————-

Similar to using read.FASTA, the commands are the same, except there is no ‘type’ input for the read.fastq command

fastq <- read.fastq('ERR1072710.fastq')
# Reads in the fastq file
head(fastq)
## 3 DNA sequences in binary format stored in a list.
## 
## Mean sequence length: 183.667 
##    Shortest sequence: 146 
##     Longest sequence: 259 
## 
## Labels:
## ERR1072710.1 10317.000001315_0 length=151
## ERR1072710.2 10317.000001315_1 length=116
## ERR1072710.4 10317.000001315_3 length=151
## 
## Base composition:
##     a     c     g     t 
## 0.318 0.208 0.254 0.219 
## (Total: 551 bases)
# Output that describes the fastq file including data type and format as well as labels
str(fastq)
## List of 3
##  $ ERR1072710.1 10317.000001315_0 length=151: raw [1:146] 18 18 88 88 ...
##  $ ERR1072710.2 10317.000001315_1 length=116: raw [1:259] 18 28 18 28 ...
##  $ ERR1072710.4 10317.000001315_3 length=151: raw [1:146] 28 28 88 28 ...
##  - attr(*, "class")= chr "DNAbin"
##  - attr(*, "QUAL")=List of 7
##   ..$ ERR1072710.1 10317.000001315_0 length=151: num [1:11] 32 38 51 34 32 34 32 34 32 38 ...
##   ..$ ERR1072710.2 10317.000001315_1 length=116: num [1:11] 30 30 30 30 30 30 30 30 30 30 ...
##   ..$ ERR1072710.4 10317.000001315_3 length=151: num [1:42] 10 36 49 49 16 15 22 17 22 16 ...
##   ..$ NA                                       : num [1:70] 51 32 34 38 38 32 38 38 38 51 ...
##   ..$ NA                                       : num [1:67] 30 30 30 30 30 30 30 30 30 30 ...
##   ..$ NA                                       : num [1:11] 32 51 51 32 38 32 38 34 34 51 ...
##   ..$ NA                                       : num [1:11] 30 30 30 30 30 30 30 30 30 30 ...
# Output that shows the values of each string within the list with associated labels

Step 4: Reading in a vcf file as a dataframe ——————————

Reading in vcf file can be done in different ways and this method uses read.csv. Using sep = ‘ separates each line of the VCF into its own row, and setting comment.char =’#’ will remove the metadata comments in the beginning of the file. Setting header = FALSE prevents the first row of variant data from being used as the column headers, as it is the first row that is not preceded by a #. colnames(df) was used to manually enter in the column headers that could be found in the VCF. The ‘Results’ column header was created to connect the ‘Format’ column to the results it displays in the final column of the dataframe. df[df] == ‘.’ <- changes missing data in the VCF dataframe to NAs to clean up the data.

variants <- read.csv('TwoVariants.vcf', sep = '\t', header = FALSE, comment.char = '#')
# Read vcf in as a csv, separating each line into a row with sep = \t. By setting comment.char = #, will remove metadeta comments in the beginning of the vcf file, and header = FALSE prevents the first variant currently in row 1 as being the header line.
colnames(variants) <- c('Chromosome', 'Position', 'ID', 'Ref', 'Alt', 'Qual', 'Filter', 'Info', 'Format', 'Results' )
# Manually changing column headers per the vcf file, with the 'Results' column reflecting the order in which the 'Format' column is displayed
variants[variants == '.'] <- NA
# Data clean up step to replace '.' where there is missing data to 'NA' to represent data better