# BI278 Lab #2 - NCBI data, kmers, and prokka

Excerise 1 - Download public genomic data from NCBI
1.1 - How to navigate NCBI for genomic data

- You can search thru projects to find data sets based on various attributes (thru "filters" button):
- types of organisms
- types of sequences (imited to various types of next gen sequencing reads)

- usually "genome sequencing" or "raw sequence reads" in Data type filter
- can then sort by Data Type (click on arrow in column headings)

- once you find a dataset of interest (data + publication is good), you click on link
- clicking "Genome" will take you to the link --> it's the fastest way to get to a genome that's already put together
- second genome link ("genome") will download the chromosome level genome sequence for you onto desktop
- if you want individual protein sequences, you can click on the number (usually provided in table)

1.1b - Finding a bacteria that has both a genome and raw sequence reads available thru NCBI

Bacteria I picked: Burkholderia pseudomallei
Registration date: 2015-10-16
Link: https://www.ncbi.nlm.nih.gov/genome/476

GC % - 4836611
Genome Size - 7085397

Remember to unzip files first! And check for the correct path!

Path: /home2/sdivit25/lab_02
Unzip file: gunzip "filename"

1.2 - Download Raw Read via SRA Toolkit

NCBI SARS CoV2 genome sequencing example
BioProject: https://www.ncbi.nlm.nih.gov/bioproject/?term=NYU+Langone+SARS-CoV-2+Sequencing+team

Then click # of links for BioSample.
Link: https://www.ncbi.nlm.nih.gov/biosample?Db=biosample&DbFrom=bioproject&Cmd=Link&LinkName=bioproject_biosample&LinkReadableName=BioSample&ordinalpos=1&IdsFromResult=650245

Then search for ID # & click on it.
Link: https://www.ncbi.nlm.nih.gov/sra/SRX8878313/

From there, click on "Run".
Link: https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR12379546&display=metadata

Make sure to checksize of file & length of sequencing read. The ID you need for next part is the SRR run # (usually starts w/ "SSR").
- SSR # = under "Runs" column
- for this COVID file: SRR12379546

- for the genomic (Burkholderia pseudomallei) file: SRR2134242

LINK:
raw seq - https://www.ncbi.nlm.nih.gov/bioproject/PRJNA290528

SRA exp tells you that raq seq data is available:
https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=290528
click the first one:
https://www.ncbi.nlm.nih.gov/sra/SRX1124256[accn]
GIVES YOU RUN NUMBER:
https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR2134242&display=metadata

SRAtoolkit = commmand line program we're using to download some genomic data from NCBI
- sysadmin installed this on bi278

**Command to download NCBI SRA (Sequence read archive) database:

fastq-dump

In order to run this^:

vdb-config --interactive

BUT it's imp to know exactly what that's doing so everytime you run a new software, find out more info:

fastq-dump --help

To check what the downloaded data may lok like before downloading the whole file:

fastq-dump -X 3 -Z SRR12379546
#sequence data of the first 3 spots will
#be written to our screen instead of a
#file b/c of the other optional flags -X
#and -Z used here

^Doing this step, you will notice:
- these seqs are supposed to be 2 101bp length seqs per "spot"
- BUT output --> these 2 101bp length seqs are being combined
- so we want to download this data w/ the pairs separated:

fastq-dump -X 3 --split-3 --skip-technical --readids --read-filter pass --dumpbase --clip -Z SRR12379546

You know it worked if the length is roughly half. So before for the first fastq-dump, length was 302 and after the second fastq-dump, length was 151!

Briefly:
Format of FASTQ files = units of 4 lines at a time indicate for a single read the DNA sequence (marked by @) and individual base quality scores (marked by +)

**Important to do these types of checks b/c there is lot of data in SRA that are not deposited/described properly.
^Easy way of avoiding these bad data is by checking whether the traces and the metadata of a sample of the files match up.

Now, you have a genome for in the prev steps, so download raw sequencing data reads but w/out quality scores:

fastq-dump --split-3 --skip-technical --readids --read-filter pass --dumpbase --clip -v --fasta default --outdir ~/lab_02 SRR12379546

To look at first 20 lines of FASTA files:

