Physalia Speciation Genomics Course

https://speciationgenomics.github.io/
Admin
This project provides tutorials and resources for the Physalia Speciation Genomics Course, covering...
Tokens:68,164
Snippets:919
Trust Score:9.7
Update:3 months ago
Show doc for...
Tokens:
Raw
### Setup Working Directory and Folders

Source: https://speciationgenomics.github.io/orthofinder

Creates a dedicated directory for the Orthofinder exercise and organizes subfolders for data and results. This is the initial setup step for the tutorial.

```bash
# We make a folder and move inside
cd ~
mkdir Orthofinder_exercise
cd Orthofinder_exercise

# We make two folders for us to work in
mkdir 01_data 02_orthofinder_run

```

--------------------------------

### Setup R Environment and Load Packages

Source: https://speciationgenomics.github.io/candidate_genes

Clears the R environment and loads the tidyverse package. This is the initial setup for the analysis.

```R
# clear environment
rm(list = ls())
# load packages
library(tidyverse)

```

--------------------------------

### Basic 'parallel' Syntax Example

Source: https://speciationgenomics.github.io/mapping_reference

Demonstrates the fundamental syntax of the 'parallel' command, showing how to read input from a file and use a placeholder for input values.

```bash
parallel 'echo {}' :::: inds
```

--------------------------------

### Create Directories

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `mkdir` to create new directories. This example creates a main directory and two subdirectories within the home directory.

```bash
cd
mkdir unix_test
mkdir unix_test/dir1 unix_test/dir2

```

--------------------------------

### Set up Environment for Genome Scan

Source: https://speciationgenomics.github.io/per_site_Fst

Sets up the working directory and an environmental variable for the input VCF file before starting genome scan analyses.

```bash
# move to home directory
cd ~
# make a genome scan directory
mkdir genome_scan
# make an environmental variable for the input vcf
VCF=~/vcf/chr20.vcf.gz

```

--------------------------------

### Start Screen Session for Indexing

Source: https://speciationgenomics.github.io/mapping_reference

Starts a new 'screen' session named 'index'. This allows background processing of long-running commands like genome indexing.

```bash
screen -S index

```

--------------------------------

### Pseudo-code For Loop Example

Source: https://speciationgenomics.github.io/more_advanced_R

Illustrates the concept of a for loop for iterating through a set of items and performing an action on each.

```pseudo-code
# a pseudo-code exmaple of a for loop
for(i in 1:10 tanks){
  stop in front of tank
  feed 20 g of fish food
  move to next tank
}

```

--------------------------------

### Get help for a command

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `man` or `--help` flags to access the manual pages or help information for a command.

```bash
man rm

```

```bash
rm --help

```

--------------------------------

### Install and Load Tidyverse Package

Source: https://speciationgenomics.github.io/data_visualisation_handling

Installs the tidyverse package if not already present and then loads it into the R session. This is a prerequisite for using many data manipulation and visualization functions.

```r
install.packages("tidyverse")
library(tidyverse)

```

--------------------------------

### Get FastQC Help

Source: https://speciationgenomics.github.io/readsExploration

Display the help message for the FastQC tool, which is used for assessing read quality. The output is piped to 'less' for easier viewing.

```bash
fastqc -h | less

```

--------------------------------

### Get Help for an R Function

Source: https://speciationgenomics.github.io/R_introduction

Access help documentation for R functions by prefixing the function name with a question mark. This opens a help dialogue with usage examples.

```r
?getwd

```

--------------------------------

### Unix Terminal Prompt Example

Source: https://speciationgenomics.github.io/getting_used_to_unix

This is an example of a typical Unix terminal prompt. The '$' symbol indicates that the system is ready to accept commands.

```bash
[mark@ourserver ~]$

```

--------------------------------

### Pseudo-code If-Else Statement Example

Source: https://speciationgenomics.github.io/more_advanced_R

Demonstrates a basic if-else structure to control code execution based on a condition.

```pseudo-code
# a pseudo-code exmaple of an if statement
if(cup != full){
  pour kettle
} else {
  drink tea
}

```

--------------------------------

### Install an R Package

Source: https://speciationgenomics.github.io/R_introduction

Installs the 'ggplot2' package from the CRAN repository. This function only needs to be run once per package.

```R
install.packages("ggplot2")
```

--------------------------------

### Skip Lines from Start of File with tail

Source: https://speciationgenomics.github.io/going_further_unix

Displays the content of a file starting from a specified line number. The -n+ flag followed by the line number indicates the starting point.

```bash
tail -n+3 my_file.txt
```

--------------------------------

### List Files Starting or Containing a Character

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `ls` with patterns to list files that start with a specific character or contain a specific character anywhere in their name.

```bash
ls c*

```

```bash
ls *c*

```

--------------------------------

### Navigate Using Absolute and Relative Paths

Source: https://speciationgenomics.github.io/getting_used_to_unix

Navigate to a specific directory using its absolute path (starting from root) or a relative path from the current location. The `$HOME` variable can substitute for `~`.

```bash
cd ~/unix_test/dir1

```

```bash
cd ~/unix_test/dir2

```

```bash
cd $HOME/unix_test/dir2

```

```bash
cd ../dir1

```

--------------------------------

### Load R Packages

Source: https://speciationgenomics.github.io/haplotypes

Clears the R environment and loads the 'rehh' and 'tidyverse' packages. Ensure these packages are installed before running.

```R
# clear environment
rm(list = ls())
# load packages
library(rehh)
library(tidyverse)

```

--------------------------------

### List Available Datasets in R

Source: https://speciationgenomics.github.io/R_introduction

Use the `data()` function to display a list of all datasets available in R. This is a good starting point for exploring built-in data.

```R
data()

```

--------------------------------

### Complete R Script for Data Visualization

