Welcome to the third day of the [FAS Informatics](https://informatics.fas.harvard.edu/) [Bioinformatics Tips and Tricks Workshop](https://harvardinformatics.github.io/workshops/2023-spring/biotips/)!

If you're viewing this file on the website, you are viewing the final, formatted version of the workshop. The workshop itself will take place in the RStudio program and you will *edit and execute the code in this file*. Please download the raw file [here](https://harvardinformatics.github.io/workshops/2023-spring/biotips/Biotips-workshop-2023_Day3-student.Rmd)

Today we're going to continue our tour and explanation of common genomics file formats and their associated tools by talking about variant files, that is files which indicate where variants in some sample relative to a reference are.

We'll also introduce some basic concepts about bash scripting today, in preparation for tomorrow when we'll talk about using Slurm and submitting jobs to the cluster.

# Setting up our link to the data files

```{bash}

mkdir -p data3
ln -s -f /n/holylfs05/LABS/informatics/Everyone/workshop-data/biotips-2023/day3/* data3
# ln: The Unix link command, which can create shortcuts to folders and files at the provided path to the second provided path
# -s: This option tells ln to create a symbolic link rather than a hard link (original files are not changed)
# -f: This option forces ln to create the link

ls -l data3
# Show the details of the files in the new linked directory

```

# VCF files

Another important type of genomic coordinate file is the **VCF** file, or *variant call format* file. This is another tab-delimited file that contains **variants** between a single genome or a sample of genomes and a reference genome. These variants can include single nucleotide changes, small insertions and deletions, or even larger structural variants. The most common type of variation you'll see in a **VCF** file are the single nucleotide ones, so the content of a **VCF** file is commonly referred to as *SNP calls*.

We'll now be working with a **VCF** file that contains SNP data from some Amazon molly (*Poecilia formosa*), a freshwater fish. The original data consists of hundreds of thousands of SNPs across the molly genome among 16 individuals, but we'll be looking only at SNPs on a single scaffold today (about 100,000).

Let's take a look at our file and explain a bit about the format.

> Run the code block below to view the beginning of our **VCF** file:

```{bash}

head -n 43 data3/poeFor_NW_006799939.vcf 
# head: A Unix command to display the first few lines of a file
# -n 43: This tells head to display exactly the first 43 lines of the provided file

```

So what we see here at the top of the file is the **header** of the VCF. Much like headers in **BAM** and **SAM** files, headers in **VCF** files contain information about how this file was created and the contents of the file. For example, the lines that begin with `##FORMAT` describe the contents of a particular column of the file, while the lines begin with `##FILTER` describe what the filter labels in the *Filter* column mean. Unlike the headers in **BAM** and **SAM** files, which begin with the `@` character, the headers in a **VCF** file begin with the `##` string (like GFF files).

Now that we know a bit about the header, let's look at the actual content of the file, the tab-delimited rows containing SNP calls.

> Run the code block below to view the beginning of our **VCF** file while skipping header lines:

```{bash}

grep -v "##" data3/poeFor_NW_006799939.vcf | head -n 2
# grep: The Unix string search command
# "##": The string to search for in the provided file - we just don't want to display these for this demonstration
# head: A Unix command to display the first few lines of a file
# | : The Unix pipe operator to pass output from one command as input to another command
# -n 2: This tells head to display exactly the first 2 lines of the provided file
## NOTE: ignore the "Broken pipe" and "error closing standard output" errors! This is just an artifact of our setup for the workshop today.

```

The first line here, beginning with a single `#` shows us the names of the columns, which may be referred to in the file header like we saw above. Each subsequent line represents one **variant** relative to the reference genome, with the following column definitions:

1.  CHROM: The chromosome or assembly scaffold in the reference genome.
2.  POS: The position of the variant in the reference genome.
3.  ID: A name for the variant. If the variant has no name, just the `.` character.
4.  REF: The allele in the reference genome.
5.  ALT: The alternate allele in the the sample.
6.  QUAL: The quality score of the variant call. How this is calculated depends on the variant calling program.
7.  FILTER: The filter label for the current variant. If the variant isn't filter, this will be `PASS` (or sometimes `.`).
8.  INFO: Additional info for the current variant. Usually a semi-colon (`;`) delimited list of information.
9.  FORMAT: A colon (`:`) delimited list of the keys for the information in the subsequent columns. These keys are defined in the file header. 10+. SAMPLE GENOTYPES: Each subsequent column contains a colon (`:`) delimited list that follows the keys in the FORMAT column for the genotypes of the samples. Each column corresponds to one sample with the name being the column header.

The genotypes are encoded in the sample columns with the `GT` key. In this file, it is the first entry. For example, with the keys defined in column 9, the genotype of the first sample:

    ./.:0,0:0:.:.:.:0,0,0:.

means that the GT (GENOTYPE from the header) is `./.`, the AD (ALLELIC DEPTH) is 0,0, the DP (DEPTH) is 0, and so on.

Genotypes use the following syntax, with `0` being the reference allele and `1` being the alternate allele:

-   `0/0` means this sample is homozygous for the reference allele.
-   `0/1` or `1/0` means this sample is heterozygous for the reference allele (0) and the alternate allele (1).
-   `1/1` means this sample is homozygous for the alternate allele.
-   `.` is used for missing; this site was not genotyped for this individual, which happens to be the case for the sample outlined above.

For **VCF** files with multiple samples (like this one), sites can have more than two alleles, in which case the genotypes `0/2` or `0/3` or `2/2` may be seen.

There are many quirks to **VCF** files, but they are also highly standardized, so if you find yourself working with them be sure to read the very detailed specifications here:

[VCF file specs](https://samtools.github.io/hts-specs/VCFv4.4.pdf)

### Note that the binary version of **VCF** files are **BCF** files, usually with the `.bcf` file extension. BCF files are more efficient to read than VCF files, and are typically compressed.

There are a myriad of ways to use and analyze a **VCF**, but for general processing of them, we will use `bcftools`. [bcftools](https://samtools.github.io/bcftools/bcftools.html) falls under the same development umbrella as `samtools` and thus shares much of the same usage and underlying philosophy to act like a native Unix command. It was preceded by [vcftools](https://vcftools.github.io/), which still retains some useful functionality, but is no longer developed.

# bcftools

Much like `bedtools`, `bcftools` has a wide range of functions so we won't be able to get to everything today. Be sure to read [the documentation](https://samtools.github.io/bcftools/bcftools.html) to find out if it can accomplish the task you need with your **VCF** or **BCF** file.

One of the most basic things you can do to a file is to **view** it. Indeed, *looking at your data* is an important part of each step of a bioinformatics workflow.

Like `samtools`, `bcftools` has a dedicated `view` command.

> **Exercise**: 
> 1. Run the code block below to `view` the top of the Amazon molly **VCF** file. 
> 2. Edit the `bcftools view` command to exclude the header from the output.

```{bash}

## Use the view command to display the VCF file, then edit it to exclude the header
## data3/poeFor_NW_006799939.vcf
bcftools view data3/poeFor_NW_006799939.vcf | head -n 10
## Use the view command to display the VCF file, then edit it to exclude the header
## NOTE: ignore the "cannot write" error! This is just an artifact of our setup for the workshop today.

```

So `view` acts much the same way as the Unix command `cat`, but has specific functionality for the **VCF** format, and also handles compressed **BCF** files natively.

One useful function when `view`ing files is to be able to subset or filter them. We've provided a file that lists the sample IDs for a made-up population and can use `bcftools view` to extract genotypes from specific samples.

> Run the code block below to extract SNPs only from a subset of samples:

```{bash}

cat data3/pop1.txt
echo "-------"
# cat: A Unix command that displays the contents of a file in the Terminal
# echo: A Unix command that simply prints the provided input to the screen
# These two commands are just to see the pop1.txt file

bcftools view -S data3/pop1.txt -H data3/poeFor_NW_006799939.vcf | head
# bcftools: A suite of programs to process VCF files
# view: The sub-program of bcftools to execute
# -S : A bcftools view option that specifies to only extract genotype columns from the samples listed in the provided file
# -H : A bcftools view option that specifies to exclude the VCF header from the output
## NOTE: ignore the "cannot write" error! This is just an artifact of our setup for the workshop today.

```

### A note about subsetting

This can be very handy, but be warned - the INFO fields are not automatically recalculated, and will no longer be correct once you remove some individuals. In general, the INFO column and its information displays information about the *site* while the FORMAT fields display information about each *sample*.

# SNP density with `bcftools` and `bedtools`

Many of the tools we've demonstrated, `samtools`, `bedtools`, `bcftools`, are general purpose tools for basic processing and manipulation of specific file formats used in genomics. Many of the tasks they do could be done with native `bash` commands, but that would require a lot of piping and probably some custom scripting on uncompressed files. However, if you really start to get to know these tools and their versatility, you can do some more advanced analyses purely with them. Let's demonstrate by calculating SNP density across 100k windows the specific 7 Mb scaffold we are working with from the Amazon molly genome.

To do this, we'll need:

1.  The **.fai index** file of the Amazon molly genome, which we have already generated.
2.  A **bed** file with the coordinates of the 100 kilobase windows in the scaffold
3.  A **VCF** file with Amazon molly SNPs, which we have provided.

Since we've already downloaded the genome and generated the index file, let's start with the second one, getting the window coordinates. We can do this easily with the `bedtools getwindows` command:

> Run the code block below to generate a bed file with 100kb windows in the Amazon molly genome:

```{bash}

bedtools makewindows -g data3/poeFor_NW_006799939.fasta.fai -w 100000 > poeFor-windows-100k.bed
# bedtools: A suite of programs to process bed files
# makewindows: The sub-program of bedtools to execute
# -g : A genome size file, which needs two columns, the scaffold/chromosome name and the size of that feature
# -w : The size of the windows to partition each scaffold in -g into
# -b : The second interval file to check overlaps
# > : The Unix redirect operator to write the output of the command to the following file

head poeFor-windows-100k.bed
# Display the first few lines of the bed file containing windows

```

Now, we can use `bedtools intersect`, which also works with **VCF** files, to count how many SNPs appear in each window.

> **Exercise**: Use `bedtools intersect` to count the number of SNPs in each 100kb window in the Amazon molly genome. You will need to look in the help menu of `bedtools intersect` to find the right option to output a column with the number of overlaps per region! Save this file as `poeFor-windows-100k-snps.bed`.

```{bash}

## Write a bedtools intersect command to count how many SNPs overlap each 100kb window in the molly genome
## data3/poeFor_NW_006799939.vcf
## poeFor-windows-100k.bed

## Write a bedtools intersect command to count how many SNPs overlap each 100kb window in the molly genome

head -n 10 poeFor-windows-100k-snps.bed
# Display the first few lines of the bed file containing windows and overlapping SNPs

```

> **BONUS Exercise**: If you attended the R workshop last week or know R already, use the code block below to read in the file we just generated with SNPs per 100k window and plot the density along this scaffold.

```{r}

library(tidyverse)
# Load the tidyverse

## Write code in R to plot SNP density across the scaffold in our bed file


## Write code in R to plot SNP density across the scaffold in our bed file


# Display the plot

```

# Scripts

So far, we've been focused on interactive analysis of small files. However, in real bioinformatics work, you will likely find yourself with many samples that you want to run the same commands on. And often these commands can require significant compute time, such that you want to be able to run these commands non-interactively. As nice as these markdown notebooks are for teaching and data exploration, they are not the ideal format for this task. Instead, what we want are scripts.

A bash script (like an R script, if you attended the R workshop earlier), is just a text file that executes a series of commands in order. Bash scripts can get exceedingly complicated, but today we are just going to introduce a few basic concepts.

We are ultimately going to make a script that takes as input a **VCF** file and one or more **bed** files, and reports as output the average SNP density in regions of the genome inside bed intervals. Let's start simple.

> Using the file menu, go to New file --\> Shell script and create a new file. Call it `snp-density.sh`.This is going to be our script in a new tab of our **text editor**, that we will modify as we go.

The first thing we need to do is a little bit of shell magic to get the shell to treat the script we are writing as a command. This is called a shebang (`#!`) and basically says "run the script with the following **interpreter**". Type the following line at the top of your script:

```
#!/bin/bash
```

This means that when we run `snp-density.sh` on the command line, the commands in the script will be interpreted as shell commands.

We need to do one more thing, which is make our script **executable**. We'll use the `chmod` command for this.

> **In the Terminal below**, make sure you are in the same directory as your script file, and then type the following command.

```
chmod +x snp-density.sh
```

We now have a script we can run.

> In the code block below, try running your (empty) script. Note that the `./` specifies the *relative path* to the script, and means "execute the script snp-density.sh in the current working directory". If omitted, you'll get a "command not found" error.

```{bash}

./snp-density.sh
# ./ : Execute the following file as a shell script

echo "done"
# echo: A Unix command that simply prints the provided input to the screen

```

However, nothing will happen if we do, since there are currently no commands in the file to run.

> **Exercise**:
> In order to have our shell script do anything, we need to put some commands in it.
> 1. Copy the `bedtools intersect` command from above into the script.
> 2. Instead of writing the output of the command to the file, **pipe** the output of the command to `awk`.
> 3. Use `awk`, to calculate the *average number of SNPs per base across intervals*. To do this, you'll need to increment two counters for each row of the output: one with the number of snps in that row, the other with the length of the element in that row, and then divide the number of snps by the total length.
> *Warning: `awk` has a built-in function called `length`, so you'll need to call your length variable something else*.

Now, try running your script again.

> In the code block below, try running your script.

```{bash}

./snp-density.sh
# ./ : Execute the following file as a shell script

```

We should get a number now.

However, we have a problem. Ideally, we'd like to be able to do this calculation on any **VCF** file and any **bed** file. Right now those files are *hard-coded*, which means to change what file the script operates on requires changing the script.

To make our script flexible, we need to introduce the concept of variables. We've seen this in `awk` commands; we can also create variables in the shell. We do this just by assigning something (text or the output of a command, usually) to a string.

> Run the code block below to create a variable called TEXT, and print it out with the `echo` command.

```{bash}

TEXT="Hello World!"
# = : Assign the variable name on the left to the value on the right

echo $TEXT
# echo: A Unix command that simply prints the provided input to the screen
# $ : A symbol that the shell understands to mean that the following is the name of a variable

```

Note that when we access what is stored in the variable we use the `$` notation, but we don't use this when we assign to a variable. Note also we cannot have spaces before or after the equal sign or the shell will thing we are trying to run the `TEXT` command, not create a variable called `TEXT`.

> **Exercise**:
> 1. Modify our shell script, `snp-density.sh`, so that it includes two variables, one called `VCF` and one called `BED`. For now, define these in the second and third lines of the script with the paths to the **BED** file and the **VCF** file we currently have.
> 2. Replace the hard-coded file names in the `bedtools intersect` command options `-a` and `-b` with the `BED` and `VCF` variables.
> 3. In the code block below run your script again and make sure it still works.

```{bash}

./snp-density.sh
# ./ : Execute the following file as a shell script

```

You should get the same answer as before. 

Note that right now, we aren't actually changing anything, we've just moved the files we are going to operate on to the top of the script, but we still would need to change the script each time we run it on different files. To make our script more *general* (able to take various inputs more easily), we should have a way to change the input file without having to edit the script. The way to do this is with **command line options**. You've seen this a lot in the tools we've run up to this point. For instance, when we run almost any shell command, it requires a filename, e.g.:

```

head poeFor-windows-100k.bed

```

In this case, the file name `poeFor-windows-100k.bed`, is a **command line option**, also known as a **command line argument** or just **option**, **argument**, or **parameter**. Commands and scripts can take multiple arguments. Some are **positional**, meaning that the way they are read into the program depends on their placement in the command. Others have specific **flags**, meaning the way they are read into the program depends on the flag they come after in the command. For instance:

```

head -n 20 poeFor-windows-100k.bed

```

This is a command that hopefully looks pretty understandable to you by now. This command, `head` has two **command line options** as we've written it:

1. A **flag**, `-n 20`, which changes the number of lines the command displays. Internally, this sets the value of some variable to a different number.
2. A **positional argument**, the file name `poeFor-windows-100k.bed`, which is the last listed argument. Again this internally sets the value of some variable to the path to that file (actually, for most bash commands, the first un-flagged option is usually reserved for the file name).

Well, we can modify our script to take input from the command line as well. Shell scripts have access to special variables that can be used to create simple, **positional**, command line options. Specifically, the variables `$1`, `$2`, `$3`, etc will contain the first, second, third, etc. **arguments** specified after the script is called.

> **Exercise**:
> 1. Modify your script to replace the path to the **bed** file with `$1`, and the path to the **VCF** file with `$2` using **positional command line arguments**.
> 2. Run the following code block to call your script again, this time taking arguments, with the bed file first and the vcf second.

```{bash}

./snp-density.sh poeFor-windows-100k.bed data3/poeFor_NW_006799939.vcf 
# ./ : Execute the following file as a shell script

```

If everything is working, you should get the same answer as before. But now we can give our script any **bed** file and any **VCF** file, right from the command line. Note that we don't do any error checking - if you accidentally specify a **VCF** file that uses a different reference genome than your **bed** file, nothing in the script will prevent this.

This is better, but it still requires a lot of typing if we want to run this on 10 files. If, for example, we wanted to compute SNP density separately for each chromosome, or for different interval types (e.g., genes, introns, exons), we'd have to type out each bed file separately.

We can get around this by using loops.

The most common kind of loop we'll use is a `for` loop, which iterates over elements.

> Run the following code block to see a for loop in action

```{bash}

# a simple loop
for i in 1 2 3;
do
  echo $i
done

```

Importantly, a loop can iterate over the output of another command.

> Run the following code block to iterate over the results of the `ls` command in a for loop

```{bash}

# use the $() syntax to capture the output of a command into a loop; by default each newline will be a separate iteration

for FILE in $(ls)
do
  echo $FILE
done
```

All this does right now is print out the files in our current working directory, just like running the `ls` command would. 

Loops are very powerful. For example, let's say we want to make a few different subsets of the annotation we have for the Amazon molly: coding exons, introns, intergenic regions, UTRs. 

Let's start by writing an `awk` command that will read through the **GFF** file associated with the Amazon molly reference genome, find all the exons in it, and make a **bed** file of just exons intervals. Since all the information for a **bed** file is also contained in a **GFF** file (scaffold, start, and end positions) we can convert a **GFF** to a **bed** with `awk`: 

```
BEGIN{OFS="\t"}; {print $1, $4-1, $5}
``` 

The `BEGIN{OFS="\t"}` tells `awk` to use a tab separator in the output (required for `bedtools`), and then we print the first, fourth, and fifth columns, remembering that **bed** and **GFF** files have different coordinate systems. 

> Run the code block below to convert the exon coordinates in the Amazon molly **GFF** file to **bed** format:

```{bash}

awk 'BEGIN{OFS="\t"}; $3=="exon"{print $1, $4-1, $5}' data3/poeFor_NW_006799939.gff > exon.bed
# awk: A command line scripting language command
# '' : Within the single quotes is the user defined script for awk to run on the provided file
# > : The Unix redirect operator to write the output of the command to the following file

head exon.bed

```

We could just copy and paste this command line to run this again to get **bed** files for other features like genes, mRNAs, CDS, etc. But we can also use a loop! Note that `awk`, by default doesn't have access to shell variables, but we can pass them to an `awk` script with the `-v` option.

```{bash}

mkdir -p poeFor_beds
# mkdir: The Unix command to make a new directory
# -p: This option tells mkdir to make any parent directories that do not exist

GFF="data3/poeFor_NW_006799939.gff"
# = : Assign the variable name on the left to the value on the right

for TYPE in exon CDS mRNA gene;
do
  awk -v type=$TYPE 'BEGIN{OFS="\t"}; $3==type {print $1, $4-1, $5}' $GFF | bedtools sort -i - | bedtools merge > poeFor_beds/$TYPE.bed
done
# This loop runs an awk script on several different inputs and pipes it to bedtools

ls poeFor_beds
# ls: The Unix list directory contents command

```

We are almost done; let's make a few more using `bedtools complement`, which reports all regions of the genome NOT in an interval, and `bedtools subtract`, which removes overlaps. For example, we don't have an `intergenic` feature in our GFF, but we can create one by getting everything that is not a gene. And we don't have a `UTR` feature or an `intron` feature, but we can get them by taking our exons and removing coding sequences (CDS) or by taking genes and removing exons.

```{bash}

bedtools complement -i poeFor_beds/gene.bed -g data3/poeFor_NW_006799939.fasta.fai > poeFor_beds/intergenic.bed
# bedtools: A suite of programs to process bed files
# complement: The sub-program of bedtools to execute
# -i: The input bed file
# -g:  A genome size file, which needs two columns, the scaffold/chromosome name and the size of that feature
# > : The Unix redirect operator to write the output of the command to the following file

bedtools subtract -a poeFor_beds/gene.bed -b poeFor_beds/exon.bed | bedtools sort -i - | bedtools merge > poeFor_beds/intron.bed
bedtools subtract -a poeFor_beds/exon.bed -b poeFor_beds/CDS.bed | bedtools sort -i - | bedtools merge > poeFor_beds/UTR.bed
# Two commands that subtract bed intervals in one file (-b) from another (-a) and pipe both with | and - into other bedtools commands and
# finally redirect to an output bed file.

wc -l poeFor_beds/*.bed
# Display the number of lines in each bed file in the poeFor_beds directory using wc -l and a wildcard (*)
# Note that here we have run a command on multiple files without an explicit loop, thanks to regex!

```

Now we have a directory of **bed** files, and we might want to run our script over each one, getting the average SNP density by feature type. We can do this with a loop.

> **Exercise**:
> 1. Modify the code block below to add a for loop that iterates over the **bed** files in the poeFor_beds directory, and for each bed file calls the snp density script.

```{bash}

VCF="data3/poeFor_NW_006799939.vcf"
# = : Assign the variable name on the left to the value on the right

# add a for loop that iterates over each bed file in the poeFor_beds directory

  echo "SNP density for $BED:"
  ./snp-density.sh $BED $VCF
  echo "---"

# add a for loop that iterates over each bed file in the poeFor_beds directory

# Loop over each bed file in our directory and run our snp-density script on it and the VCF file

```

We can actually make our script a bit better by putting the loop in the script, and passing a **file of files**, that is a single file that has all the bed files we want to operate on, one per line, as an option. 

> **Exercise**:
> 1. Re-write the `snp-density.sh` script to take as input a single file that contains a list of **bed** files (one per line) in addition to the **VCF** file and write a *loop* to run the SNP density `awk` command on each file in that file. *Hint: the `cat` command will return the lines in a file, so $(cat <FILENAME>) will let you iterate over the lines of a file in a for loop*
> 2. Run the code block below to generate the file containing all of our **bed** files and run the `snp-density.sh` script

```{bash}

VCF="data3/poeFor_NW_006799939.vcf"
# = : Assign the variable name on the left to the value on the right

ls poeFor_beds/*.bed > bed_files
# Redirect the output of ls to a file called bed_files

./snp-density.sh bed_files $VCF
# ./ : Execute the following file as a shell script

```

Progress! We can now calculate SNP density for any arbitrary set of **bed** elements, given a **VCF** file (or, we could even do this over multiple **VCF**s by using another loop).

However, this script is kind of brittle. What if we included a non-bed file in our list of files? Our script will fail cryptically. We might want to add some basic **error checking**. In particular, we are going to check that the files in the **file of bed files** actually exist. Note that it is possible, but complicated, to do more sophisticated error handling. In this case, we just want to make sure we don't have typos or path errors in our file list. 

To do this, we need to learn about one more bit of *control flow* code: **conditional statements**, or **if statements**. **If statements** are pretty intuitive - they just test `if` something is `True` (in the logical sense), and only run the following code if it is. We've seen a lot of things like this in `awk`, when only doing actions on line that match a certain pattern. Those are *implicit* `if` statements. Here is an example of an *explicit* if statement in **bash** syntax.

> Run the code block below to see how an `if` statement works in **bash**. Change the value of `SPEAK` to see how it affects the behavior o the statement.

```{bash}

SPEAK="no"
# = : Assign the variable name on the left to the value on the right

if [[ $SPEAK == "yes" ]];
then
  echo "Hello, World!"
else
  echo "sshhhh"
fi

```

The `[[ ]]` syntax here is a *bash test operator*, that returns `True` or `False.` The commands after `then` are run if the condition(s) inside the `[[]]` are `True`, and the commands after `else` are run if the condition(s) are `False`. 

Importantly, `else` statements are not always necessary if you don't want to perform any actions if the **conditional** (`if`) statement is `False`.

> Run the code block below to see an `if` statement without an `else` statement:

```{bash}

SPEAK="no"
# = : Assign the variable name on the left to the value on the right

if [[ $SPEAK == "yes" ]];
then
  echo "Hello, World!"
fi

```

There are some particularly useful conditional tests built in to bash that let you check file properties:

`-e filename` returns *True* if `filename` exists
`-s filename` returns *True* if `filename` exists and has a size > 0
`-d filename` returns *True* if `filename` is a directory

You can also check string properties:

`-z string` returns *True* if `string` is empty or null (0-length)
`-n string` returns *True* if `string` is not empty/null

> **Exercise**:
> Write an `if` statement that will check if a file exists and has data, printing (with `echo`) "file exists" if it exists and "file is missing" if it doesn't exist. Try changing the FILE variable to point to something else.

```{bash}

FILE="snp-densitsy.sh"
# = : Assign the variable name on the left to the value on the right

## Write an if statement to check if FILE exists and has data

## Write an if statement to check if FILE exists and has data

```

> **Exercise**:
> 1. In our `snp-density.sh` script, use an `if` statement to check if each file in the input file of bed files exists and only run our `bedtools` command on the file if it does. Note that you don't need to include an else statement if you don't want.
> 2. **BONUS**: Another user-friendly behavior for scripts is to tell the user how to run the program. A full list of **command line options** is warranted for some programs, but for a simply one like this only a simple message is required. To do this, write another `if` statement that checks if both of the **command line options** (`$1` and `$2`) have been provided and if not prints a user-friendly message about how to run the program, e.g:

```
echo "Usage: ./snp-density.sh <file of bed files to process> <vcf file>"
```

> For this you will have to check two conditions with an AND operator using the `&&` syntax, e.g. `if [[ test1 && test2 ]]`.

> 3. Run the code block below to test your script. If you did the bonus, remove either `bed_files` or `$VCF` from the command that calls your script to test the usage message:

```{bash}

VCF="data3/poeFor_NW_006799939.vcf"
# = : Assign the variable name on the left to the value on the right

echo "bed_doesnt_exist.bed" >> bed_files
# Add a non-existant file name to our list of bed files to check our error-handling if statement

./snp-density.sh bed_files $VCF
# ./ : Execute the following file as a shell script

```

# End of Day 3

Tomorrow, we'll learn how to run scripts like this using the job scheduling software **SLURM**, so we can submit them to run non-interactively on the cluster.

