readGENBANK Tutorial

Example Workflow for Gasterosteus (Stickleback) Sequences

Create character vector of accession numbers of Gasterosteus sequences

# Use paste() function to create a chr vector of accession numbers for Gasterosteus sequences
# These sequences all belong to one genus of sticklebacks
# Change in the tutorial to be sequences previous Endicott
# Bioinformatics students uploaded MT103163-MT103183
seq1 <- paste("JQ", seq(983161, 983255), sep = "") #paste is similar to c(), but output is a string instead of vector

Download all sequential sequences from Genbank

# This would be really hard to do my hand
# Note that the downloaded sequences are stored in a single variable called a list
sequences <- read.GenBank(seq1,
                          seq.names = seq1,
                          species.names = TRUE,
                          as.character = TRUE)

Write sequences to fasta file

# Write fasta file
write.dna(sequences, "fish.fasta", format = "fasta")

Workflow for Bicyclus anisops MT-CO1 Sequences

# Read in the downloaded NCBI accession list 
accList <- read.table("sequence.seq")

# See how data frame is structured to understand what to call in following code
## Leading letters are KM, but accession numbers are out of order
accList

##            V1
## 1  KM984276.1
## 2  KM984274.1
## 3  KM984271.1
## 4  KM984268.1
## 5  KM984267.1
## 6  KM984266.1
## 7  KM984265.1
## 8  KM984264.1
## 9  KM984256.1
## 10 KM984255.1
## 11 KM984254.1
## 12 KM984253.1
## 13 KM984252.1
## 14 KM984231.1
## 15 KM984230.1
## 16 KM984229.1
## 17 KM984228.1
## 18 KM984227.1
## 19 KM984226.1
## 20 KM984225.1
## 21 KM984223.1
## 22 KM984222.1
## 23 KM984221.1
## 24 KM984220.1
## 25 KM984219.1
## 26 KM984218.1
## 27 KM984217.1
## 28 KM984215.1
## 29 KM984213.1
## 30 KM984211.1
## 31 KM984207.1
## 32 KM984205.1
## 33 KM984202.1
## 34 KM984190.1
## 35 KM984188.1
## 36 KM984185.1
## 37 KM984182.1
## 38 KM984181.1
## 39 KM984180.1
## 40 KM984179.1
## 41 KM984178.1
## 42 KM984170.1
## 43 KM984169.1
## 44 KM984168.1
## 45 KM984167.1
## 46 KM984166.1
## 47 KM984145.1
## 48 KM984144.1
## 49 KM984143.1
## 50 KM984142.1
## 51 KM984141.1
## 52 KM984140.1
## 53 KM984139.1
## 54 KM984137.1
## 55 KM984136.1
## 56 KM984135.1
## 57 KM984134.1
## 58 KM984133.1
## 59 KM984132.1
## 60 KM984131.1
## 61 KM984129.1
## 62 KM984127.1
## 63 KM984125.1
## 64 KM984121.1
## 65 KM984119.1
## 66 KM984116.1

Create character vector of accession numbers of sequences

# Use paste() function to create a chr vector of accession numbers for Bicyclus anisops MT-CO1 sequences
BAseq <- paste("KM", accList, sep = "") 

# The string looks messy, so need to double check
print(head(BAseq))

## [1] "KMc(\"KM984276.1\", \"KM984274.1\", \"KM984271.1\", \"KM984268.1\", \"KM984267.1\", \"KM984266.1\", \"KM984265.1\", \"KM984264.1\", \"KM984256.1\", \"KM984255.1\", \"KM984254.1\", \"KM984253.1\", \"KM984252.1\", \"KM984231.1\", \"KM984230.1\", \"KM984229.1\", \"KM984228.1\", \"KM984227.1\", \"KM984226.1\", \"KM984225.1\", \"KM984223.1\", \"KM984222.1\", \"KM984221.1\", \"KM984220.1\", \"KM984219.1\", \"KM984218.1\", \"KM984217.1\", \"KM984215.1\", \"KM984213.1\", \"KM984211.1\", \"KM984207.1\", \"KM984205.1\", \"KM984202.1\", \"KM984190.1\", \"KM984188.1\", \"KM984185.1\", \n\"KM984182.1\", \"KM984181.1\", \"KM984180.1\", \"KM984179.1\", \"KM984178.1\", \"KM984170.1\", \"KM984169.1\", \"KM984168.1\", \"KM984167.1\", \"KM984166.1\", \"KM984145.1\", \"KM984144.1\", \"KM984143.1\", \"KM984142.1\", \"KM984141.1\", \"KM984140.1\", \"KM984139.1\", \"KM984137.1\", \"KM984136.1\", \"KM984135.1\", \"KM984134.1\", \"KM984133.1\", \"KM984132.1\", \"KM984131.1\", \"KM984129.1\", \"KM984127.1\", \"KM984125.1\", \"KM984121.1\", \"KM984119.1\", \"KM984116.1\")"

# Separate "KM" from digits and the numbers that follow after the dot in each accession number
clean_BAseq <- str_extract_all(BAseq, "KM\\d+\\.\\d+")[[1]]

# Check to see if the string is tidy
print(clean_BAseq)

##  [1] "KM984276.1" "KM984274.1" "KM984271.1" "KM984268.1" "KM984267.1"
##  [6] "KM984266.1" "KM984265.1" "KM984264.1" "KM984256.1" "KM984255.1"
## [11] "KM984254.1" "KM984253.1" "KM984252.1" "KM984231.1" "KM984230.1"
## [16] "KM984229.1" "KM984228.1" "KM984227.1" "KM984226.1" "KM984225.1"
## [21] "KM984223.1" "KM984222.1" "KM984221.1" "KM984220.1" "KM984219.1"
## [26] "KM984218.1" "KM984217.1" "KM984215.1" "KM984213.1" "KM984211.1"
## [31] "KM984207.1" "KM984205.1" "KM984202.1" "KM984190.1" "KM984188.1"
## [36] "KM984185.1" "KM984182.1" "KM984181.1" "KM984180.1" "KM984179.1"
## [41] "KM984178.1" "KM984170.1" "KM984169.1" "KM984168.1" "KM984167.1"
## [46] "KM984166.1" "KM984145.1" "KM984144.1" "KM984143.1" "KM984142.1"
## [51] "KM984141.1" "KM984140.1" "KM984139.1" "KM984137.1" "KM984136.1"
## [56] "KM984135.1" "KM984134.1" "KM984133.1" "KM984132.1" "KM984131.1"
## [61] "KM984129.1" "KM984127.1" "KM984125.1" "KM984121.1" "KM984119.1"
## [66] "KM984116.1"

Download sequences from GenBank using clean string

BAsequences <- read.GenBank(clean_BAseq,
                          seq.names = clean_BAseq,
                          species.names = TRUE,
                          as.character = TRUE)

Write sequences to fasta file

# Write fasta file
write.dna(BAsequences, "BA.fasta", format = "fasta")