Source: https://speciationgenomics.github.io/data_visualisation_handling

A comprehensive R script demonstrating environment setup, package loading, data filtering, plotting, and saving a plot to a PDF file. This serves as a template for standalone analysis scripts.

```r
### A simple R script to filter and plot the star wars data ###
rm(list = ls())

# load packages
library(tidyverse)

# load the star wars data
data(starwars)

# remove Jabba
starwars2 <- filter(starwars, name != "Jabba Desilijic Tiure")

# create a height/mass scatterplot
a <- ggplot(starwars2, aes(height, mass))
a <- a + geom_point(colour = "red")
a <- a + geom_smooth(method = lm, se = FALSE)
a <- a + theme_light() + ggtitle("Height vs Mass in Star Wars")
a

# write out the plot as a pdf
ggsave("./starwars_height_mass.pdf", plot = last_plot(), device = "pdf")
```

--------------------------------

### Run fastsimcoal2 with specific parameters

Source: https://speciationgenomics.github.io/fastsimcoal2

Execute fastsimcoal2 with a MAF, ignoring monomorphic sites, pooling low-SNP SFS entries, and specifying simulation and optimization cycles. Use this command to start a new simulation run.

```bash
PREFIX="early_geneflow"
mkdir fastsimcoal
cd fastsimcoal
mkdir early_geneflow
cd early_geneflow
cp /home/data/fastsimcoal2/early_geneflow/* ./
fastsimcoal2 -t ${PREFIX}.tpl -e ${PREFIX}.est -m -0 -C 10 -n 10000 -L 40 -s 0 -M

```

--------------------------------

### Count Header Lines in FASTQ File

Source: https://speciationgenomics.github.io/readsExploration

Count all lines starting with '@' in a gzipped FASTQ file using 'zgrep -c'. Note that this method can overcount due to '@' also appearing in quality score lines.

```bash
zgrep "@" $FILE -c

```

--------------------------------

### Set up VCF and working directory

Source: https://speciationgenomics.github.io/easysfs

Prepare the working directory and define the VCF file path for easySFS analysis. This involves creating a directory, changing into it, and copying the relevant VCF file.

```bash
# As always, we work in a separate directory for this analysis
mkdir SFS
cd SFS

# We will work with a subset of sites that are already filtered for low missing data proportions:
cp /home/data/vcf/cichlid_subset_filtered.vcf.gz ./ 
VCF="cichlid_subset_filtered.vcf.gz"

```

--------------------------------

### Create and Navigate to Phasing Directory

Source: https://speciationgenomics.github.io/phasing

Sets up the working directory for phasing analysis by creating a 'phasing' folder in the home directory and moving into it.

```bash
# move into home directory
cd ~
# make phasing directory
mkdir phasing
# move inside
cd phasing

```

--------------------------------

### Create a Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the 'mkdir' command to create a new directory. You can create multiple directories at once.

```bash
mkdir stuff1

```

```bash
mkdir stuff2 stuff3

```

--------------------------------

### Select Columns Starting with a Prefix using dplyr

Source: https://speciationgenomics.github.io/data_visualisation_handling

Selects columns from the 'starwars' dataset whose names start with 's' using dplyr.

```r
# choose only columns beginning with "s"
starwars %>% select(starts_with("s"))
```

--------------------------------

### Create Files

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the `touch` command to create empty files.

```bash
cd ~/unix_test/dir1
touch abc.txt abc.jpg xyz.txt xyz.jpg cat.txt car.txt

```

--------------------------------

### Change to Home Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use 'cd ~' to navigate to your home directory. This is a shorthand for the home directory.

```bash
cd ~

```

--------------------------------

### Set up working directory for LD decay analysis

Source: https://speciationgenomics.github.io/ld_decay

Navigate to your home directory and create a new directory for LD decay calculations. This ensures a clean working environment.

```bash
cd ~
mkdir ld_decay
cd ld_decay
```

--------------------------------

### Get Vector Length

Source: https://speciationgenomics.github.io/R_introduction

Determines the number of elements in a vector using the length() function.

```R
length(x)
```

```R
length(y)
```

--------------------------------

### Display fastp help information

Source: https://speciationgenomics.github.io/fastp

Check the available options and parameters for the fastp command.

```bash
fastp -h
```

--------------------------------

### Prepare ADMIXTOOLS Directory and Convert VCF

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

This snippet shows how to create a directory for ADMIXTOOLS analyses, copy a VCF file and a conversion script into it, and then execute the script to convert the VCF to Eigenstrat format.

```bash
# move to your home directory
cd ~
# make a directory for ADMIXTOOLS analyses
mkdir admixtools

# move into that new folder
cd admixtools

# copy the dogs dataset to the new folder
cp /home/data/sparrows/sparrows_chr8_unphased.vcf.gz ./

# copy the conversion script to your directory
cp /home/scripts/convertVCFtoEigenstrat.sh .

# Convert the vcf to the eigenstrat format required by AMDIXTOOLS
sh convertVCFtoEigenstrat.sh sparrows_chr8_unphased.vcf.gz

```

--------------------------------

### Create and List Text Files

Source: https://speciationgenomics.github.io/going_further_unix

Creates three text files in the current directory. Use `ls *.txt` to verify their creation.

```bash
cd ~/home/dir2
touch file1.txt file2.txt file3.txt
```

```bash
ls *.txt
```

--------------------------------

### Clear R Environment

Source: https://speciationgenomics.github.io/data_visualisation_handling

Start R scripts with this command to remove all objects and data from the environment, preventing potential conflicts.

```r
rm(list = ls())
```

--------------------------------

### Access Single Vector Element

Source: https://speciationgenomics.github.io/R_introduction

Retrieves a specific element from a vector using its index. R indices start at 1.

```R
x[5]
```

```R
y[5]
```

--------------------------------