head -20 PATH/*.fasta

#but no path needed if you're already in
#that folder/directory

so then:
head - 20 *.fasta

**CHECK ENDING OF ".fasta" to make sure it
matches with what you're looking for

**Practice** --> download one of other organisms other than P. fungorum (overlap w/ last week)

|Organism|SRA Instrument [type] record|SRA Run Number|Genome Size (bp)|Estimated size|
| -------- | -------- | --- | --- | -------- |
|P. Fungorum|Illumina NovaSeq 6000|SRR11022347|9058983|9185185|
|P. sprentiae|Illumina GA IIx|SRR3927471|7829542| |
|**P. terrae**|PacBio RS II|DRR322713|10062489| **COULDN'T CALCULATE BC OF MEMORY ISSUE|
|P. xenovorans|Illumina HiSeq 2000|SRR2889773|9702951| |
|Your organism: Burkholderia pseudomallei| Illumina|SRR2134242|7085397|7800419|

Link for P. terrae: https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=776439

Exercise 2 - Count k-mers & estimate genome size based on k-mer frequencies
2.1 Count Kmers w/ Jellyfish

- program used to count k-mers in our genome is called jellyfish

- commands you want to use:

jellyfish count
jellyfish histo

typing " --help" for additional usage information

Steps:
1) Choose a k-mer size between 19-31
My chosen size: 20
2) run "jellyfish count" for genome you will compare

jellyfish count -t 2 -C -s 1G -m kmerSize -o name.m29.count input_file

jellyfish count -t 2 -C -s 1G -m 20 -o B.PseudomalleiIllumina.count B.PseudomalleiIllumina.fasta

3) run "jellyfish histo" on the count file to get the frequency distribution across all kmers
a) get rid of count file after checking that it worked b/c it takes up a lot of memory and you don't need it anymore

jellyfish histo -o name.m29.histo name.m29.count
jellyfish histo -o B.PseudomalleiIllumina.histo B.PseudomalleiIllumina.count

4) Look at histo file. 2 choices:

cat ==> will print whole file to screen
less ==> will show you one screen & a line at a time, type 'q' to exit

**YOUR FILE WILL SHOW TWO COLUMNS OCCURENCE & DENSITY.
- occurence = frequency, how often a k-mer appears in that file
- density = how many kmers show at that level of occurence

**UNDERSTANDING IT:
There are 4^20 kmers available of size 20.
So 23 254926 --> there are 254926 of the 4^20 kmers that appear 23 times
- NOTICE: lot of kmers that occur only once
- these are most likely to be result of errors ==> thus, used to remove reads from datasets that are likely errors
- Quake is a software that does that

5) Now that your code is working: write a shell script to do the same for other SRA downloaded file.

jellyfish count -t 2 -C -s 1G -m 20 -o $1.count $1
jellyfish histo -o $1.histo $1.count
rm $1.count

named: ColumnOcurrence.sh

2.2 - Visualize your Kmer counts & estimate genome size in R
1) open colby server: bi278.colby.edu
2) set working directory your histo files are in

getwd()
setwd("/home2/sdivit25/lab_02")
3) read in one by one the result files you created in the prev step

name <- read.table("name.m29.histo", h=F, sep=" ")
B.PseudomalleiIllumina <- read.table ("B.PseudomalleiIllumina.histo", h=F, sep=" ")

4) check out shape of histogram

plot(name, type ="l")
plot(B.PseudomalleiIllumina, type="l")

Since, we want to find the midpoint of
the main peak that represents the single
copy regions of he genome, we can narrow
down the range of the histogram to try
to figure out what the x-axis value is
at the main peak:

plot(name[5:250,], type="l")
plot(B.PseudomalleiIllumina[5:250,], type="l")

5) we can check out the data directly to determine the midpoint of the main peak based on where it was located in your plot. ***THIS MIDPOINT NUMBER REPRESENTS A ROUGH ESTIMATE OF GENOME COVERAGE.

name[150:180,]
B.PseudomalleiIllumina[0:50,]

My midpoint with the highest density was about 23 (254926)

6) Now, we want to area under the curve & divide it by this coverage estimate to get an estimate of genome size. This estimation is made possible because the distribution of k-mers in a genome follows a Poisson distribution. We're ignoring the part of the curve (inital peak) that contains the results of sequencing error:

sum(name[5:nrow(name), 1]*name[5:nrow(name),2])/"midpoint with highest density"

sum(B.PseudomalleiIllumina[5:nrow(B.PseudomalleiIllumina), 1]*B.PseudomalleiIllumina[5:nrow(B.PseudomalleiIllumina),2])/23

Explanation:
nrow(name)
head(name)
head(name[,1]) #square brackets are used to limit [row, column]
head(name[,2])
head(name[5:nrow(name),])
head(name[5:nrow(name),1]*name[5:nrow(name),2])
sum(name[5:nrow(name),1]*name[5:nrow(name),2])

***NOTE:

Professor Noh told us that we didn't have to do Exercise 3.