### Navigate to Root, Home, or Current Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `cd` to change directories. `cd /` goes to the root, `cd ~` or `cd` goes to the home directory.

```bash
cd /
cd ~
cd

```

--------------------------------

### Display Sample Information

Source: https://speciationgenomics.github.io/phasing

Shows the first few lines of the '.samples' file generated by shapeit, which lists sample IDs and missing data proportions.

```bash
head sparrows_chr8_phased.samples

```

--------------------------------

### Load admixr Package in R

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

Loads the admixr package into the R session. Ensure the package is installed before running this command.

```r
library(admixr)

```

--------------------------------

### Copy Keyfile to Home Directory

Source: https://speciationgenomics.github.io/logging_on

Copies the provided keyfile to your home directory. Replace 'c123.pem' with your actual keyfile name and '123' with your user number.

```bash
cp c123.pem ~

```

--------------------------------

### Run FastQC on a Single FASTQ File

Source: https://speciationgenomics.github.io/readsExploration

Run FastQC on a specified gzipped FASTQ file to assess read quality. The file name is provided as an argument.

```bash
fastqc $FILE

```

--------------------------------

### Get Italian population names

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

Filters a list of populations to extract only those containing 'italian' in their name, preparing for batch processing.

```R
italians <- pops %>% grep("italian", ., value = T)

```

--------------------------------

### View First Read in FASTQ File

Source: https://speciationgenomics.github.io/readsExploration

This snippet shows how to view the first read (4 lines) from a gzipped FASTQ file. It requires creating a directory, copying the file, and then using 'zcat' piped to 'head'.

```bash
# First, we will make a folder to work in
mkdir fastqc
cd fastqc

# Now let's specify FILE as the name of the file containing the forward reads
FILE="wgs.R1.fastq.gz"
cp /home/data/fastq/$FILE ./

# Let's have a look at the first read:
zcat $FILE | head -4

```

--------------------------------

### Verify and Inspect Matrix Dimensions in R

Source: https://speciationgenomics.github.io/R_introduction

Use is.matrix to check if an object is a matrix and dim to get its dimensions (rows and columns).

```r
is.matrix(z)
dim(z)
```

--------------------------------

### Create a Basic Bash Script

Source: https://speciationgenomics.github.io/going_further_unix

This snippet shows the structure of a minimal Bash script, including the shebang line, a comment, and a command to print 'Hello world'.

```bash
#!/bin/sh

# a simple bash script
echo "Hello world"

exit

```

--------------------------------

### Split VCF by Sample List

Source: https://speciationgenomics.github.io/phasing

Utilize bcftools view with the -S flag to create new VCF files for each population. Use -O z for compressed output and -o for specifying the output file.

```bash
bcftools view -S house -O z -o house_chr8.vcf.gz $VCF
bcftools view -S bac -O z -o bac_chr8.vcf.gz $VCF
```

--------------------------------

### Navigate Up One or More Directories

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `cd ..` to move up one directory level. Multiple `..` can be chained or combined.

```bash
cd ..
cd ..

```

```bash
cd ../..

```

--------------------------------

### Compute f4 Ratio with admixr

Source: https://speciationgenomics.github.io/_pages/ADMIXTOOLS

Example of computing an f4 ratio test using the admixr package. Specify the taxa and the data source.

```r
test4 <- f4ratio(X = "Basenji", A = "ChineseWolf", B = "IsraeliWolf", C = "Dingo", O = "GoldenJackal", data = snps)
```

--------------------------------

### Get Matrix Row and Column Counts in R

Source: https://speciationgenomics.github.io/R_introduction

Retrieve the number of rows and columns of a matrix using nrow and ncol functions, respectively.

```r
nrow(z)
col(z)
```

--------------------------------

### Perform Block-Bootstrapping and SFS Calculation

Source: https://speciationgenomics.github.io/fastsimcoal2

This loop performs 50 iterations of block-bootstrapping. In each iteration, it creates a new directory, constructs a bootstrapped VCF file by randomly concatenating blocks of sites, compresses it, and then computes the Site Frequency Spectrum (SFS) using easySFS.py. It also copies the observed SFS file for comparison.

```bash
# Generate 50 files each with randomly concatenated blocks and compute the SFS for each:
for i in {1..50}
do
  # Make a new folder for each bootstrapping iteration:
  mkdir bs$i
  cd bs$i

  # Add the header to our new bootstrapped vcf file
  cat ../header > $PREFIX.bs.$i.vcf
  # Randomly add 100 blocks
  for r in {1..100}
  do
    cat `shuf -n1 -e ../$PREFIX.sites.*` >> ${PREFIX}.bs.$i.vcf
  done
  # Compress the vcf file again
  gzip ${PREFIX}.bs.$i.vcf

  # Make an SFS from the new bootstrapped file
  easySFS.py -i ${PREFIX}.bs.$i.vcf.gz -p pop_file -a -f --proj 8,8

  # Copy the observed SFS file into this folder renaming it to match the .tpl prefix
  cp ../${PREFIX}_jointDAFpop1_0.obs  ${PREFIX}.bs.${i}_jointDAFpop1_0.obs

  # Say that it is finished with iteration $i
  echo bs$i" ready"

  cd ..
done

```

--------------------------------

### Get Current Working Directory in R

Source: https://speciationgenomics.github.io/R_introduction

Use the getwd() function to determine the directory R is currently operating from. This is important for managing data files.

```r
getwd()

```

--------------------------------

### Change to Home Directory using Environmental Variable

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use 'cd $HOME' to navigate to your home directory using the HOME environmental variable.

```bash
cd $HOME

```

--------------------------------

### Create and Navigate to PLINK Directory

Source: https://speciationgenomics.github.io/pca

Sets up a dedicated directory for PLINK operations and navigates into it. This is a preparatory step for subsequent PLINK commands.

```bash
# move to your home directory
cd ~
# make a plink directory
mkdir plink
# move into it
cd plink

```

--------------------------------

### Clone Unix Exercises Repository

Source: https://speciationgenomics.github.io/going_further_unix

Clones a GitHub repository containing Unix exercises and data files into the home directory. This is useful for obtaining practice files and ensuring you have the necessary data for exercises.

```bash
# move into your home directory
cd ~
# clone the appropriate repository
git clone https://github.com/speciationgenomics/unix_exercises.git
```

--------------------------------

### Create Numeric Vectors

Source: https://speciationgenomics.github.io/R_introduction

Creates two numeric vectors, 'x' using a sequence of integers and 'y' using the seq() function with specified start, stop, and increment values.

```R
x <- 1:10
y <- seq(from = 10, to = 100, by = 10)
```

--------------------------------

### Run fastsimcoal2 with Bootstrapped SFS

Source: https://speciationgenomics.github.io/fastsimcoal2

This script automates running fastsimcoal2 for parameter estimation on bootstrapped SFS. It iterates through each bootstrap replicate, runs fastsimcoal2 100 times, and then selects the best run for each replicate. Ensure you are in the correct directory and have the necessary input files.

```bash
for bs in {1..50}
do
  cd bs$bs
  # Run fastsimcoal 100 times:
  for i in {1..100}
  do
    mkdir run$i
    cd run$i
    cp ${PREFIX}.bs.$bs.* ./
    fastsimcoal2 -t ${PREFIX}.bs.$bs.tpl -e ${PREFIX}.bs.$bs.est -m -0 -C 10 -n 10000 -L 40 -s0 -M -q
    cd ..
  done
  # Find the best run:
  fsc-selectbestrun.sh

  cd ..
done

```

--------------------------------

### Display Line Range with sed

Source: https://speciationgenomics.github.io/going_further_unix

Prints a range of lines from a file. The -n flag suppresses default output, and 'N,Mp' specifies the start and end line numbers.

```bash
sed -n 3,5p ishmael.txt
```

--------------------------------

### Navigate to a Subdirectory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Change directories sequentially to reach a nested directory.

```bash
cd unix_test
cd dir1

```

--------------------------------

### List Files Excluding Patterns with Character Sets

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use character set negation `[^...]` with `ls` for more flexible pattern matching, excluding files that start with specified characters.

```bash
ls [^x]*

```

```bash
ls [^xa]*

```

--------------------------------

### List Directory Contents

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the 'ls' command to list the files and directories in the current directory.

```bash
ls

```

--------------------------------

### Count Sequences in FASTQ File

Source: https://speciationgenomics.github.io/readsExploration

Count the number of sequences in a gzipped FASTQ file by counting lines starting with '@J00' using 'zgrep -c'. This is a more reliable method than counting all '@' symbols.

```bash
zgrep "@J00" $FILE -c

```

--------------------------------

### Modify Column Values with awk

Source: https://speciationgenomics.github.io/going_further_unix

Perform arithmetic operations on column values using awk. This example adds 1 to each value in the third column and prints it alongside the original column.

```bash
awk '{print $3,"\t"$3+1}' iris_data.tsv | head

```

--------------------------------

### Prepare Data for Block-Bootstrapping

Source: https://speciationgenomics.github.io/fastsimcoal2

This script prepares genomic data for block-bootstrapping by extracting all sites and the header from a VCF file, then splitting the sites into smaller files. It's a prerequisite for the bootstrapping process.

```bash
# Get all lines with genomic data
zgrep -v "^#" $PREFIX.vcf > $PREFIX.allSites

# Get the header
zgrep "^#" $PREFIX.vcf > header

# get 100 files with 4338 sites each (number 101 removed due to only 90 sites)
split -l 4338 $PREFIX.allSites $PREFIX.sites.

```

--------------------------------

### Run FastQC on Reverse Reads with Output Directory

Source: https://speciationgenomics.github.io/readsExploration

Run FastQC on the reverse reads FASTQ file, specifying an output directory using the '-o' flag. This helps organize the generated quality reports.

```bash
# Reverse reads
FILE=wgs.R2.fastq.gz
fastqc -o ./ /home/data/fastq/$FILE

```

--------------------------------

### Calculate D statistic for single population set

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

Example of calculating the D statistic for a specific set of four taxa using the `d()` function in admixr. Ensure the 'snps' data and population labels are correctly defined.

```r
d1 <- d(W = "house-european_house-uk", X = "bactrianus-bactrianus-kazahkstan", Y =  "spanish-spanish-cv", Z = "tree", data = snps)
```

--------------------------------

### Select Best fastsimcoal2 Run

Source: https://speciationgenomics.github.io/fastsimcoal2

Use a helper script to automatically identify the best performing fastsimcoal2 run based on likelihood and copy its files to a 'bestrun' directory.

```bash
fsc-selectbestrun.sh

```

--------------------------------

### Combine Pattern and Column Filters

Source: https://speciationgenomics.github.io/awk

Combine multiple conditions using logical operators like '&&' (AND) to refine line selection. This example selects lines containing 'Sa' and having more than 0.01 missing data.

```awk
awk '/Sa/ && $5>0.01' cichlids.imiss
```

--------------------------------

### Download HTML Files from Server

Source: https://speciationgenomics.github.io/readsExploration

Download all HTML files from the 'fastqc' directory on a remote server to the local current directory using the 'scp' command. Requires specifying the private key, user, IP address, and remote path.

```bash
scp -i c1.pem user1@<IP address>:/home/user1/fastqc/*html ./

```

--------------------------------

### Filter Lines by Column Value

Source: https://speciationgenomics.github.io/awk

Select lines where a specific column's value meets a condition. $N refers to the Nth column. This example filters for missing data proportion greater than 0.01.

```awk
awk '$5>0.01' cichlids.imiss
```

--------------------------------

### Load an R Package

Source: https://speciationgenomics.github.io/R_introduction

Loads the 'ggplot2' package into the current R session, making its functions and data available for use. This must be run each time you start a new R session and want to use the package.

```R
library("ggplot2")
```

--------------------------------

### Print Working Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the 'pwd' command to display the full path of the current working directory.

```bash
pwd

```

--------------------------------

### List Files in Current or Parent Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `ls .` to list files in the current directory and `ls ..` to list files in the parent directory.

```bash
ls .
ls ..

```

--------------------------------

### Plotting Fst using ggplot2 in R

Source: https://speciationgenomics.github.io/sliding_windows

Reads sliding window statistics from a CSV file and plots the Fst values between two specified populations using ggplot2. Ensure the 'ggplot2' library is installed and loaded.

```r
rm(list = ls())

# Prepare input files:
# Read in the file with sliding window estimates of FST, pi and dxy
windowStats<-read.csv("genomeScans/Pundamilia_kivu_div_stats.csv",header=T)

# Read in the fd estimates of 20 kb windows for NyerMak into NyerPyt (P1=PundPyt, P2=NyerPyt, P3=NyerMak, outgroup=Kivu cichlid)
fd<-read.csv("genomeScans/Pundamilia_ABBABABA.csv",header=T,na.strings = "NaN")

# Let's have a look at the FST and fd datasets
head(windowStats)
head(fd)

# Let's plot FST between the two younger species
require(ggplot2)
ggplot(windowStats,aes(mid,Fst_PundPyt_NyerPyt))+geom_point()
# Annoyingly, it plots all the FST values on top of each other

```

--------------------------------

### Create and Navigate VCF Directory

Source: https://speciationgenomics.github.io/variant_calling

Basic shell commands to create a directory for VCF files and move a gzipped VCF file into it.

```bash
mkdir vcf
mv cichlid.vcf.gz ./vcf
cd vcf
```

--------------------------------

### Convert VCF to Geno.gz Format

Source: https://speciationgenomics.github.io/sliding_windows

Converts a VCF file into the geno.gz format required by Simon Martin's population genetics scripts. Ensure you have Python 2 installed and adjust the shebang line if necessary.

```bash
cd ~
mkdir genome_scans
cd genome_scans

# Copy the vcf file and the information file and a small subset file
cp /home/data/genomeScans/* ./

# Convert the vcf file to geno.gz which is the format that Simons script requires
parseVCF.py -i subset.vcf.gz | bgzip > subset.geno.gz

```

--------------------------------

### SLiM: Simulate Stable Population

Source: https://speciationgenomics.github.io/slim

Initializes a neutral simulation with a constant population size of 1000 diploids over 5000 generations. Outputs VCF files for 20 sampled individuals.

```slim
// Constant popiulation size of 1000 with no selection
initialize()
{
	// set the overall mutation rate
	initializeMutationRate(1e-5);

	// m1 mutation type: neutral
	initializeMutationType("m1", 0.5, "f", 0.0);

	// g1 genomic element type: uses m1 for all mutations
	initializeGenomicElementType("g1", m1, 1.0);

	// Chromosome of length 10 kb, homogenous (all of type g1)
	initializeGenomicElement(g1, 0, 9999);

	// uniform recombination along the chromosome
	initializeRecombinationRate(1e-8);
}

// create a population of 1000 individuals
1 early()
{
	sim.addSubpop("p1", 1000);
}

// run to generation 5000
5000 late() {

// randomly subsample 20 individuals for output
allIndividuals = sim.subpopulations.individuals;
 sampledIndividuals = sample(allIndividuals, 20);
 sampledIndividuals.genomes.outputVCF();

}


```

--------------------------------

### Calculate D statistic for multiple population sets

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

Example of calculating the D statistic for multiple populations by passing a vector of population labels to the `W` argument in the `d()` function. This allows for efficient testing across many population combinations.

```r
d2 <- d(W = pops, X = "bactrianus-bactrianus-kazahkstan", Y = "spanish-spanish-cv", Z = "tree", data = snps)
```

--------------------------------

### Create Empty Files

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the 'touch' command to create empty files. It is also used to update file access times.

```bash
touch file1 file2

```

--------------------------------

### Print Entire Line Based on Multiple Conditions

Source: https://speciationgenomics.github.io/awk

Print the entire line ($0) when multiple conditions are met. This example prints lines where the first column contains 'Sa' and the fifth column exceeds 0.005.

```awk
awk '{if($1~/Sa/ && $5>0.005) print $0}' cichlids.imiss
```

--------------------------------

### Change Directory

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use the 'cd' command followed by the directory name to move into that directory.

```bash
cd stuff1

```

--------------------------------

### Run FastQC on Multiple FASTQ Files using a Loop

Source: https://speciationgenomics.github.io/readsExploration

Use a for loop to efficiently run FastQC on multiple FASTQ files, including RAD and Novaseq datasets. The '-o' flag specifies the output directory.

```bash
# Let's also have a look at some RAD data  and a Novaseq whole genome dataset
# Because we do not want to type the fastqc command many times, we will use a for loop
for PREFIX in RAD1 RAD2 Novaseq.R1 Novaseq.R2
do
fastqc -o ./ /home/data/fastq/$PREFIX.fastq.gz
done

```

--------------------------------

### Edit File with Nano

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use 'nano <filename>' to open a file in the nano text editor for editing. Save and exit using on-screen prompts.

```bash
nano file1

```

--------------------------------

### Extract Candidate Gene Attributes

Source: https://speciationgenomics.github.io/candidate_genes

Selects relevant columns (chromosome, start, end, attribute, distance to hit) for the candidate genes and extracts the 'attribute' column to display gene identifiers. This provides the final list of candidate genes.

```R
# what are these genes?
gene_hits <- gene_hits %>% select(chr, start,end, attribute,hit_dist)
# separate out the attribute column
gene_hits %>% pull(attribute)

```

--------------------------------

### Trigger RStudio Syntax Checking

Source: https://speciationgenomics.github.io/more_advanced_R

Examples demonstrating RStudio's code syntax checking capabilities. These snippets are designed to trigger visual cues in the RStudio script viewer, highlighting potential errors like missing arguments or mismatched brackets.

```R
# Try these to trigger RStudio syntax checking
plot(x, )
result <- c(1,2,3))
```

--------------------------------

### View FASTA Index

Source: https://speciationgenomics.github.io/variant_calling

Displays the first few lines of a FASTA index file created by samtools faidx. Useful for understanding the index structure.

```bash
head P_nyererei_v2.fasta.fai

```

--------------------------------

### Calculate f4 ratio for admixture proportion

Source: https://speciationgenomics.github.io/ADMIXTOOLS_admixr

Computes the admixture proportion (alpha) using an f4 ratio test, which requires five taxa: A, B, X, C, and an outgroup O. This specific example estimates Spanish sparrow ancestry in Italian sparrows.

```R
f4r1 <- f4ratio(X = "italian-italian_mainland-crotone", A =  "spanish-spanish-cv", B = "spanish-spanish-kazahkstan", C = "house-european_house-uk", O = "tree", data = snps)

```

```R
f4r1

```

--------------------------------

### Query VCF Header for INFO and FORMAT Descriptions

Source: https://speciationgenomics.github.io/variant_calling

Search the VCF header to understand the meaning of INFO and FORMAT fields.

```bash
bcftools view -h cichlid.vcf.gz | grep "INFO"
bcftools view -h cichlid.vcf.gz | grep "##FORMAT"
```

--------------------------------

### Create Directory for vcftools Output

Source: https://speciationgenomics.github.io/filtering_vcfs

Creates a directory to store the output files generated by vcftools.

```bash
mkdir ~/vcftools

```

--------------------------------

### Display Head of a File

Source: https://speciationgenomics.github.io/going_further_unix

Displays the first 10 lines of a file by default. Useful for quickly inspecting the beginning of large files.

```bash
head udhr.txt
```

--------------------------------

### Define Input VCF and Output Prefix Variables

Source: https://speciationgenomics.github.io/filtering_vcfs

Sets variables for the input gzipped VCF file and the output file prefix to simplify commands.

```bash
SUBSET_VCF=~/vcf/cichlid_subset.vcf.gz
OUT=~/vcftools/cichlid_subset

```

--------------------------------

### Visualize SFS Data with R Scripts

Source: https://speciationgenomics.github.io/fastsimcoal2

Utilize R scripts to visualize the simulated Site Frequency Spectrum (SFS) against observed data and to plot the best model fit. Requires downloading SFStools.r and plotModel.r.

```bash
SFStools.r -t print2D -i early_geneflow
plotModel.r -p early_geneflow -l NyerMak,PundMak

```

--------------------------------

### Shell Script for Multiple SLiM Iterations and Summary Statistics

Source: https://speciationgenomics.github.io/slim

This script runs 100 iterations of a SLiM simulation for a stable population size. For each iteration, it generates a VCF file, calculates FST, nucleotide diversity (pi), and Tajima's D using vcftools, and cleans up intermediate files. Finally, it collates the summary statistics from all iterations into single files.

```bash
for (( i = 1; i < 101; i++ )); do

  echo "running simulation ${i}"

  ### 1: Run the simulation
  slim neutral_stable.slim > stable_${i}.out # run the simulation
  tail -n+15 stable_${i}.out | bgzip -c > stable_${i}.vcf.gz # zip and tidy up the vcf output
  rm stable_${i}.out # remove a temporary file (raw SLIM output)

  ### 2: Calculate summary statistics

  ## Calculate Fst from the VCF
  vcftools --gzvcf stable_${i}.vcf.gz --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt --out stable_fst_${i} # produces a log file
  grep "weighted" stable_fst_${i}.log | awk '{print $7}' > fst_stable_${i}.txt # Extract the fst value from the vcftools log file
  rm stable_fst_${i}.weir.fst # remove a temporary file (genome-wide FST; we don't need this)
  rm stable_fst_${i}.log # remove a temporary file (the vcftools log file)


  ## Calculate nucleotide diversity from the VCF
  vcftools --gzvcf stable_${i}.vcf.gz --site-pi --out stable_pi_${i}
  awk 'NR > 1 {sum+=$3} END {print sum / (NR - 1)}' stable_pi_${i}.sites.pi > stable_pi_mean_${i}.txt  # calculate the mean
  rm stable_pi_${i}.log # remove intermediate logfile
  rm stable_pi_${i}.sites.pi # Remove the intermediate file

  ## Calculate Tajima's D from the VCF
  vcftools --gzvcf stable_${i}.vcf.gz --out stable_${i} --TajimaD 10000
  rm stable_${i}.log # remove the intermediate logfile
  awk 'NR==2' stable_${i}.Tajima.D | cut -f 4 > stable_tajD_${i}.txt # extract out the value of D
  rm stable_${i}.Tajima.D # remove the intermediate file

  rm stable_${i}.vcf.gz # remove a temporary file (vcf; we don't need this any more)

done


### Collate Fst data
cat fst_stable_*.txt >> all_fst_stable.txt # Put all 100 fst values into one text output
rm fst_stable_*.txt # remove the individual fst text files

# As above for pi and Tajima's D
cat stable_pi_mean_*.txt >> all_pi_stable.txt
rm stable_pi_mean_*.txt

cat stable_tajD_*.txt >> all_tajD_stable.txt
rm stable_tajD_*.txt

```

--------------------------------

### Log in to the Cluster via SSH

Source: https://speciationgenomics.github.io/logging_on

Connects to the cluster using SSH with the specified keyfile. Replace '~/c123.pem' with the path to your keyfile and 'user123@54.245.175.86' with your username and the cluster's IP address.

```bash
ssh -i "~/c123.pem" user123@54.245.175.86

```

--------------------------------

### Create population file

Source: https://speciationgenomics.github.io/easysfs

Generate a population file required by easySFS using bcftools, grep, and awk. This file maps sample IDs to population labels based on the VCF.

```bash
# make pops
bcftools query -l $VCF | grep "Mak" | awk '{split($0,a,"."); print $1,a[2]}' > pop_file

```

--------------------------------

### Copy Pre-built Index Files

Source: https://speciationgenomics.github.io/mapping_reference

Copies pre-built bwa-mem2 index files for the reference genome. This is an alternative to running the indexing command yourself, saving time.

```bash
cp /home/data/reference/P_nyererei_v2.fasta.* .

```

--------------------------------

### Upload a File to the Cluster using SCP

Source: https://speciationgenomics.github.io/logging_on

Uploads a local file ('test_file') to the home directory of the cluster. Ensure you have the correct keyfile path, username, and cluster IP address.

```bash
# make a file
touch test_file
# upload a file from your computer to the cluster
scp -i "~/c123.pem" test_file user123@54.245.175.86:~

```

--------------------------------

### Prepare Sequence Data with Species Prefixes

Source: https://speciationgenomics.github.io/orthofinder

Copies input data files into the data directory and modifies protein FASTA headers by adding species prefixes. This ensures clear identification of gene origins in downstream analysis.

```bash
cd ~/Orthofinder_exercise/01_data
cp /home/data/orthofinder/* .
# Take a minute or two to get acquainted with the data.
# Specifically - what do the .faa files contain?
# What does the general feature format file (gff) contain?

# Now, for the 'cleaning'.
# Orthofinder will give you tons of results. When you run it with more than two species it often gets confusing - which gene comes from species A? Is this gene from species C or D?
# I have a solution. Add the species name to each gene :)

# Run this code to extract all house sparrow genes and all chichen genes separately:
sed -i "/^>/ s/>/>House_sparrow_/" house_sparrow_proteins_reduced.faa
sed -i "/^>/ s/>/>chicken_/" Gallus_gallus.GRCg6a.pep.reduced.faa

# Check the files using less to see what output files we have now.

```

--------------------------------

### Run Treemix with Varying Migration Edges

Source: https://speciationgenomics.github.io/Treemix

Executes Treemix for a range of migration edges (0 to 5), specifying the input frequency file, output prefix, outgroup, and flags for bootstrapping and no small sample size correction. Each run is logged.

```bash
for i in {0..5}
do
 treemix -i $FILE.treemix.frq.gz -m $i -o $FILE.$i -root GoldenJackal -bootstrap -k 500 -noss > treemix_${i}_log &
done
```

--------------------------------

### Copy a file and rename it

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `cp` to copy a file to the current directory and give it a new name.

```bash
cd ../stuff2
cp file1 file1_copy

```

--------------------------------

### Estimate optimal SFS projections

Source: https://speciationgenomics.github.io/easysfs

Use easySFS.py with the --preview flag to estimate the best projection values for downsampling populations. This helps maximize the number of sites available for SFS calculation by accounting for missing data.

```bash
# estimate projections
easySFS.py -i $VCF -p pop_file -a -f --preview

```

--------------------------------

### Set Permissions for Keyfile

Source: https://speciationgenomics.github.io/logging_on

Adjusts permissions for the keyfile to ensure it is readable and executable by the owner. Use this if you encounter permission errors.

```bash
chmod go= c123.pem
chmod -R u+x c123.pem

```

--------------------------------

### Run ADMIXTURE in a Loop (K=3 to K=5)

Source: https://speciationgenomics.github.io/ADMIXTURE

This loop iterates from K=3 to K=5, running ADMIXTURE with cross-validation for each K value. The output for each run, including log information, is saved to separate files named 'log3.out', 'log4.out', and 'log5.out'.

```bash
for i in {3..5}
do
 admixture --cv $FILE.bed $i > log${i}.out
done
```

--------------------------------

### Extract Sample Names from VCF

Source: https://speciationgenomics.github.io/phasing

Use bcftools query to list all sample names in a VCF file and save them to a file. This is the first step to splitting populations.

```bash
# set a new variable
VCF=~/phasing/sparrows_chr8_phased.vcf.gz
# look at the sample names
bcftools query -l $VCF
# extract sample names
bcftools query -l $VCF > samples
```

--------------------------------

### SLiM: Simulate Population Expansion

Source: https://speciationgenomics.github.io/slim

Simulates a population that expands instantaneously from 1000 to 7500 individuals at generation 4000. This neutral simulation runs for 5000 generations and outputs VCF files for 20 sampled individuals.

```slim
// Instantaneous population expansion from 1000 to 7500 at generation 4000 with no selection
initialize()
{
	// set the overall mutation rate
	initializeMutationRate(1e-5);

	// m1 mutation type: neutral
	initializeMutationType("m1", 0.5, "f", 0.0);

	// g1 genomic element type: uses m1 for all mutations
	initializeGenomicElementType("g1", m1, 1.0);

	// Chromosome of length 10 kb, homogenous (all of type g1)
	initializeGenomicElement(g1, 0, 9999);

	// uniform recombination along the chromosome
	initializeRecombinationRate(1e-8);
}

// create a population of 1000 individuals
1 early()
{
	sim.addSubpop("p1", 1000);
}
// at generation 4000 increase the population size to 7500 immediately

4000 early() { p1.setSubpopulationSize(7500); }

// run to generation 5000
5000 late() {

// randomly subsample 20 individuals for output
allIndividuals = sim.subpopulations.individuals;
 sampledIndividuals = sample(allIndividuals, 20);
 sampledIndividuals.genomes.outputVCF();

}


```

--------------------------------

### View BAM File Content

Source: https://speciationgenomics.github.io/mapping_reference

Displays the beginning of a BAM file's content. This is useful for inspecting the alignment data. Note that direct 'head' command on a BAM file will not work due to its binary nature.

```bash
samtools view 10558.PunPundMak.bam | head
```

--------------------------------

### View First Few Alignments in SAM File

Source: https://speciationgenomics.github.io/mapping_reference

Displays the first few alignment records from a SAM file using samtools view. This provides a glimpse into the mapped read data.

```bash
samtools view 10558.PunPundMak.sam | head
```

--------------------------------

### List directory contents

Source: https://speciationgenomics.github.io/getting_used_to_unix

Use `ls` with different path specifiers to view the contents of the current directory, parent directory, home directory, and a specific subdirectory.

```bash
ls ./

```

```bash
ls ../

```

```bash
ls $HOME

```

```bash
ls $HOME/stuff2/

```

--------------------------------

### Run an Interactive Bash Script with Arguments

Source: https://speciationgenomics.github.io/going_further_unix

This command executes the 'name_script.sh' Bash script, passing 'Mark' as the first argument and 'Milo' as the second.

```bash
sh name_script.sh Mark Milo

```

--------------------------------

### View SAM File Header and First Alignments

Source: https://speciationgenomics.github.io/mapping_reference

Displays the header and the first few lines of a SAM alignment file using samtools view with the -h flag. This helps in understanding the SAM file structure.

```bash
cd align
samtools view -h 10558.PunPundMak.sam | head
```

--------------------------------

### Full Alignment and Sorting For Loop

Source: https://speciationgenomics.github.io/mapping_reference

Executes the BWA-MEM2 alignment and Samtools sorting for each individual within the loop. Includes echoing progress.

```bash
for IND in ${INDS[@]};
do
	# declare variables
	FORWARD=/home/data/wgs_raw/${IND}.R1.fastq.gz
	REVERSE=/home/data/wgs_raw/${IND}.R2.fastq.gz
	OUTPUT=~/align/${IND}_sort.bam

	# then align and sort
	echo "Aligning $IND with bwa"
	bwa-mem2 mem -t 4 $REF $FORWARD \
	$REVERSE | samtools view -b | \
	samtools sort -T ${IND} > $OUTPUT
done

```

--------------------------------

### Create and Inspect a Factor in R

Source: https://speciationgenomics.github.io/R_introduction

Demonstrates creating a character vector and converting it into a factor using `as.factor()`. The `levels()` function is then used to view the distinct categories within the factor.

```r
# create a character vector
myFactor <- c("male", "female", "male", "female", "female", "male")
# turn it into a factor
myFactor <- as.factor(myFactor)
# view the available levels
levels(myFactor)
```

--------------------------------

### Display Specific Number of Lines with Head

Source: https://speciationgenomics.github.io/going_further_unix

Displays a specified number of lines from the beginning of a file. Use the `-n` flag followed by the desired number of lines.

```bash
head -20 udhr.txt
```

--------------------------------

### Create Alignment Directory

Source: https://speciationgenomics.github.io/mapping_reference

Creates a new directory named 'align' in the home directory for storing alignment data. This is a preparatory step before running alignment commands.

```bash
cd ~
mkdir align
cd align
```

--------------------------------

### Download a File from the Cluster using SCP

Source: https://speciationgenomics.github.io/logging_on

Downloads a file ('test_file') from the cluster to your current local directory. Specify the cluster's file path and your local destination.

```bash
# download a file from the cluster to your computer
scp -i "~/c123.pem" user123@54.245.175.86:~/test_file ./

```

--------------------------------

### Shell Script for Running SLiM Simulations and Processing Output

Source: https://speciationgenomics.github.io/slim

This shell script automates the process of running SLiM simulations, extracting VCF data, and calculating summary statistics (FST, pi, Tajima's D) using vcftools. It is designed to run multiple iterations of a stable population model and collate the results.

```bash
# Make and go into a working directory called SLiM
cd ~
mkdir SLiM
cd SLiM/

# Get the shell scripts and slim files
cp /home/scripts/SLiM/*.sh .
cp /home/scripts/SLiM/*.slim .

# Let's check the version of slim
slim -version

# Have a look at the slim file
cat neutral_stable.slim

# Running slim is very easy, just specify the slim file
slim neutral_stable.slim > stable.out

# Familiarise yourself with the output file
head stable.out

# The last lines contain the genotype information, make a vcf file of those lines
tail -n+15 stable.out | bgzip -c > stable.vcf.gz

rm stable.out

# Compute FST with vcftools
vcftools --gzvcf stable.vcf.gz --weir-fst-pop pop1.txt --weir-fst-pop pop2.txt --out stable_fst

# have a look at the pop text files that are used to assign individuals to populations
cat pop1.txt
cat pop2.txt

# compute pi (nucleotide diversity) with vcftools
vcftools --gzvcf stable.vcf.gz --site-pi --out stable_pi

# calculate the mean and remove the intermediate file
awk 'NR > 1 {sum+=$3} END {print sum / (NR - 1)}' stable_pi.sites.pi > stable_pi_mean.txt
rm stable_pi.sites.pi

# Compute Tajima's D with vcftools
vcftools --gzvcf stable.vcf.gz --out stable --TajimaD 10000

cat stable.Tajima.D | column -t

awk 'NR==2' stable.Tajima.D | cut -f 4 > stable_tajD.txt
rm stable.Tajima.D

```
Tokens:
Raw